/**
* 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 <typename NumberType> friend class ScaLAPACKMatrix;
+ /**
+ * Declare class ScaLAPACK as friend to provide access to private members, e.g. the MPI Communicator
+ */
+ template <typename NumberType> friend class ScaLAPACKMatrix;
/*
* Constructor for a process grid for a given @p mpi_communicator
* @image html scalapack_invert.png
*
* @ingroup Matrix1
- * @author Denis Davydov, 2017
+ * @author Denis Davydov, Benjamin Brands, 2017
*/
template <typename NumberType>
class ScaLAPACKMatrix : protected TransposeTable<NumberType>
/**
* 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<NumberType> &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<NumberType> &eigenvalues);
*/
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);
}
template <typename NumberType>
ScaLAPACKMatrix<NumberType>::ScaLAPACKMatrix(const std::pair<size_type,size_type> &sizes,
- std::shared_ptr<ProcessGrid> process_grid,
- const std::pair<size_type,size_type> &block_sizes,
- const LAPACKSupport::Property property)
-:
-ScaLAPACKMatrix<NumberType>(sizes.first,sizes.second,process_grid,block_sizes.first,block_sizes.first,property)
+ std::shared_ptr<ProcessGrid> process_grid,
+ const std::pair<size_type,size_type> &block_sizes,
+ const LAPACKSupport::Property property)
+ :
+ ScaLAPACKMatrix<NumberType>(sizes.first,sizes.second,process_grid,block_sizes.first,block_sizes.first,property)
{}
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
Assert (properties == LAPACKSupport::symmetric,
ExcMessage("Matrix has to be symmetric for this operation."));
- ScaLAPACKMatrix<NumberType> eigenvectors (n_rows, n_columns, grid, row_block_size, column_block_size, LAPACKSupport::Property::general);
- eigenvectors.descriptor[1]=0;
+ ScaLAPACKMatrix<NumberType> 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));
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;
}
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;
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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <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 <boost/random/mersenne_twister.hpp>
+#include <boost/random/uniform_01.hpp>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+
+#include <deal.II/lac/scalapack.h>
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <memory>
+
+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 <typename NumberType>
+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<NumberType> full_A(size);
+ std::vector<NumberType> lapack_A(size*size);
+
+ std::pair<int,int> sizes = std::make_pair(size,size), block_sizes = std::make_pair(block_size,block_size);
+ std::shared_ptr<ProcessGrid> grid = std::make_shared<ProcessGrid>(mpi_communicator,sizes,block_sizes);
+
+ ScaLAPACKMatrix<NumberType> 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<NumberType> 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<double> 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<NumberType> 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<max_n_eigenvalues; ++i)
+ pcout << eigenvalues_ScaLapack[i] << " ";
+
+ pcout << std::endl << "First " << max_n_eigenvalues << " Lapack eigenvalues" << std::endl;
+
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
+ pcout << eigenvalues_Lapack[i] << " ";
+ pcout << std::endl;
+ pcout << std::endl;
+
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
+ AssertThrow ( std::fabs(eigenvalues_ScaLapack[i]-eigenvalues_Lapack[i]) < std::fabs(eigenvalues_Lapack[i])*1e-10, dealii::ExcInternalError());
+
+
+ FullMatrix<NumberType> s_eigenvectors (size,size);
+ for (int i=0; i<size; ++i)
+ for (int j=0; j<size; ++j)
+ s_eigenvectors(i,j) = lapack_A[i*size+j];
+
+ //product of eigenvectors computed using Lapack and ScaLapack has to be either 1 or -1
+ for (unsigned int i=0; i<size; ++i)
+ {
+ Vector<double> p_ev(size), s_ev(size);
+ for (unsigned int j=0; j<size; ++j)
+ {
+ p_ev[j] = p_eigenvectors(j,i);
+ s_ev[j] = s_eigenvectors(i,j);
+ }
+ double product = p_ev * s_ev;
+ AssertThrow (std::fabs(std::fabs(product)-1) < 1e-6, dealii::ExcInternalError());
+ }
+}
+
+
+
+int main (int argc,char **argv)
+{
+
+ try
+ {
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+ const std::vector<unsigned int> sizes = {{32,64,120,320,640}};
+ const std::vector<unsigned int> blocks = {{32,64}};
+
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<double>(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;
+ };
+}
--- /dev/null
+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
+
--- /dev/null
+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
+