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.
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,
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,
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,
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,
{
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,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+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
+
+
--- /dev/null
+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
+
+
--- /dev/null
+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
+
+