From 953b803fac8d48166ab1201eb72cd112983a5ade Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Tue, 23 Jan 2018 10:18:49 +0100 Subject: [PATCH] add routine copy_to to copy content of a submatrix of A to a submatrix of B --- include/deal.II/lac/scalapack.h | 16 +++++ source/lac/scalapack.cc | 108 +++++++++++++++++++++++++++++++- tests/scalapack/scalapack_07.cc | 17 ++++- 3 files changed, 139 insertions(+), 2 deletions(-) diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index f1182355e0..584e714b1c 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -166,6 +166,22 @@ public: void copy_to (ScaLAPACKMatrix &dest) const; /** + * Copy a submatrix (subset) of the distributed matrix A to a submatrix of the distributed matrix @p B. + * + * - The global row and column index of the first element of the submatrix A is provided by @p offset_A + * with row index=offset_A.first and column index=offset_A.second + * + * - The global row and column index of the first element of the submatrix B is provided by @p offset_B. + * with row index=offset_B.first and column index=offset_B.second + * + * - The dimension of the submatrix to be copied is given by @p submatrix_size + * with number of rows=submatrix_size.first and number of columns=submatrix_size.second + * . + */ + void copy_to(ScaLAPACKMatrix &B, + const std::pair &offset_A, + const std::pair &offset_B, + const std::pair &submatrix_size) const; /** * Stores the distributed matrix in @p filename using HDF5 diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index d2bebe2d12..6879cec96c 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -214,6 +214,112 @@ ScaLAPACKMatrix::copy_to (FullMatrix &matrix) const +template +void +ScaLAPACKMatrix::copy_to(ScaLAPACKMatrix &B, + const std::pair &offset_A, + const std::pair &offset_B, + const std::pair &submatrix_size) const +{ + //submatrix is empty + if (submatrix_size.first == 0 || submatrix_size.second == 0) + return; + + //range checking for matrix A + Assert (offset_A.first<(unsigned int)(n_rows-submatrix_size.first+1), + ExcIndexRange(offset_A.first,0,n_rows-submatrix_size.first+1)); + Assert (offset_A.second<(unsigned int)(n_columns-submatrix_size.second+1), + ExcIndexRange(offset_A.second,0,n_columns-submatrix_size.second+1)); + + //range checking for matrix B + Assert (offset_B.first<(unsigned int)(B.n_rows-submatrix_size.first+1), + ExcIndexRange(offset_B.first,0,B.n_rows-submatrix_size.first+1)); + Assert (offset_B.second<(unsigned int)(B.n_columns-submatrix_size.second+1), + ExcIndexRange(offset_B.second,0,B.n_columns-submatrix_size.second+1)); + + //Currently, copying of matrices will only be supported if A and B share the same MPI communicator + int ierr, comparison; + ierr = MPI_Comm_compare(grid->mpi_communicator,B.grid->mpi_communicator,&comparison); + Assert (comparison == MPI_IDENT,ExcMessage("Matrix A and B must have a common MPI Communicator")); + + /* + * The routine pgemr2d requires a BLACS context resembling at least the union of process grids + * described by the BLACS contexts of matrix A and B + */ + int union_blacs_context = Csys2blacs_handle(this->grid->mpi_communicator); + const char *order = "Col"; + int union_n_process_rows = Utilities::MPI::n_mpi_processes(this->grid->mpi_communicator); + int union_n_process_columns = 1; + Cblacs_gridinit(&union_blacs_context, order, union_n_process_rows, union_n_process_columns); + + + int n_grid_rows_A,n_grid_columns_A,my_row_A,my_column_A; + Cblacs_gridinfo(this->grid->blacs_context,&n_grid_rows_A,&n_grid_columns_A,&my_row_A,&my_column_A); + + //check whether process is in the BLACS context of matrix A + const bool in_context_A = (my_row_A>=0 && my_row_A=0 && my_column_Ablacs_context,&n_grid_rows_B,&n_grid_columns_B,&my_row_B,&my_column_B); + + //check whether process is in the BLACS context of matrix B + const bool in_context_B = (my_row_B>=0 && my_row_B=0 && my_column_B desc_A, desc_B; + + const NumberType *loc_vals_A = nullptr; + NumberType *loc_vals_B = nullptr; + + // Note: the function pgemr2d has to be called for all processes in the union BLACS context + // If the calling process is not part of the BLACS context of A, desc_A[1] has to be -1 + // and all other parameters do not have to be set + // If the calling process is not part of the BLACS context of B, desc_B[1] has to be -1 + // and all other parameters do not have to be set + if (in_context_A) + { + if (this->values.size() != 0) + loc_vals_A = & this->values[0]; + + for (unsigned int i=0; idescriptor[i]; + } + else + desc_A[1] =-1; + + if (in_context_B) + { + if (B.values.size() != 0) + loc_vals_B = & B.values[0]; + + for (unsigned int i=0; i void ScaLAPACKMatrix::copy_to (ScaLAPACKMatrix &dest) const @@ -252,7 +358,7 @@ ScaLAPACKMatrix::copy_to (ScaLAPACKMatrix &dest) const /* * The routine pgemr2d requires a BLACS context resembling at least the union of process grids - * described by the BLACS contexts of the source and destination matrix + * described by the BLACS contexts of matrix A and B */ int union_blacs_context = Csys2blacs_handle(mpi_communicator_union); const char *order = "Col"; diff --git a/tests/scalapack/scalapack_07.cc b/tests/scalapack/scalapack_07.cc index 2b2e54355a..19f508613b 100644 --- a/tests/scalapack/scalapack_07.cc +++ b/tests/scalapack/scalapack_07.cc @@ -16,7 +16,8 @@ #include "../tests.h" #include "../lapack/create_matrix.h" -// test copying of distributed ScaLAPACKMatrices using ScaLAPACK routine p_gemr2d +// test copying of distributed ScaLAPACKMatrices using ScaLAPACK routine p_gemr2d +// test copying submatrices of distributed ScaLAPACKMatrices using ScaLAPACK routine p_gemr2d #include #include @@ -89,6 +90,20 @@ void test(const unsigned int block_size_i, const unsigned int block_size_j) AssertThrow(test_2d.frobenius_norm() < 1e-12,ExcInternalError()); AssertThrow(test_1d.frobenius_norm() < 1e-12,ExcInternalError()); AssertThrow(test_one.frobenius_norm() < 1e-12,ExcInternalError()); + + //copying submatrices + unsigned int sub_size=100; + std::pair offset_A = std::make_pair(49,99); + std::pair offset_B = std::make_pair(0,0); + std::pair submatrix_size = std::make_pair(sub_size,sub_size); + ScaLAPACKMatrix scalapack_matrix_dest(sub_size,sub_size,grid_2d,block_size_j,block_size_i); + scalapack_matrix_2d.copy_to(scalapack_matrix_dest,offset_A,offset_B,submatrix_size); + FullMatrix dest (sub_size,sub_size); + scalapack_matrix_dest.copy_to(dest); + for (unsigned int i=0; i