From: Benjamin Brands Date: Tue, 6 Feb 2018 15:47:04 +0000 (+0100) Subject: add routines for multiplication of ScaLAPACKMatrices + tests X-Git-Tag: v9.0.0-rc1~416^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3572a1eeff9344cd99074058207c1736f159a8e4;p=dealii.git add routines for multiplication of ScaLAPACKMatrices + tests --- diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index 934b2fea7b..788355e4a1 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -246,6 +246,108 @@ public: void Tadd(const NumberType b, const ScaLAPACKMatrix &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
NB_A = MB_B
NB_B = NB_C | C = a A * B + b C | + * | false | true | MB_A = MB_C
NB_A = NB_B
MB_B = NB_C | C = a A * BT + b C | + * | true | false | MB_A = MB_B
NB_A = MB_C
NB_B = NB_C | C = a AT * B + b C | + * | true | true | MB_A = NB_B
NB_A = MB_C
MB_B = NB_C | C = a AT * BT + b C | + * + * It is assumed that A and B have compatible sizes and that + * C already has the right size. + * + * The matrices A, B and C must have the same process grid. + */ + void mult(ScaLAPACKMatrix &C, + const ScaLAPACKMatrix &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 adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) C += A*B + * + * else C = A*B + * + * Assumes that A and B have compatible sizes and that + * C 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 &C, + const ScaLAPACKMatrix &B, + const bool adding=false) const; + + /** + * Matrix-matrix-multiplication using transpose of this. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) C += AT*B + * + * else C = AT*B + * + * Assumes that A and B have compatible sizes and that + * C 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 &C, + const ScaLAPACKMatrix &B, + const bool adding=false) const; + + /** + * Matrix-matrix-multiplication using transpose of B. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) C += A*BT + * + * else C = A*BT + * + * Assumes that A and B have compatible sizes and that + * C 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 &C, + const ScaLAPACKMatrix &B, + const bool adding=false) const; + + /** + * Matrix-matrix-multiplication using transpose of this and + * B. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) C += AT*BT + * + * else C = AT*BT + * + * Assumes that A and B have compatible sizes and that + * C 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 &C, + const ScaLAPACKMatrix &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. diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index a711a99f9e..bd02083a49 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -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, diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index be2cd3391f..1898d38c92 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -499,6 +499,135 @@ void ScaLAPACKMatrix::Tadd(const NumberType a, { add(B,1,a,true); } + + + +template +void ScaLAPACKMatrix::mult(ScaLAPACKMatrix &C, + const ScaLAPACKMatrix &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 +void ScaLAPACKMatrix::mmult(ScaLAPACKMatrix &C, + const ScaLAPACKMatrix &B, + const bool adding) const +{ + if (adding) + mult(C,B,1.,1.,false,false); + else + mult(C,B,1.,0.,false,false); +} + + + +template +void ScaLAPACKMatrix::Tmmult(ScaLAPACKMatrix &C, + const ScaLAPACKMatrix &B, + const bool adding) const +{ + if (adding) + mult(C,B,1.,1.,true,false); + else + mult(C,B,1.,0.,true,false); +} + + + +template +void ScaLAPACKMatrix::mTmult(ScaLAPACKMatrix &C, + const ScaLAPACKMatrix &B, + const bool adding) const +{ + if (adding) + mult(C,B,1.,1.,false,true); + else + mult(C,B,1.,0.,false,true); +} + + + +template +void ScaLAPACKMatrix::TmTmult(ScaLAPACKMatrix &C, + const ScaLAPACKMatrix &B, + const bool adding) const +{ + if (adding) + mult(C,B,1.,1.,true,true); + else + mult(C,B,1.,0.,true,true); +} + + + +template void ScaLAPACKMatrix::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 index 0000000000..a832eea9a5 --- /dev/null +++ b/tests/scalapack/scalapack_12.cc @@ -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 +#include +#include +#include +#include + +#include + +#include +#include +#include + + +template +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 grid = std::make_shared(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 sizes = {{300,400,500}}; + + // test C = alpha A*B + beta C + { + FullMatrix full_A(sizes[0],sizes[2]); + FullMatrix full_B(sizes[2],sizes[1]); + FullMatrix 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 scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A); + ScaLAPACKMatrix scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B); + ScaLAPACKMatrix 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 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 full_A(sizes[2],sizes[0]); + FullMatrix full_B(sizes[2],sizes[1]); + FullMatrix 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 scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A); + ScaLAPACKMatrix scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B); + ScaLAPACKMatrix 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 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 full_A(sizes[0],sizes[2]); + FullMatrix full_B(sizes[1],sizes[2]); + FullMatrix 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 scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A); + ScaLAPACKMatrix scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B); + ScaLAPACKMatrix 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 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 full_A(sizes[2],sizes[0]); + FullMatrix full_B(sizes[1],sizes[2]); + FullMatrix 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 scalapack_A (full_A.m(),full_A.n(),grid,mb_A,nb_A); + ScaLAPACKMatrix scalapack_B (full_B.m(),full_B.n(),grid,mb_B,nb_B); + ScaLAPACKMatrix 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 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(); +} diff --git a/tests/scalapack/scalapack_12.mpirun=1.output b/tests/scalapack/scalapack_12.mpirun=1.output new file mode 100644 index 0000000000..a65b495305 --- /dev/null +++ b/tests/scalapack/scalapack_12.mpirun=1.output @@ -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 index 0000000000..9646ef5461 --- /dev/null +++ b/tests/scalapack/scalapack_12.mpirun=11.output @@ -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 index 0000000000..9646ef5461 --- /dev/null +++ b/tests/scalapack/scalapack_12.mpirun=9.output @@ -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 + +