From 12426c29d9bf12726efc0d3ae5123643b2459dd9 Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Mon, 20 Nov 2017 08:39:51 +0100 Subject: [PATCH] Added routine and tests to compute all eigenpairs of s.p.d matrices --- include/deal.II/lac/scalapack.h | 15 +- source/lac/scalapack.cc | 54 +++-- tests/scalapack/scalapack_06.cc | 3 - tests/scalapack/scalapack_06.mpirun=1.output | 36 ++-- tests/scalapack/scalapack_06.mpirun=4.output | 36 ++-- tests/scalapack/scalapack_07.cc | 199 +++++++++++++++++++ tests/scalapack/scalapack_07.mpirun=1.output | 54 +++++ tests/scalapack/scalapack_07.mpirun=4.output | 54 +++++ 8 files changed, 376 insertions(+), 75 deletions(-) create mode 100644 tests/scalapack/scalapack_07.cc create mode 100644 tests/scalapack/scalapack_07.mpirun=1.output create mode 100644 tests/scalapack/scalapack_07.mpirun=4.output diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index d9bce6af67..48fb3d6c22 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -40,15 +40,17 @@ template class ScaLAPACKMatrix; /** * A class taking care of setting up the process grid for BLACS * ScaLapack matrices will have shared pointer of this class to perform block-cyclic distribution + * + * @author Benjamin Brands, 2017 */ class ProcessGrid { public: - /** - * Declare class ScaLAPACK as friend to provide access to private members, e.g. the MPI Communicator - */ - template friend class ScaLAPACKMatrix; + /** + * Declare class ScaLAPACK as friend to provide access to private members, e.g. the MPI Communicator + */ + template friend class ScaLAPACKMatrix; /* * Constructor for a process grid for a given @p mpi_communicator @@ -193,7 +195,7 @@ protected: * @image html scalapack_invert.png * * @ingroup Matrix1 - * @author Denis Davydov, 2017 + * @author Denis Davydov, Benjamin Brands, 2017 */ template class ScaLAPACKMatrix : protected TransposeTable @@ -267,11 +269,14 @@ public: /** * Compute all eigenvalues of real symmetric matrix using pdsyev + * After successful computation the eigenvalues are stored in @p eigenvalues in ascending order */ void eigenvalues_symmetric (std::vector &eigenvalues); /** * Compute all eigenpairs of real symmetric matrix using pdsyev + * After successful computation the eigenvalues are stored in @p eigenvalues in ascending order + * The eigenvectors are stored in the columns of the matrix, therefore overwriting the original content of the matrix */ void eigenpairs_symmetric (std::vector &eigenvalues); diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index f9ccbc0526..a2283156de 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -235,6 +235,16 @@ extern "C" */ void pdsyev_(const char *jobz, const char *uplo, const int *m, double *A, const int *ia, const int *ja, int *desca, double *w, double *z, const int *iz, const int *jz, int *descz, double *work, const int *lwork, int *info); + + /* + * pdlacpy copies all or a part of a distributed matrix A to another + * distributed matrix B. No communication is performed, pdlacpy + * performs a local copy sub(A) := sub(B), where sub(A) denotes + * A(ia:ia+m-1,ja:ja+n-1) and sub(B) denotes B(ib:ib+m-1,jb:jb+n-1) + * + */ + void pdlacpy_(const char *uplo, const int *m, const int *n, double *A, const int *ia, const int *ja, int *desca, + double *B, const int *ib, const int *jb, int *descb); } @@ -489,11 +499,11 @@ ScaLAPACKMatrix::ScaLAPACKMatrix(const size_type rows, const size_ty template ScaLAPACKMatrix::ScaLAPACKMatrix(const std::pair &sizes, - std::shared_ptr process_grid, - const std::pair &block_sizes, - const LAPACKSupport::Property property) -: -ScaLAPACKMatrix(sizes.first,sizes.second,process_grid,block_sizes.first,block_sizes.first,property) + std::shared_ptr process_grid, + const std::pair &block_sizes, + const LAPACKSupport::Property property) + : + ScaLAPACKMatrix(sizes.first,sizes.second,process_grid,block_sizes.first,block_sizes.first,property) {} @@ -726,9 +736,6 @@ void ScaLAPACKMatrix::eigenvalues_symmetric(std::vector Z_loc, &Z.submatrix_row, &Z.submatrix_column, Z.descriptor, &work[0], &lwork, &info); AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pdsyev", info)); - - //Scalapack puts eigenvalues in ascending order --> reversing to obtain descending order - std::reverse (ev.begin(),ev.end()); } /* * send the eigenvalues to processors not being part of the process grid @@ -751,8 +758,8 @@ void ScaLAPACKMatrix::eigenpairs_symmetric(std::vector & Assert (properties == LAPACKSupport::symmetric, ExcMessage("Matrix has to be symmetric for this operation.")); - ScaLAPACKMatrix eigenvectors (n_rows, n_columns, grid, row_block_size, column_block_size, LAPACKSupport::Property::general); - eigenvectors.descriptor[1]=0; + ScaLAPACKMatrix eigenvectors (n_rows, grid, row_block_size); + eigenvectors.properties = properties; ev.resize (n_rows); const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(grid->mpi_communicator)); @@ -776,45 +783,30 @@ void ScaLAPACKMatrix::eigenpairs_symmetric(std::vector & NumberType *eigenvectors_loc = &eigenvectors.values[0]; work.resize(1); - pcout << "Starting workspace query" << std::endl; - - pcout << "Descriptor A: "; - for (unsigned int i=0; i<9; ++i) - pcout << " " << descriptor[i]; - pcout << std::endl; - pcout << "Descriptor Z: "; - for (unsigned int i=0; i<9; ++i) - pcout << " " << eigenvectors.descriptor[i]; - pcout << std::endl; - pdsyev_(&jobz, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, &ev[0], eigenvectors_loc, &eigenvectors.submatrix_row, &eigenvectors.submatrix_column, eigenvectors.descriptor, &work[0], &lwork, &info); - pcout << "info = " << info << std::endl << std::endl; - lwork=work[0]; work.resize (lwork); - pcout << "Starting computation" << std::endl; - pdsyev_(&jobz, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, &ev[0], eigenvectors_loc, &eigenvectors.submatrix_row, &eigenvectors.submatrix_column, eigenvectors.descriptor, &work[0], &lwork, &info); - pcout << "info = " << info << std::endl << std::endl; - AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pdsyev", info)); - //Scalapack puts eigenvalues in ascending order --> reversing to obtain descending order - std::reverse (ev.begin(),ev.end()); + //copy eigenvectors to original matrix + //as the temporary matrix eigenvectors has identical dimensions and block-cyclic distribution we simply swap the local array + this->values.swap(eigenvectors.values); } /* * send the eigenvalues to processors not being part of the process grid */ MPI_Bcast(&ev.front(),ev.size(),MPI_DOUBLE, 0/*from root*/, grid->mpi_communicator_inactive_with_root); - /* On exit, the lower triangle (if uplo='L') or the upper triangle (if uplo='U') of A, - * including the diagonal, is destroyed. Therefore, the matrix is unusable + /* + * On exit matrix A stores the eigenvectors in the columns */ + properties = LAPACKSupport::Property::general; state = LAPACKSupport::eigenvalues; } diff --git a/tests/scalapack/scalapack_06.cc b/tests/scalapack/scalapack_06.cc index cb9fca3838..352cfa80d4 100644 --- a/tests/scalapack/scalapack_06.cc +++ b/tests/scalapack/scalapack_06.cc @@ -110,9 +110,6 @@ void test(const unsigned int size, const unsigned int block_size) dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info); AssertThrow (info==0, LAPACKSupport::ExcErrorCode("syev", info)); - - //save eigenvalues_Lapack in descending order instead of ascending order - std::reverse (eigenvalues_Lapack.begin(),eigenvalues_Lapack.end()); } // Scalapack: scalapack_A = full_A; diff --git a/tests/scalapack/scalapack_06.mpirun=1.output b/tests/scalapack/scalapack_06.mpirun=1.output index 34374f1a96..3269f0428f 100644 --- a/tests/scalapack/scalapack_06.mpirun=1.output +++ b/tests/scalapack/scalapack_06.mpirun=1.output @@ -1,54 +1,54 @@ 32 32 1 1 First 5 ScaLapack eigenvalues -48.3154 34.7121 34.4886 34.4474 34.1038 +29.2445 29.432 29.5416 29.7204 29.8717 First 5 Lapack eigenvalues -48.3154 34.7121 34.4886 34.4474 34.1038 +29.2445 29.432 29.5416 29.7204 29.8717 64 32 1 1 First 5 ScaLapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 First 5 Lapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 64 64 1 1 First 5 ScaLapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 First 5 Lapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 120 32 1 1 First 5 ScaLapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 First 5 Lapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 120 64 1 1 First 5 ScaLapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 First 5 Lapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 320 32 1 1 First 5 ScaLapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 First 5 Lapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 320 64 1 1 First 5 ScaLapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 First 5 Lapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 640 32 1 1 First 5 ScaLapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 First 5 Lapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 640 64 1 1 First 5 ScaLapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 First 5 Lapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 diff --git a/tests/scalapack/scalapack_06.mpirun=4.output b/tests/scalapack/scalapack_06.mpirun=4.output index 034152e1f2..ff89f7ffdf 100644 --- a/tests/scalapack/scalapack_06.mpirun=4.output +++ b/tests/scalapack/scalapack_06.mpirun=4.output @@ -1,54 +1,54 @@ 32 32 1 1 First 5 ScaLapack eigenvalues -48.3154 34.7121 34.4886 34.4474 34.1038 +29.2445 29.432 29.5416 29.7204 29.8717 First 5 Lapack eigenvalues -48.3154 34.7121 34.4886 34.4474 34.1038 +29.2445 29.432 29.5416 29.7204 29.8717 64 32 2 2 First 5 ScaLapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 First 5 Lapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 64 64 1 1 First 5 ScaLapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 First 5 Lapack eigenvalues -95.9029 68.4041 68.061 67.7781 67.7537 +59.9263 60.0671 60.1596 60.3945 60.5197 120 32 2 2 First 5 ScaLapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 First 5 Lapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 120 64 2 2 First 5 ScaLapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 First 5 Lapack eigenvalues -180.662 126.026 125.642 125.511 125.285 +113.975 114.221 114.571 114.636 114.811 320 32 2 2 First 5 ScaLapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 First 5 Lapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 320 64 2 2 First 5 ScaLapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 First 5 Lapack eigenvalues -479.554 330.118 329.968 329.781 329.646 +309.946 310.139 310.275 310.403 310.451 640 32 2 2 First 5 ScaLapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 First 5 Lapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 640 64 2 2 First 5 ScaLapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 First 5 Lapack eigenvalues -959.748 654.219 654.077 654.023 654.015 +625.685 625.75 625.877 626.063 626.217 diff --git a/tests/scalapack/scalapack_07.cc b/tests/scalapack/scalapack_07.cc new file mode 100644 index 0000000000..a40fd9d41e --- /dev/null +++ b/tests/scalapack/scalapack_07.cc @@ -0,0 +1,199 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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" + +// test eigenvalues() + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +extern "C" //Some Lapack routines +{ + void dsyev_(char *jobz, char *uplo, int *n, double *A, int *lda, double *w, double *work, int *lwork, int *info); + void ssyev_(char *jobz, char *uplo, int *n, float *A, int *lda, float *w, float *work, int *lwork, int *info); +} + +template +void test(const unsigned int size, const unsigned int block_size) +{ + 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)); + + ConditionalOStream pcout (std::cout, (this_mpi_process ==0)); + + // test multiplication with random vectors + boost::random::mt19937 gen; + boost::random::uniform_01<> dist; + + // Create SPD matrices of requested size: + FullMatrix full_A(size); + std::vector lapack_A(size*size); + + std::pair sizes = std::make_pair(size,size), block_sizes = std::make_pair(block_size,block_size); + std::shared_ptr grid = std::make_shared(mpi_communicator,sizes,block_sizes); + + ScaLAPACKMatrix scalapack_A (sizes.first, grid, block_sizes.first); + + pcout << size << " " << block_size << " " << scalapack_A.get_process_grid_rows() << " " << scalapack_A.get_process_grid_columns() << std::endl; + { + full_A = 0.; + boost::random::mt19937 gen; + boost::random::uniform_01<> dist; + + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = i; j < size; ++j) + { + const double val = dist(gen); + Assert (val >= 0. && val <= 1., + ExcInternalError()); + if (i==j) + { + // since A(i,j) < 1 and + // a symmetric diagonally dominant matrix is SPD + full_A(i,j) = val + size; + lapack_A[i*size+j] = val+size; + } + else + { + full_A(i,j) = val; + full_A(j,i) = val; + lapack_A[i*size+j] = val; + lapack_A[j*size+i] = val; + } + } + } + std::vector eigenvalues_ScaLapack, eigenvalues_Lapack(size); + //Lapack as reference + { + int info; //Variable containing information about the successfull exit of the lapack routine + char jobz = 'V'; //'V': all eigenpairs of A are computed + char uplo = 'U'; //storage format of the matrix A; not so important as matrix is symmetric + int LDA = size; //leading dimension of the matrix A + int lwork; //length of vector/array work + std::vector work (1); + + //by setting lwork to -1 a workspace query for work is done + //as matrix is symmetric: LDA == size of matrix + lwork = -1; + dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info); + lwork=work[0]; + work.resize (lwork); + dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info); + + AssertThrow (info==0, LAPACKSupport::ExcErrorCode("syev", info)); + } + // Scalapack: + scalapack_A = full_A; + scalapack_A.eigenpairs_symmetric(eigenvalues_ScaLapack); + FullMatrix p_eigenvectors (size,size); + scalapack_A.copy_to(p_eigenvectors); + unsigned int n_eigenvalues = eigenvalues_ScaLapack.size(), max_n_eigenvalues=5; + + pcout << "First " << max_n_eigenvalues << " ScaLapack eigenvalues" << std::endl; + + for (unsigned int i=0; i s_eigenvectors (size,size); + for (int i=0; i p_ev(size), s_ev(size); + for (unsigned int j=0; j sizes = {{32,64,120,320,640}}; + const std::vector blocks = {{32,64}}; + + for (const auto &s : sizes) + for (const auto &b : blocks) + if (b <= s) + test(s,b); + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/scalapack/scalapack_07.mpirun=1.output b/tests/scalapack/scalapack_07.mpirun=1.output new file mode 100644 index 0000000000..3269f0428f --- /dev/null +++ b/tests/scalapack/scalapack_07.mpirun=1.output @@ -0,0 +1,54 @@ +32 32 1 1 +First 5 ScaLapack eigenvalues +29.2445 29.432 29.5416 29.7204 29.8717 +First 5 Lapack eigenvalues +29.2445 29.432 29.5416 29.7204 29.8717 + +64 32 1 1 +First 5 ScaLapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 +First 5 Lapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 + +64 64 1 1 +First 5 ScaLapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 +First 5 Lapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 + +120 32 1 1 +First 5 ScaLapack eigenvalues +113.975 114.221 114.571 114.636 114.811 +First 5 Lapack eigenvalues +113.975 114.221 114.571 114.636 114.811 + +120 64 1 1 +First 5 ScaLapack eigenvalues +113.975 114.221 114.571 114.636 114.811 +First 5 Lapack eigenvalues +113.975 114.221 114.571 114.636 114.811 + +320 32 1 1 +First 5 ScaLapack eigenvalues +309.946 310.139 310.275 310.403 310.451 +First 5 Lapack eigenvalues +309.946 310.139 310.275 310.403 310.451 + +320 64 1 1 +First 5 ScaLapack eigenvalues +309.946 310.139 310.275 310.403 310.451 +First 5 Lapack eigenvalues +309.946 310.139 310.275 310.403 310.451 + +640 32 1 1 +First 5 ScaLapack eigenvalues +625.685 625.75 625.877 626.063 626.217 +First 5 Lapack eigenvalues +625.685 625.75 625.877 626.063 626.217 + +640 64 1 1 +First 5 ScaLapack eigenvalues +625.685 625.75 625.877 626.063 626.217 +First 5 Lapack eigenvalues +625.685 625.75 625.877 626.063 626.217 + diff --git a/tests/scalapack/scalapack_07.mpirun=4.output b/tests/scalapack/scalapack_07.mpirun=4.output new file mode 100644 index 0000000000..ff89f7ffdf --- /dev/null +++ b/tests/scalapack/scalapack_07.mpirun=4.output @@ -0,0 +1,54 @@ +32 32 1 1 +First 5 ScaLapack eigenvalues +29.2445 29.432 29.5416 29.7204 29.8717 +First 5 Lapack eigenvalues +29.2445 29.432 29.5416 29.7204 29.8717 + +64 32 2 2 +First 5 ScaLapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 +First 5 Lapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 + +64 64 1 1 +First 5 ScaLapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 +First 5 Lapack eigenvalues +59.9263 60.0671 60.1596 60.3945 60.5197 + +120 32 2 2 +First 5 ScaLapack eigenvalues +113.975 114.221 114.571 114.636 114.811 +First 5 Lapack eigenvalues +113.975 114.221 114.571 114.636 114.811 + +120 64 2 2 +First 5 ScaLapack eigenvalues +113.975 114.221 114.571 114.636 114.811 +First 5 Lapack eigenvalues +113.975 114.221 114.571 114.636 114.811 + +320 32 2 2 +First 5 ScaLapack eigenvalues +309.946 310.139 310.275 310.403 310.451 +First 5 Lapack eigenvalues +309.946 310.139 310.275 310.403 310.451 + +320 64 2 2 +First 5 ScaLapack eigenvalues +309.946 310.139 310.275 310.403 310.451 +First 5 Lapack eigenvalues +309.946 310.139 310.275 310.403 310.451 + +640 32 2 2 +First 5 ScaLapack eigenvalues +625.685 625.75 625.877 626.063 626.217 +First 5 Lapack eigenvalues +625.685 625.75 625.877 626.063 626.217 + +640 64 2 2 +First 5 ScaLapack eigenvalues +625.685 625.75 625.877 626.063 626.217 +First 5 Lapack eigenvalues +625.685 625.75 625.877 626.063 626.217 + -- 2.39.5