]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add routines for multiplication of ScaLAPACKMatrices + tests
authorBenjamin Brands <benjamin.brands@fau.de>
Tue, 6 Feb 2018 15:47:04 +0000 (16:47 +0100)
committerBenjamin Brands <benjamin.brands@fau.de>
Wed, 14 Feb 2018 20:10:16 +0000 (21:10 +0100)
include/deal.II/lac/scalapack.h
include/deal.II/lac/scalapack.templates.h
source/lac/scalapack.cc
tests/scalapack/scalapack_12.cc [new file with mode: 0644]
tests/scalapack/scalapack_12.mpirun=1.output [new file with mode: 0644]
tests/scalapack/scalapack_12.mpirun=11.output [new file with mode: 0644]
tests/scalapack/scalapack_12.mpirun=9.output [new file with mode: 0644]

index 934b2fea7bca7c7845d7559aee7edf0d45786c24..788355e4a11352d7e3913f0616a03c4d89d04f4a 100644 (file)
@@ -246,6 +246,108 @@ public:
   void Tadd(const NumberType b,
             const ScaLAPACKMatrix<NumberType> &B);
 
+  /**
+   * Matrix-matrix-multiplication:
+   *
+   * The operations based on the input parameters and the alignment conditions are summarized in the following table:
+   *
+   * | transpose_A | transpose_B |                   Block Sizes                 |                     Operation                    |
+   * | :---------: | :---------: | :-------------------------------------------: | :----------------------------------------------: |
+   * | false       |   false     | MB_A = MB_C <br> NB_A = MB_B <br> NB_B = NB_C |             <i>C = a A * B + b C</i>             |
+   * | false       |   true      | MB_A = MB_C <br> NB_A = NB_B <br> MB_B = NB_C |      <i>C = a A * B<sup>T</sup> + b C</i>        |
+   * | true        |   false     | MB_A = MB_B <br> NB_A = MB_C <br> NB_B = NB_C |        <i>C = a A<sup>T</sup> * B + b C</i>      |
+   * | true        |   true      | MB_A = NB_B <br> NB_A = MB_C <br> MB_B = NB_C | <i>C = a A<sup>T</sup> * B<sup>T</sup> + b C</i> |
+   *
+   * It is assumed that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+   * <tt>C</tt> already has the right size.
+   *
+   * The matrices <tt>A</tt>, <tt>B</tt> and <tt>C</tt> must have the same process grid.
+   */
+  void mult(ScaLAPACKMatrix<NumberType> &C,
+            const ScaLAPACKMatrix<NumberType> &B,
+            const NumberType a=1,
+            const NumberType b=0,
+            const bool transpose_A=false,
+            const bool transpose_B=false) const;
+
+  /**
+   * Matrix-matrix-multiplication.
+   *
+   * The optional parameter <tt>adding</tt> determines, whether the result is
+   * stored in <tt>C</tt> or added to <tt>C</tt>.
+   *
+   * if (adding) <i>C += A*B</i>
+   *
+   * else <i>C = A*B</i>
+   *
+   * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+   * <tt>C</tt> already has the right size.
+   *
+   * Following alignment conditions have to be fulfilled: MB_A = MB_C, NB_A = MB_B and NB_B = NB_C
+   */
+  void mmult(ScaLAPACKMatrix<NumberType> &C,
+             const ScaLAPACKMatrix<NumberType> &B,
+             const bool adding=false) const;
+
+  /**
+   * Matrix-matrix-multiplication using transpose of <tt>this</tt>.
+   *
+   * The optional parameter <tt>adding</tt> determines, whether the result is
+   * stored in <tt>C</tt> or added to <tt>C</tt>.
+   *
+   * if (adding) <i>C += A<sup>T</sup>*B</i>
+   *
+   * else <i>C = A<sup>T</sup>*B</i>
+   *
+   * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+   * <tt>C</tt> already has the right size.
+   *
+   * Following alignment conditions have to be fulfilled: MB_A = MB_B, NB_A = MB_C and NB_B = NB_C
+   */
+  void Tmmult (ScaLAPACKMatrix<NumberType> &C,
+               const ScaLAPACKMatrix<NumberType> &B,
+               const bool adding=false) const;
+
+  /**
+   * Matrix-matrix-multiplication using transpose of <tt>B</tt>.
+   *
+   * The optional parameter <tt>adding</tt> determines, whether the result is
+   * stored in <tt>C</tt> or added to <tt>C</tt>.
+   *
+   * if (adding) <i>C += A*B<sup>T</sup></i>
+   *
+   * else <i>C = A*B<sup>T</sup></i>
+   *
+   * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+   * <tt>C</tt> already has the right size.
+   *
+   * Following alignment conditions have to be fulfilled: MB_A = MB_C, NB_A = NB_B and MB_B = NB_C
+   */
+  void mTmult (ScaLAPACKMatrix<NumberType> &C,
+               const ScaLAPACKMatrix<NumberType> &B,
+               const bool adding=false) const;
+
+  /**
+   * Matrix-matrix-multiplication using transpose of <tt>this</tt> and
+   * <tt>B</tt>.
+   *
+   * The optional parameter <tt>adding</tt> determines, whether the result is
+   * stored in <tt>C</tt> or added to <tt>C</tt>.
+   *
+   * if (adding) <i>C += A<sup>T</sup>*B<sup>T</sup></i>
+   *
+   * else <i>C = A<sup>T</sup>*B<sup>T</sup></i>
+   *
+   * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+   * <tt>C</tt> already has the right size.
+   *
+   * Following alignment conditions have to be fulfilled: MB_A = NB_B, NB_A = MB_C and MB_B = NB_C
+   */
+  void TmTmult (ScaLAPACKMatrix<NumberType> &C,
+                const ScaLAPACKMatrix<NumberType> &B,
+                const bool adding=false) const;
+
+  /**
    * Stores the distributed matrix in @p filename using HDF5.
    * In case that deal.II was built without HDF5
    * a call to this function will cause an exception to be thrown.
index a711a99f9e99a2ea88665113e7577ac48b4b6d6d..bd02083a49f26bbbee43d4f2346953c9cfcf2407 100644 (file)
@@ -365,11 +365,11 @@ extern "C"
                const int *n,
                const int *k,
                const double *alpha,
-               double *A,
+               const double *A,
                const int *IA,
                const int *JA,
                const int *DESCA,
-               double *B,
+               const double *B,
                const int *IB,
                const int *JB,
                const int *DESCB,
@@ -384,11 +384,11 @@ extern "C"
                const int *n,
                const int *k,
                const float *alpha,
-               float *A,
+               const float *A,
                const int *IA,
                const int *JA,
                const int *DESCA,
-               float *B,
+               const float *B,
                const int *IB,
                const int *JB,
                const int *DESCB,
@@ -1062,11 +1062,11 @@ inline void pgemm(const char *transa,
                   const int *n,
                   const int *k,
                   const double *alpha,
-                  double *A,
+                  const double *A,
                   const int *IA,
                   const int *JA,
                   const int *DESCA,
-                  double *B,
+                  const double *B,
                   const int *IB,
                   const int *JB,
                   const int *DESCB,
@@ -1085,11 +1085,11 @@ inline void pgemm(const char *transa,
                   const int *n,
                   const int *k,
                   const float *alpha,
-                  float *A,
+                  const float *A,
                   const int *IA,
                   const int *JA,
                   const int *DESCA,
-                  float *B,
+                  const float *B,
                   const int *IB,
                   const int *JB,
                   const int *DESCB,
index be2cd3391ff16a7c303315d51b2cfa0cdd49c0ac..1898d38c92c79e69edb6c40812623599638d4508 100644 (file)
@@ -499,6 +499,135 @@ void ScaLAPACKMatrix<NumberType>::Tadd(const NumberType a,
 {
   add(B,1,a,true);
 }
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::mult(ScaLAPACKMatrix<NumberType> &C,
+                                       const ScaLAPACKMatrix<NumberType> &B,
+                                       const NumberType alpha,
+                                       const NumberType beta,
+                                       const bool transpose_A,
+                                       const bool transpose_B) const
+{
+  Assert(this->grid==B.grid,ExcMessage("The matrices A and B need to have the same process grid"));
+  Assert(C.grid==B.grid,ExcMessage("The matrices B and C need to have the same process grid"));
+
+  // see for further info:
+  // https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgemm.htm
+  if (!transpose_A && !transpose_B)
+    {
+      Assert(this->n_columns==B.n_rows,ExcDimensionMismatch(this->n_columns,B.n_rows));
+      Assert(this->n_rows==C.n_rows,ExcDimensionMismatch(this->n_rows,C.n_rows));
+      Assert(B.n_columns==C.n_columns,ExcDimensionMismatch(B.n_columns,C.n_columns));
+      Assert(this->row_block_size==C.row_block_size,ExcDimensionMismatch(this->row_block_size,C.row_block_size));
+      Assert(this->column_block_size==B.row_block_size,ExcDimensionMismatch(this->column_block_size,B.row_block_size));
+      Assert(B.column_block_size==C.column_block_size,ExcDimensionMismatch(B.column_block_size,C.column_block_size));
+    }
+  else if (transpose_A && !transpose_B)
+    {
+      Assert(this->n_rows==B.n_rows,ExcDimensionMismatch(this->n_rows,B.n_rows));
+      Assert(this->n_columns==C.n_rows,ExcDimensionMismatch(this->n_columns,C.n_rows));
+      Assert(B.n_columns==C.n_columns,ExcDimensionMismatch(B.n_columns,C.n_columns));
+      Assert(this->column_block_size==C.row_block_size,ExcDimensionMismatch(this->column_block_size,C.row_block_size));
+      Assert(this->row_block_size==B.row_block_size,ExcDimensionMismatch(this->row_block_size,B.row_block_size));
+      Assert(B.column_block_size==C.column_block_size,ExcDimensionMismatch(B.column_block_size,C.column_block_size));
+    }
+  else if (!transpose_A && transpose_B)
+    {
+      Assert(this->n_columns==B.n_columns,ExcDimensionMismatch(this->n_columns,B.n_columns));
+      Assert(this->n_rows==C.n_rows,ExcDimensionMismatch(this->n_rows,C.n_rows));
+      Assert(B.n_rows==C.n_columns,ExcDimensionMismatch(B.n_rows,C.n_columns));
+      Assert(this->row_block_size==C.row_block_size,ExcDimensionMismatch(this->row_block_size,C.row_block_size));
+      Assert(this->column_block_size==B.column_block_size,ExcDimensionMismatch(this->column_block_size,B.column_block_size));
+      Assert(B.row_block_size==C.column_block_size,ExcDimensionMismatch(B.row_block_size,C.column_block_size));
+    }
+  else // if (transpose_A && transpose_B)
+    {
+      Assert(this->n_rows==B.n_columns,ExcDimensionMismatch(this->n_rows,B.n_columns));
+      Assert(this->n_columns==C.n_rows,ExcDimensionMismatch(this->n_columns,C.n_rows));
+      Assert(B.n_rows==C.n_columns,ExcDimensionMismatch(B.n_rows,C.n_columns));
+
+      Assert(this->column_block_size==C.row_block_size,ExcDimensionMismatch(this->row_block_size,C.row_block_size));
+      Assert(this->row_block_size==B.column_block_size,ExcDimensionMismatch(this->column_block_size,B.row_block_size));
+      Assert(B.row_block_size==C.column_block_size,ExcDimensionMismatch(B.column_block_size,C.column_block_size));
+    }
+  Threads::Mutex::ScopedLock lock (mutex);
+
+  if (this->grid->mpi_process_is_active)
+    {
+      char trans_a = transpose_A ? 'T' : 'N';
+      char trans_b = transpose_B ? 'T' : 'N';
+
+      const NumberType *A_loc = (this->values.size()>0) ? (&(this->values[0])) : nullptr;
+      const NumberType *B_loc = (B.values.size()>0) ? (&(B.values[0])) : nullptr;
+      NumberType *C_loc = (C.values.size()>0) ? (&(C.values[0])) : nullptr;
+      int m = C.n_rows;
+      int n = C.n_columns;
+      int k = transpose_A ? this->n_rows : this->n_columns;
+
+      pgemm(&trans_a,&trans_b,&m,&n,&k,
+            &alpha,A_loc,&(this->submatrix_row),&(this->submatrix_column),this->descriptor,
+            B_loc,&B.submatrix_row,&B.submatrix_column,B.descriptor,
+            &beta,C_loc,&C.submatrix_row,&C.submatrix_column,C.descriptor);
+    }
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::mmult(ScaLAPACKMatrix<NumberType> &C,
+                                        const ScaLAPACKMatrix<NumberType> &B,
+                                        const bool adding) const
+{
+  if (adding)
+    mult(C,B,1.,1.,false,false);
+  else
+    mult(C,B,1.,0.,false,false);
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::Tmmult(ScaLAPACKMatrix<NumberType> &C,
+                                         const ScaLAPACKMatrix<NumberType> &B,
+                                         const bool adding) const
+{
+  if (adding)
+    mult(C,B,1.,1.,true,false);
+  else
+    mult(C,B,1.,0.,true,false);
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::mTmult(ScaLAPACKMatrix<NumberType> &C,
+                                         const ScaLAPACKMatrix<NumberType> &B,
+                                         const bool adding) const
+{
+  if (adding)
+    mult(C,B,1.,1.,false,true);
+  else
+    mult(C,B,1.,0.,false,true);
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::TmTmult(ScaLAPACKMatrix<NumberType> &C,
+                                          const ScaLAPACKMatrix<NumberType> &B,
+                                          const bool adding) const
+{
+  if (adding)
+    mult(C,B,1.,1.,true,true);
+  else
+    mult(C,B,1.,0.,true,true);
+}
+
+
+
+template <typename NumberType>
 void ScaLAPACKMatrix<NumberType>::compute_cholesky_factorization()
 {
   Assert (n_columns == n_rows,
diff --git a/tests/scalapack/scalapack_12.cc b/tests/scalapack/scalapack_12.cc
new file mode 100644 (file)
index 0000000..a832eea
--- /dev/null
@@ -0,0 +1,218 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include "../tests.h"
+#include "../lapack/create_matrix.h"
+
+// test multiplication of distributed ScaLAPACKMatrices
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/multithread_info.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <typeinfo>
+
+
+template <typename NumberType>
+void test()
+{
+  MPI_Comm mpi_communicator(MPI_COMM_WORLD);
+  const unsigned int n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator));
+  const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));
+
+  std::cout << std::setprecision(10);
+  ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
+
+  const unsigned int proc_rows = std::floor(std::sqrt(n_mpi_processes));
+  const unsigned int proc_columns = std::floor(n_mpi_processes/proc_rows);
+  //create 2d process grid
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,proc_rows,proc_columns);
+  pcout << "2D process grid: " << grid->get_process_grid_rows() << "x" << grid->get_process_grid_columns() << std::endl << std::endl;
+
+  const std::vector<unsigned int> sizes = {{300,400,500}};
+
+  // test C = alpha A*B + beta C
+  {
+    FullMatrix<NumberType> full_A(sizes[0],sizes[2]);
+    FullMatrix<NumberType> full_B(sizes[2],sizes[1]);
+    FullMatrix<NumberType> full_C(sizes[0],sizes[1]);
+    create_random(full_A);
+    create_random(full_B);
+    create_random(full_C);
+
+    // conditions for block sizes: mb_A=mb_C, nb_B=nb_C, nb_A=mb_B
+    const unsigned int mb_A=32, nb_A=64, nb_B=16;
+    const unsigned int mb_C=mb_A, mb_B=nb_A;
+    const unsigned int nb_C=nb_B;
+
+    ScaLAPACKMatrix<NumberType> scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A);
+    ScaLAPACKMatrix<NumberType> scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B);
+    ScaLAPACKMatrix<NumberType> scalapack_C (full_C.m(),full_C.n(),grid,mb_C,nb_C);
+    scalapack_A = full_A;
+    scalapack_B = full_B;
+    scalapack_C = full_C;
+
+    const NumberType alpha = 1.4, beta = 0.1;
+
+    full_A *= alpha;
+    full_C *= beta;
+    full_A.mmult(full_C,full_B,true);
+
+    scalapack_A.mult(scalapack_C,scalapack_B,alpha,beta,false,false);
+    FullMatrix<NumberType> tmp_full_C(full_C.m(),full_C.n());
+    scalapack_C.copy_to(tmp_full_C);
+
+    pcout << "   computing C = alpha A * B + beta C with"
+          << " A in R^(" << scalapack_A.m() << "x" << scalapack_A.n() << "),"
+          << " B in R^(" << scalapack_B.m() << "x" << scalapack_B.n() << ") and"
+          << " C in R^(" << scalapack_C.m() << "x" << scalapack_C.n() << ")" << std::endl;
+    pcout << "   norms: " << tmp_full_C.frobenius_norm()<< " & "
+          << full_C.frobenius_norm() << "  for "
+          << typeid(NumberType).name() << std::endl << std::endl;
+  }
+  // test C = alpha A^T*B + beta C
+  {
+    FullMatrix<NumberType> full_A(sizes[2],sizes[0]);
+    FullMatrix<NumberType> full_B(sizes[2],sizes[1]);
+    FullMatrix<NumberType> full_C(sizes[0],sizes[1]);
+    create_random(full_A);
+    create_random(full_B);
+    create_random(full_C);
+
+    // conditions for block sizes: nb_A=mb_C, nb_B=nb_C, mb_A=mb_B
+    const unsigned int mb_A=32, nb_A=64, nb_B=16;
+    const unsigned int mb_B=mb_A, mb_C=nb_A;
+    const unsigned int nb_C=nb_B;
+
+    ScaLAPACKMatrix<NumberType> scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A);
+    ScaLAPACKMatrix<NumberType> scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B);
+    ScaLAPACKMatrix<NumberType> scalapack_C (full_C.m(),full_C.n(),grid,mb_C,nb_C);
+    scalapack_A = full_A;
+    scalapack_B = full_B;
+    scalapack_C = full_C;
+
+    const NumberType alpha = 1.4, beta = 0.1;
+
+    full_A *= alpha;
+    full_C *= beta;
+    full_A.Tmmult(full_C,full_B,true);
+
+    scalapack_A.mult(scalapack_C,scalapack_B,alpha,beta,true,false);
+    FullMatrix<NumberType> tmp_full_C(full_C.m(),full_C.n());
+    scalapack_C.copy_to(tmp_full_C);
+
+    pcout << "   computing C = alpha A^T * B + beta C with"
+          << " A in R^(" << scalapack_A.m() << "x" << scalapack_A.n() << "),"
+          << " B in R^(" << scalapack_B.m() << "x" << scalapack_B.n() << ") and"
+          << " C in R^(" << scalapack_C.m() << "x" << scalapack_C.n() << ")" << std::endl;
+    pcout << "   norms: " << tmp_full_C.frobenius_norm()<< " & "
+          << full_C.frobenius_norm() << "  for "
+          << typeid(NumberType).name() << std::endl << std::endl;
+  }
+  // test C = alpha A * B^T + beta C
+  {
+    FullMatrix<NumberType> full_A(sizes[0],sizes[2]);
+    FullMatrix<NumberType> full_B(sizes[1],sizes[2]);
+    FullMatrix<NumberType> full_C(sizes[0],sizes[1]);
+    create_random(full_A);
+    create_random(full_B);
+    create_random(full_C);
+
+    // conditions for block sizes: mb_A=mb_C, mb_B=nb_C, nb_A=nb_B
+    const unsigned int mb_A=32, nb_A=64, mb_B=16;
+    const unsigned int nb_B=nb_A, mb_C=mb_A;
+    const unsigned int nb_C=mb_B;
+
+    ScaLAPACKMatrix<NumberType> scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A);
+    ScaLAPACKMatrix<NumberType> scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B);
+    ScaLAPACKMatrix<NumberType> scalapack_C (full_C.m(),full_C.n(),grid,mb_C,nb_C);
+    scalapack_A = full_A;
+    scalapack_B = full_B;
+    scalapack_C = full_C;
+
+    const NumberType alpha = 1.4, beta = 0.1;
+
+    full_A *= alpha;
+    full_C *= beta;
+    full_A.mTmult(full_C,full_B,true);
+
+    scalapack_A.mult(scalapack_C,scalapack_B,alpha,beta,false,true);
+    FullMatrix<NumberType> tmp_full_C(full_C.m(),full_C.n());
+    scalapack_C.copy_to(tmp_full_C);
+
+    pcout << "   computing C = alpha A * B^T + beta C with"
+          << " A in R^(" << scalapack_A.m() << "x" << scalapack_A.n() << "),"
+          << " B in R^(" << scalapack_B.m() << "x" << scalapack_B.n() << ") and"
+          << " C in R^(" << scalapack_C.m() << "x" << scalapack_C.n() << ")" << std::endl;
+    pcout << "   norms: " << tmp_full_C.frobenius_norm()<< " & "
+          << full_C.frobenius_norm() << "  for "
+          << typeid(NumberType).name() << std::endl << std::endl;
+  }
+  // test C = alpha A^T * B^T + beta C
+  {
+    FullMatrix<NumberType> full_A(sizes[2],sizes[0]);
+    FullMatrix<NumberType> full_B(sizes[1],sizes[2]);
+    FullMatrix<NumberType> full_C(sizes[0],sizes[1]);
+    create_random(full_A);
+    create_random(full_B);
+    create_random(full_C);
+
+    // conditions for block sizes: nb_A=mb_C, mb_B=nb_C, mb_A=nb_B
+    const unsigned int mb_A=32, nb_A=64, mb_B=16;
+    const unsigned int nb_B=mb_A, mb_C=nb_A;
+    const unsigned int nb_C=mb_B;
+
+    ScaLAPACKMatrix<NumberType> scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A);
+    ScaLAPACKMatrix<NumberType> scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B);
+    ScaLAPACKMatrix<NumberType> scalapack_C (full_C.m(),full_C.n(),grid,mb_C,nb_C);
+    scalapack_A = full_A;
+    scalapack_B = full_B;
+    scalapack_C = full_C;
+
+    const NumberType alpha = 1.4, beta = 0.1;
+
+    full_A *= alpha;
+    full_C *= beta;
+    full_A.TmTmult(full_C,full_B,true);
+
+    scalapack_A.mult(scalapack_C,scalapack_B,alpha,beta,true,true);
+    FullMatrix<NumberType> tmp_full_C(full_C.m(),full_C.n());
+    scalapack_C.copy_to(tmp_full_C);
+
+    pcout << "   computing C = alpha A^T * B^T + beta C with"
+          << " A in R^(" << scalapack_A.m() << "x" << scalapack_A.n() << "),"
+          << " B in R^(" << scalapack_B.m() << "x" << scalapack_B.n() << ") and"
+          << " C in R^(" << scalapack_C.m() << "x" << scalapack_C.n() << ")" << std::endl;
+    pcout << "   norms: " << tmp_full_C.frobenius_norm()<< " & "
+          << full_C.frobenius_norm() << "  for "
+          << typeid(NumberType).name() << std::endl << std::endl;
+  }
+  pcout << std::endl;
+}
+
+
+
+int main (int argc,char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+  test<double>();
+}
diff --git a/tests/scalapack/scalapack_12.mpirun=1.output b/tests/scalapack/scalapack_12.mpirun=1.output
new file mode 100644 (file)
index 0000000..a65b495
--- /dev/null
@@ -0,0 +1,15 @@
+2D process grid: 1x1
+
+   computing C = alpha A * B + beta C with A in R^(300x500), B in R^(500x400) and C in R^(300x400)
+   norms: 60787.8449 & 60787.8449  for d
+
+   computing C = alpha A^T * B + beta C with A in R^(500x300), B in R^(500x400) and C in R^(300x400)
+   norms: 60655.07764 & 60655.07764  for d
+
+   computing C = alpha A * B^T + beta C with A in R^(300x500), B in R^(400x500) and C in R^(300x400)
+   norms: 60707.53954 & 60707.53954  for d
+
+   computing C = alpha A^T * B^T + beta C with A in R^(500x300), B in R^(400x500) and C in R^(300x400)
+   norms: 60757.09659 & 60757.09659  for d
+
+
diff --git a/tests/scalapack/scalapack_12.mpirun=11.output b/tests/scalapack/scalapack_12.mpirun=11.output
new file mode 100644 (file)
index 0000000..9646ef5
--- /dev/null
@@ -0,0 +1,15 @@
+2D process grid: 3x3
+
+   computing C = alpha A * B + beta C with A in R^(300x500), B in R^(500x400) and C in R^(300x400)
+   norms: 60787.8449 & 60787.8449  for d
+
+   computing C = alpha A^T * B + beta C with A in R^(500x300), B in R^(500x400) and C in R^(300x400)
+   norms: 60655.07764 & 60655.07764  for d
+
+   computing C = alpha A * B^T + beta C with A in R^(300x500), B in R^(400x500) and C in R^(300x400)
+   norms: 60707.53954 & 60707.53954  for d
+
+   computing C = alpha A^T * B^T + beta C with A in R^(500x300), B in R^(400x500) and C in R^(300x400)
+   norms: 60757.09659 & 60757.09659  for d
+
+
diff --git a/tests/scalapack/scalapack_12.mpirun=9.output b/tests/scalapack/scalapack_12.mpirun=9.output
new file mode 100644 (file)
index 0000000..9646ef5
--- /dev/null
@@ -0,0 +1,15 @@
+2D process grid: 3x3
+
+   computing C = alpha A * B + beta C with A in R^(300x500), B in R^(500x400) and C in R^(300x400)
+   norms: 60787.8449 & 60787.8449  for d
+
+   computing C = alpha A^T * B + beta C with A in R^(500x300), B in R^(500x400) and C in R^(300x400)
+   norms: 60655.07764 & 60655.07764  for d
+
+   computing C = alpha A * B^T + beta C with A in R^(300x500), B in R^(400x500) and C in R^(300x400)
+   norms: 60707.53954 & 60707.53954  for d
+
+   computing C = alpha A^T * B^T + beta C with A in R^(500x300), B in R^(400x500) and C in R^(300x400)
+   norms: 60757.09659 & 60757.09659  for d
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.