const std::pair<int,int> &index_limits = std::make_pair(-1,-1),
const std::pair<NumberType,NumberType> &value_limits = std::make_pair(-1,-1));
+
+
+ /**
+ * Funcion to compute the singular value decomposition (SVD) of an
+ * M-by-N matrix A, optionally computing the left and/or right
+ * singular vectors. The SVD is written as A = U * SIGMA * transpose(V)
+ * where SIGMA is an M-by-N diagonal matrix, @p U is an M-by-M orthogonal matrix,
+ * and @p V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
+ * are the singular values of A and the columns of U and V are the
+ * corresponding left and right singular vectors, respectively. The
+ * singular values are returned in decreasing order and only the first min(M,N)
+ * columns of U and rows of VT = transpose(V) are computed.
+ * Upon return the content of the matrix is unusable.
+ * The matrix A must have identical block cyclic distribution for the rows and column
+ * If left singular vectors are required matrices A and U
+ * have to be constructed with the same process grid and block cyclic distribution.
+ * If right singular vectors are required matrices A and VT
+ * have to be constructed with the same process grid and block cyclic distribution.
+ */
+ std::vector<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> &U,
+ ScaLAPACKMatrix<NumberType> &VT,
+ const bool left_singluar_vectors=false,
+ const bool right_singluar_vectors=false);
+
/**
* Estimate the the condition number of a SPD matrix in the $l_1$-norm.
* The matrix has to be in the Cholesky state (see compute_cholesky_factorization()).
int *iclustr,
float *gap,
int *info);
+
+ /*
+ * PDGESVD computes the singular value decomposition (SVD) of an
+ * M-by-N matrix A, optionally computing the left and/or right
+ * singular vectors
+ */
+ //FIXME
+ void pdgesvd_(const char *jobu,
+ const char *jobvt,
+ const int *m,
+ const int *n,
+ double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ double *S,
+ double *U,
+ const int *iu,
+ const int *ju,
+ const int *descu,
+ double *VT,
+ const int *ivt,
+ const int *jvt,
+ const int *descvt,
+ double *work,
+ int *lwork,
+ int *info);
+ void psgesvd_(const char *jobu,
+ const char *jobvt,
+ const int *m,
+ const int *n,
+ float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ float *S,
+ float *U,
+ const int *iu,
+ const int *ju,
+ const int *descu,
+ float *VT,
+ const int *ivt,
+ const int *jvt,
+ const int *descvt,
+ float *work,
+ int *lwork,
+ int *info);
}
+
+
/*
* In the following we have template wrappers for the ScaLAPACK routines
* wrappers for other numeric types can be added in the future
}
+template <typename number>
+inline void pgesvd(const char *jobu,
+ const char *jobvt,
+ const int *m,
+ const int *n,
+ number *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ number *S,
+ number *U,
+ const int *iu,
+ const int *ju,
+ const int *descu,
+ number *VT,
+ const int *ivt,
+ const int *jvt,
+ const int *descvt,
+ number *work,
+ int *lwork,
+ int *info)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgesvd(const char *jobu,
+ const char *jobvt,
+ const int *m,
+ const int *n,
+ double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ double *S,
+ double *U,
+ const int *iu,
+ const int *ju,
+ const int *descu,
+ double *VT,
+ const int *ivt,
+ const int *jvt,
+ const int *descvt,
+ double *work,
+ int *lwork,
+ int *info)
+{
+ pdgesvd_(jobu,jobvt,m,n,A,ia,ja,desca,S,U,iu,ju,descu,VT,ivt,jvt,descvt,work,lwork,info);
+}
+
+inline void pgesvd(const char *jobu,
+ const char *jobvt,
+ const int *m,
+ const int *n,
+ float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ float *S,
+ float *U,
+ const int *iu,
+ const int *ju,
+ const int *descu,
+ float *VT,
+ const int *ivt,
+ const int *jvt,
+ const int *descvt,
+ float *work,
+ int *lwork,
+ int *info)
+{
+ psgesvd_(jobu,jobvt,m,n,A,ia,ja,desca,S,U,iu,ju,descu,VT,ivt,jvt,descvt,work,lwork,info);
+}
+
#endif // DEAL_II_WITH_SCALAPACK
#endif // dealii_scalapack_templates_h
template void Utilities::MPI::ProcessGrid::send_to_inactive<double>(double *, const int) const;
template void Utilities::MPI::ProcessGrid::send_to_inactive<float>(float *, const int) const;
+template void Utilities::MPI::ProcessGrid::send_to_inactive<int>(int *, const int) const;
DEAL_II_NAMESPACE_CLOSE
const std::pair<int,int> &eigenvalue_idx,
const std::pair<NumberType,NumberType> &eigenvalue_limits)
{
- //Assert(Utilities::MPI::n_mpi_processes(gird->mpi_communicator_inactive_with_root)<=1,
- // ExcMessage("For the computation of eigenpairs do not use a number of MPI processes that do not fit in a 2D process grid"));
-
Assert (state == LAPACKSupport::matrix,
ExcMessage("Matrix has to be in Matrix state before calling this function."));
Assert (property == LAPACKSupport::symmetric,
//ScaLAPACKMatrix<NumberType> eigenvectors (n_rows, grid, row_block_size);
eigenvectors->property = property;
+ // number of eigenvalues to be returned; upon successful exit ev contains the m seclected eigenvalues in ascending order
+ int m = n_rows;
std::vector<NumberType> ev(n_rows);
if (grid->mpi_process_is_active)
bool all_eigenpairs=true;
NumberType vl,vu;
int il,iu;
- // number of eigenvalues to be returned; upon successful exit ev contains the m seclected eigenvalues in ascending order
- int m = n_rows;
// number of eigenvectors to be returned;
// upon successful exit the first m=nz columns contain the selected eigenvectors (only if jobz=='V')
int nz;
while ((int)ev.size() > m)
ev.pop_back();
}
+ /*
+ * send number of computed eigenvalues to inactive processes
+ */
+ grid->send_to_inactive(&m, 1);
+
+ /*
+ * inactive processes have to resize array of eigenvalues
+ */
+ if (! grid->mpi_process_is_active)
+ ev.resize (m);
/*
* send the eigenvalues to processors not being part of the process grid
*/
}
+template <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix<NumberType> &U,
+ ScaLAPACKMatrix<NumberType> &VT,
+ const bool left_singluar_vectors,
+ const bool right_singluar_vectors)
+{
+ Assert (state == LAPACKSupport::matrix,
+ ExcMessage("Matrix has to be in Matrix state before calling this function."));
+ Assert(row_block_size==column_block_size,ExcDimensionMismatch(n_rows,U.n_rows));
+
+ if (left_singluar_vectors)
+ {
+ Assert(n_rows==U.n_rows,ExcDimensionMismatch(n_rows,U.n_rows));
+ Assert(U.n_rows==U.n_columns,ExcDimensionMismatch(U.n_rows,U.n_columns));
+ Assert(row_block_size==U.row_block_size,ExcDimensionMismatch(row_block_size,U.row_block_size));
+ Assert(column_block_size==U.column_block_size,ExcDimensionMismatch(column_block_size,U.column_block_size));
+ Assert(grid->blacs_context==U.grid->blacs_context,ExcDimensionMismatch(grid->blacs_context,U.grid->blacs_context));
+ }
+ if (right_singluar_vectors)
+ {
+ Assert(n_columns==VT.n_rows,ExcDimensionMismatch(n_columns,VT.n_rows));
+ Assert(VT.n_rows==VT.n_columns,ExcDimensionMismatch(VT.n_rows,VT.n_columns));
+ Assert(row_block_size==VT.row_block_size,ExcDimensionMismatch(row_block_size,VT.row_block_size));
+ Assert(column_block_size==VT.column_block_size,ExcDimensionMismatch(column_block_size,VT.column_block_size));
+ Assert(grid->blacs_context==VT.grid->blacs_context,ExcDimensionMismatch(grid->blacs_context,VT.grid->blacs_context));
+ }
+ Threads::Mutex::ScopedLock lock (mutex);
+
+ std::vector<NumberType> sv(std::min(n_rows,n_columns));
+
+ if (grid->mpi_process_is_active)
+ {
+ char jobu = left_singluar_vectors ? 'V' : 'N';
+ char jobvt = right_singluar_vectors ? 'V' : 'N';
+ NumberType *A_loc = &this->values[0];
+ NumberType *U_loc = left_singluar_vectors ? &(U.values[0]) : NULL;
+ NumberType *VT_loc = right_singluar_vectors ? &(VT.values[0]) : NULL;
+ int info = 0;
+ /*
+ * by setting lwork to -1 a workspace query for optimal length of work is performed
+ */
+ int lwork=-1;
+ work.resize(1);
+
+ pgesvd(&jobu,&jobvt,&n_rows,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,
+ & *sv.begin(),U_loc,&U.submatrix_row,&U.submatrix_column,U.descriptor,
+ VT_loc,&VT.submatrix_row,&VT.submatrix_column,VT.descriptor,
+ &work[0],&lwork,&info);
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgesvd", info));
+
+ lwork=work[0];
+ work.resize(lwork);
+
+ pgesvd(&jobu,&jobvt,&n_rows,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,
+ & *sv.begin(),U_loc,&U.submatrix_row,&U.submatrix_column,U.descriptor,
+ VT_loc,&VT.submatrix_row,&VT.submatrix_column,VT.descriptor,
+ &work[0],&lwork,&info);
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgesvd", info));
+ }
+
+ /*
+ * send the singular values to processors not being part of the process grid
+ */
+ grid->send_to_inactive(sv.data(), sv.size());
+
+ property = LAPACKSupport::Property::general;
+ state = LAPACKSupport::State::unusable;
+
+ return sv;
+}
+
+
template <typename NumberType>
NumberType ScaLAPACKMatrix<NumberType>::reciprocal_condition_number(const NumberType a_norm) const
if (i==j)
// since A(i,j) < 1 and
// a symmetric diagonally dominant matrix is SPD
- A(i,j) = val + size + i*i/size;
+ A(i,j) = val + size;
else
{
A(i,j) = val;
#include "../tests.h"
#include "../lapack/create_matrix.h"
-// test eigenvalues_symmetric()
+// test eigenpairs_symmetric(const bool, const std::pair<int,int>&, const std::pair<NumberType,NumberType>&)
#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/base/process_grid.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <fstream>
#include <iostream>
#include <algorithm>
-
+#include <memory>
template <typename NumberType>
ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+
+ pcout << size << " " << block_size << std::endl;
+
// Create SPD matrices of requested size:
FullMatrix<NumberType> full_A(size);
std::vector<NumberType> lapack_A(size*size);
- std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
- ScaLAPACKMatrix<NumberType> scalapack_A (size, grid, block_size,
- LAPACKSupport::Property::symmetric);
-
- pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
create_spd (full_A);
-
for (unsigned int i = 0; i < size; ++i)
for (unsigned int j = 0; j < size; ++j)
lapack_A[i*size+j] = full_A(i,j);
- std::vector<NumberType> eigenvalues_ScaLapack, eigenvalues_Lapack(size);
+ std::vector<NumberType> eigenvalues_Lapack(size);
//Lapack as reference
{
int info; //Variable containing information about the successfull exit of the lapack routine
- char jobz = 'N'; //'N': the eigenvalues_Lapack of A are computed
+ 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
AssertThrow (info==0, LAPACKSupport::ExcErrorCode("syev", info));
}
- // Scalapack:
- scalapack_A = full_A;
- eigenvalues_ScaLapack = scalapack_A.eigenvalues_symmetric();
- unsigned int n_eigenvalues = eigenvalues_ScaLapack.size(), max_n_eigenvalues=5;
+ unsigned int n_eigenvalues = eigenvalues_Lapack.size(), max_n_eigenvalues=5;
+
+ std::vector<Vector<NumberType>> s_eigenvectors_ (max_n_eigenvalues,Vector<NumberType>(size));
+ for (int i=0; i<max_n_eigenvalues; ++i)
+ for (int j=0; j<size; ++j)
+ s_eigenvectors_[i][j] = lapack_A[(size-1-i)*size+j];
+
+ pcout << "comparing " << max_n_eigenvalues << " eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:" << std::endl;
+ std::vector<NumberType> eigenvalues_psyev;
+ ScaLAPACKMatrix<NumberType> scalapack_syev (size, grid, block_size);
+ scalapack_syev.set_property(LAPACKSupport::Property::symmetric);
+ scalapack_syev = full_A;
+ eigenvalues_psyev = scalapack_syev.eigenpairs_symmetric(true);
+ FullMatrix<NumberType> p_eigenvectors (size,size);
+ scalapack_syev.copy_to(p_eigenvectors);
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
+ {
+ AssertThrow ( std::abs(eigenvalues_psyev[n_eigenvalues-i-1]-eigenvalues_Lapack[n_eigenvalues-i-1]) / std::abs(eigenvalues_Lapack[n_eigenvalues-i-1]) < tol,
+ dealii::ExcInternalError());
+ }
+ pcout << " with respect to the given tolerance the eigenvalues coincide" << std::endl;
- pcout << "comparing " << max_n_eigenvalues << " eigenvalues computed using LAPACK and ScaLAPACK:" << std::endl;
+ std::vector<Vector<NumberType>> p_eigenvectors_ (max_n_eigenvalues,Vector<NumberType>(size));
for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- AssertThrow ( std::abs(eigenvalues_ScaLapack[n_eigenvalues-i-1]-eigenvalues_Lapack[n_eigenvalues-i-1]) / std::abs(eigenvalues_Lapack[n_eigenvalues-i-1]) < tol, dealii::ExcInternalError());
- pcout << " with respect to the given tolerance the eigenvalues coincide" << std::endl << std::endl;
+ for (unsigned int j=0; j<size; ++j)
+ p_eigenvectors_[i][j] = p_eigenvectors(j,size-1-i);
+
+ //product of eigenvectors computed using Lapack and ScaLapack has to be either 1 or -1
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
+ {
+ NumberType product = p_eigenvectors_[i] * s_eigenvectors_[i];
+
+ //the requirement for alignment of the eigenvectors has to be released (primarily for floats)
+ AssertThrow (std::abs(std::abs(product)-1) < tol*10, dealii::ExcInternalError());
+ }
+ pcout << " with respect to the given tolerance also the eigenvectors coincide" << std::endl << std::endl;
+
+ pcout << "comparing " << max_n_eigenvalues << " eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:" << std::endl;
+ std::vector<NumberType> eigenvalues_psyevx_partial;
+ ScaLAPACKMatrix<NumberType> scalapack_syevx_partial (size, grid, block_size);
+ scalapack_syevx_partial.set_property(LAPACKSupport::Property::symmetric);
+ scalapack_syevx_partial = full_A;
+ eigenvalues_psyevx_partial = scalapack_syevx_partial.eigenpairs_symmetric(true, std::make_pair(size-max_n_eigenvalues+1,size));
+ scalapack_syevx_partial.copy_to(p_eigenvectors);
+ for (unsigned int i=eigenvalues_psyevx_partial.size()-1; i>0; --i)
+ {
+ if ( !(std::abs(eigenvalues_psyevx_partial[i]-eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i]) /
+ std::abs(eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i]) < tol))
+ {
+ std::cout << "process #" << this_mpi_process << ": eigenvalues do not fit: "
+ << eigenvalues_psyevx_partial[i] << " <--> " << eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i] << std::endl;
+ }
+
+ AssertThrow ( std::abs(eigenvalues_psyevx_partial[i]-eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i]) /
+ std::abs(eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i]) < tol,
+ dealii::ExcInternalError());
+ }
+ pcout << " with respect to the given tolerance the eigenvalues coincide" << std::endl;
+
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
+ for (unsigned int j=0; j<size; ++j)
+ p_eigenvectors_[i][j] = p_eigenvectors(j,max_n_eigenvalues-1-i);
+
+ //product of eigenvectors computed using Lapack and ScaLapack has to be either 1 or -1
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
+ {
+ NumberType product = p_eigenvectors_[i] * s_eigenvectors_[i];
+
+ if (!(std::abs(std::abs(product)-1) < tol*10))
+ {
+ std::cout << "process #" << this_mpi_process << ": eigenvectors do not coincide: abs(" << product << ") != 1" << std::endl;
+ }
+
+ //the requirement for alignment of the eigenvectors has to be released (primarily for floats)
+ AssertThrow (std::abs(std::abs(product)-1) < tol*10, dealii::ExcInternalError());
+ }
+ pcout << " with respect to the given tolerance also the eigenvectors coincide" << std::endl << std::endl;
+
+ pcout << std::endl;
}
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> sizes = {{200,400,600}};
const std::vector<unsigned int> blocks = {{32,64}};
- const double tol_double = 1e-10;
- const float tol_float = 1e-5;
-
- for (const auto &s : sizes)
- for (const auto &b : blocks)
- if (b <= s)
- test<float>(s,b,tol_float);
+ const double tol = 1e-10;
for (const auto &s : sizes)
for (const auto &b : blocks)
if (b <= s)
- test<double>(s,b,tol_double);
+ test<double>(s,b,tol);
}
-32 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+200 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-120 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+200 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-120 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-320 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-320 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+400 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-32 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+400 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-120 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+600 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-120 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-320 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-320 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+600 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-32 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+200 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-120 32 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+200 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-120 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-320 32 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-320 64 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+400 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 32 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 64 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-32 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+400 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-120 32 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+600 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-120 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-320 32 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-320 64 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+600 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 32 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 64 3 3
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-32 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+200 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-120 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+200 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-120 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-320 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-320 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+400 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-32 32 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+400 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-120 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+600 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-120 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-320 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-320 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+600 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 32 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
+ with respect to the given tolerance also the eigenvectors coincide
-640 64 2 2
-comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
-32 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+200 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-64 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-120 32 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+200 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-120 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-320 32 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-320 64 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+400 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-640 32 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-640 64 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-32 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+400 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-64 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-64 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-120 32 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+600 32
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-120 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-320 32 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-320 64 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+600 64
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-640 32 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
+comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:
with respect to the given tolerance the eigenvalues coincide
with respect to the given tolerance also the eigenvectors coincide
-640 64 3 3
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
#include "../tests.h"
#include "../lapack/create_matrix.h"
-// test eigenpairs_symmetric() and calls set_property()
+// test copying of distributed ScaLAPACKMatrices using ScaLAPACK routine p_gemr2d
#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/base/process_grid.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/lapack_full_matrix.h>
-#include <deal.II/lac/lapack_templates.h>
+
#include <deal.II/lac/scalapack.h>
#include <fstream>
#include <iostream>
-#include <algorithm>
-#include <memory>
template <typename NumberType>
-void test(const unsigned int size, const unsigned int block_size, const NumberType tol)
+void test(const unsigned int block_size_i, const unsigned int block_size_j)
{
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
const unsigned int n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator));
ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
- // Create SPD matrices of requested size:
- FullMatrix<NumberType> full_A(size);
- std::vector<NumberType> lapack_A(size*size);
-
- std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
- ScaLAPACKMatrix<NumberType> scalapack_A (size, grid, block_size);
- scalapack_A.set_property(LAPACKSupport::Property::symmetric);
-
- pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
-
- create_spd (full_A);
- for (unsigned int i = 0; i < size; ++i)
- for (unsigned int j = 0; j < size; ++j)
- lapack_A[i*size+j] = full_A(i,j);
-
- 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<NumberType> work (1);
-
- //by setting lwork to -1 a workspace query for work is done
- //as matrix is symmetric: LDA == size of matrix
- lwork = -1;
- syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
- lwork=work[0];
- work.resize (lwork);
- syev(&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;
- eigenvalues_ScaLapack = scalapack_A.eigenpairs_symmetric();
- FullMatrix<NumberType> p_eigenvectors (size,size);
- scalapack_A.copy_to(p_eigenvectors);
- unsigned int n_eigenvalues = eigenvalues_ScaLapack.size(), max_n_eigenvalues=5;
-
- pcout << "comparing " << max_n_eigenvalues << " eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:" << std::endl;
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- AssertThrow ( std::abs(eigenvalues_ScaLapack[n_eigenvalues-i-1]-eigenvalues_Lapack[n_eigenvalues-i-1]) / std::abs(eigenvalues_Lapack[n_eigenvalues-i-1]) < tol, dealii::ExcInternalError());
- pcout << " with respect to the given tolerance the eigenvalues coincide" << std::endl;
-
- 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];
-
- std::vector<Vector<NumberType>> s_eigenvectors_ (max_n_eigenvalues,Vector<NumberType>(size));
- for (int i=0; i<max_n_eigenvalues; ++i)
- for (int j=0; j<size; ++j)
- s_eigenvectors_[i][j] = lapack_A[i*size+j];
-
- std::vector<Vector<NumberType>> p_eigenvectors_ (max_n_eigenvalues,Vector<NumberType>(size));
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- for (unsigned int j=0; j<size; ++j)
- p_eigenvectors_[i][j] = p_eigenvectors(j,i);
-
-
- //product of eigenvectors computed using Lapack and ScaLapack has to be either 1 or -1
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- {
- NumberType product = p_eigenvectors_[i] * s_eigenvectors_[i];
-
- //the requirement for alignment of the eigenvectors has to be released (primarily for floats)
- AssertThrow (std::abs(std::abs(product)-1) < tol*10, dealii::ExcInternalError());
- }
- pcout << " with respect to the given tolerance also the eigenvectors coincide" << std::endl << std::endl;
+ const unsigned int size = 500;
+ //create FullMatrix and fill it
+ FullMatrix<NumberType> full(size);
+ unsigned int count=0;
+ for (unsigned int i=0; i<size; ++i)
+ for (unsigned int j=0; j<size; ++j, ++count)
+ full(i,j) = count;
+
+ //create 2d process grid
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid_2d = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size_i,block_size_i);
+ //create 1d process grid
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid_1d = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,n_mpi_processes,1);
+ //create process grid containing one process
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid_single = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,1,1);
+
+ ScaLAPACKMatrix<NumberType> scalapack_matrix_2d (size, size, grid_2d, block_size_i, block_size_i);
+ ScaLAPACKMatrix<NumberType> scalapack_matrix_1d (size, size, grid_1d, block_size_j, block_size_j);
+ ScaLAPACKMatrix<NumberType> scalapack_matrix_single (size, size, grid_single, block_size_i, block_size_j);
+ ScaLAPACKMatrix<NumberType> scalapack_matrix_source (size, size, grid_single, block_size_j, block_size_i);
+
+ pcout << "2D grid matrix: dim=" << scalapack_matrix_2d.m() << "x" << scalapack_matrix_2d.n()
+ << "; blocks=" << block_size_i << "x" << block_size_i
+ << "; grid=" << grid_2d->get_process_grid_rows() << "x" << grid_2d->get_process_grid_columns() << std::endl;
+
+ pcout << "1D grid matrix: " << scalapack_matrix_1d.m() << "x" << scalapack_matrix_1d.n()
+ << "; blocks=" << block_size_j << "x" << block_size_j
+ << "; grid=" << grid_1d->get_process_grid_rows() << "x" << grid_1d->get_process_grid_columns() << std::endl;
+
+ pcout << "single process matrix: " << scalapack_matrix_single.m() << "x" << scalapack_matrix_single.n()
+ << "; blocks=" << block_size_i << "x" << block_size_j
+ << "; grid=" << grid_single->get_process_grid_rows() << "x" << grid_single->get_process_grid_columns() << std::endl << std::endl;
+
+ scalapack_matrix_source = full;
+
+ scalapack_matrix_source.copy_to(scalapack_matrix_2d);
+ scalapack_matrix_source.copy_to(scalapack_matrix_1d);
+ scalapack_matrix_source.copy_to(scalapack_matrix_single);
+
+ FullMatrix<NumberType> test_2d(size), test_1d(size), test_one(size);
+ scalapack_matrix_2d.copy_to(test_2d);
+ scalapack_matrix_1d.copy_to(test_1d);
+ scalapack_matrix_single.copy_to(test_one);
+ test_2d.add (-1,full);
+ test_1d.add (-1,full);
+ test_one.add (-1,full);
+
+ AssertThrow(test_2d.frobenius_norm() < 1e-12,ExcInternalError());
+ AssertThrow(test_1d.frobenius_norm() < 1e-12,ExcInternalError());
+ AssertThrow(test_one.frobenius_norm() < 1e-12,ExcInternalError());
}
int main (int argc,char **argv)
{
-
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}};
-
- const double tol_double = 1e-10;
- const float tol_float = 1e-5;
+ const std::vector<unsigned int> blocks_i = {{16,32,64}};
+ const std::vector<unsigned int> blocks_j = {{16,32,64}};
- for (const auto &s : sizes)
- for (const auto &b : blocks)
- if (b <= s)
- test<float>(s,b,tol_float);
+ for (const auto &s : blocks_i)
+ for (const auto &b : blocks_j)
+ test<float>(s,b);
- for (const auto &s : sizes)
- for (const auto &b : blocks)
- if (b <= s)
- test<double>(s,b,tol_double);
+ for (const auto &s : blocks_i)
+ for (const auto &b : blocks_j)
+ test<double>(s,b);
}
-32 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-32 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
+2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
+1D grid matrix: 500x500; blocks=16x16; grid=1x1
+single process matrix: 500x500; blocks=16x16; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
+1D grid matrix: 500x500; blocks=32x32; grid=1x1
+single process matrix: 500x500; blocks=16x32; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
+1D grid matrix: 500x500; blocks=64x64; grid=1x1
+single process matrix: 500x500; blocks=16x64; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
+1D grid matrix: 500x500; blocks=16x16; grid=1x1
+single process matrix: 500x500; blocks=32x16; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
+1D grid matrix: 500x500; blocks=32x32; grid=1x1
+single process matrix: 500x500; blocks=32x32; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
+1D grid matrix: 500x500; blocks=64x64; grid=1x1
+single process matrix: 500x500; blocks=32x64; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
+1D grid matrix: 500x500; blocks=16x16; grid=1x1
+single process matrix: 500x500; blocks=64x16; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
+1D grid matrix: 500x500; blocks=32x32; grid=1x1
+single process matrix: 500x500; blocks=64x32; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
+1D grid matrix: 500x500; blocks=64x64; grid=1x1
+single process matrix: 500x500; blocks=64x64; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
+1D grid matrix: 500x500; blocks=16x16; grid=1x1
+single process matrix: 500x500; blocks=16x16; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
+1D grid matrix: 500x500; blocks=32x32; grid=1x1
+single process matrix: 500x500; blocks=16x32; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
+1D grid matrix: 500x500; blocks=64x64; grid=1x1
+single process matrix: 500x500; blocks=16x64; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
+1D grid matrix: 500x500; blocks=16x16; grid=1x1
+single process matrix: 500x500; blocks=32x16; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
+1D grid matrix: 500x500; blocks=32x32; grid=1x1
+single process matrix: 500x500; blocks=32x32; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
+1D grid matrix: 500x500; blocks=64x64; grid=1x1
+single process matrix: 500x500; blocks=32x64; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
+1D grid matrix: 500x500; blocks=16x16; grid=1x1
+single process matrix: 500x500; blocks=64x16; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
+1D grid matrix: 500x500; blocks=32x32; grid=1x1
+single process matrix: 500x500; blocks=64x32; grid=1x1
+
+2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
+1D grid matrix: 500x500; blocks=64x64; grid=1x1
+single process matrix: 500x500; blocks=64x64; grid=1x1
+++ /dev/null
-32 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-32 32 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-64 64 1 1
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-120 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-320 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 32 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
-640 64 2 2
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
- with respect to the given tolerance the eigenvalues coincide
- with respect to the given tolerance also the eigenvectors coincide
-
#include "../tests.h"
#include "../lapack/create_matrix.h"
-// test copying of distributed ScaLAPACKMatrices using ScaLAPACK routine p_gemr2d
+// test compute_SVD(ScaLAPACKMatrix<NumberType>&,ScaLAPACKMatrix<NumberType>&,const bool,const bool)
#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/base/process_grid.h>
-
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/lapack_templates.h>
#include <deal.II/lac/scalapack.h>
#include <fstream>
#include <iostream>
+#include <algorithm>
+#include <memory>
template <typename NumberType>
-void test(const unsigned int block_size_i, const unsigned int block_size_j)
+void test(const unsigned int size, const unsigned int block_size, const NumberType tol)
{
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
const unsigned int n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator));
ConditionalOStream pcout (std::cout, (this_mpi_process ==0));
- const unsigned int size = 500;
- //create FullMatrix and fill it
- FullMatrix<NumberType> full(size);
- unsigned int count=0;
- for (unsigned int i=0; i<size; ++i)
- for (unsigned int j=0; j<size; ++j, ++count)
- full(i,j) = count;
-
- //create 2d process grid
- std::shared_ptr<Utilities::MPI::ProcessGrid> grid_2d = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size_i,block_size_i);
- //create 1d process grid
- std::shared_ptr<Utilities::MPI::ProcessGrid> grid_1d = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,n_mpi_processes,1);
- //create process grid containing one process
- std::shared_ptr<Utilities::MPI::ProcessGrid> grid_single = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,1,1);
-
- ScaLAPACKMatrix<NumberType> scalapack_matrix_2d (size, size, grid_2d, block_size_i, block_size_i);
- ScaLAPACKMatrix<NumberType> scalapack_matrix_1d (size, size, grid_1d, block_size_j, block_size_j);
- ScaLAPACKMatrix<NumberType> scalapack_matrix_single (size, size, grid_single, block_size_i, block_size_j);
- ScaLAPACKMatrix<NumberType> scalapack_matrix_source (size, size, grid_single, block_size_j, block_size_i);
-
- pcout << "2D grid matrix: dim=" << scalapack_matrix_2d.m() << "x" << scalapack_matrix_2d.n()
- << "; blocks=" << block_size_i << "x" << block_size_i
- << "; grid=" << grid_2d->get_process_grid_rows() << "x" << grid_2d->get_process_grid_columns() << std::endl;
-
- pcout << "1D grid matrix: " << scalapack_matrix_1d.m() << "x" << scalapack_matrix_1d.n()
- << "; blocks=" << block_size_j << "x" << block_size_j
- << "; grid=" << grid_1d->get_process_grid_rows() << "x" << grid_1d->get_process_grid_columns() << std::endl;
-
- pcout << "single process matrix: " << scalapack_matrix_single.m() << "x" << scalapack_matrix_single.n()
- << "; blocks=" << block_size_i << "x" << block_size_j
- << "; grid=" << grid_single->get_process_grid_rows() << "x" << grid_single->get_process_grid_columns() << std::endl << std::endl;
-
- scalapack_matrix_source = full;
-
- scalapack_matrix_source.copy_to(scalapack_matrix_2d);
- scalapack_matrix_source.copy_to(scalapack_matrix_1d);
- scalapack_matrix_source.copy_to(scalapack_matrix_single);
-
- FullMatrix<NumberType> test_2d(size), test_1d(size), test_one(size);
- scalapack_matrix_2d.copy_to(test_2d);
- scalapack_matrix_1d.copy_to(test_1d);
- scalapack_matrix_single.copy_to(test_one);
- test_2d.add (-1,full);
- test_1d.add (-1,full);
- test_one.add (-1,full);
-
- AssertThrow(test_2d.frobenius_norm() < 1e-12,ExcInternalError());
- AssertThrow(test_1d.frobenius_norm() < 1e-12,ExcInternalError());
- AssertThrow(test_one.frobenius_norm() < 1e-12,ExcInternalError());
+ std::shared_ptr<Utilities::MPI::ProcessGrid> grid_2d = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
+
+ pcout << size << " " << block_size << std::endl;
+
+ // Create s.p.d matrices of requested size:
+ FullMatrix<NumberType> full_A(size);
+ create_spd (full_A);
+
+ //compute eigenpairs of s.p.d matrix
+ ScaLAPACKMatrix<NumberType> scalapack_A_ev (size, grid_2d, block_size);
+ scalapack_A_ev.set_property(LAPACKSupport::Property::symmetric);
+ scalapack_A_ev = full_A;
+ std::vector<NumberType> eigenvalues = scalapack_A_ev.eigenpairs_symmetric(true);
+ FullMatrix<NumberType> eigenvectors (size,size);
+ scalapack_A_ev.copy_to(eigenvectors);
+
+ //compute SVD of s.p.d matrix A = U * SIGMA * VT
+ ScaLAPACKMatrix<NumberType> scalapack_A_sv (size, grid_2d, block_size);
+ ScaLAPACKMatrix<NumberType> scalapack_U (size, grid_2d, block_size);
+ ScaLAPACKMatrix<NumberType> scalapack_VT (size, grid_2d, block_size);
+ scalapack_A_sv.set_property(LAPACKSupport::Property::symmetric);
+ scalapack_A_sv = full_A;
+ std::vector<NumberType> singular_values = scalapack_A_sv.compute_SVD(scalapack_U,scalapack_VT,true,true);
+ FullMatrix<NumberType> l_singular_vectors (size,size);
+ FullMatrix<NumberType> r_singular_vectors (size,size);
+ scalapack_U.copy_to(l_singular_vectors);
+ scalapack_VT.copy_to(r_singular_vectors);
+
+ const unsigned int max_num_values=5;
+ pcout << "comparing the SVD and Eigendecomposition of a s.p.d matrix" << std::endl;
+ for (unsigned i=0; i<max_num_values; ++i)
+ AssertThrow(std::abs(eigenvalues[size-1-i]-singular_values[i])<tol,ExcMessage("singular and eigenvalues do not match"));
+ pcout << " with respect to the given tolerance the singular and eigenvalues coincide" << std::endl;
+
+ Vector<NumberType> eigenvector(size), l_singular_vector(size), r_singular_vector(size);
+ for (unsigned int i=0; i<max_num_values; ++i)
+ {
+ for (unsigned int j=0; j<size; ++j)
+ {
+ eigenvector[j] = eigenvectors(j,size-1-i);
+ l_singular_vector[j] = l_singular_vectors(j,i);
+ r_singular_vector[j] = r_singular_vectors(i,j);
+ }
+ NumberType product_1 = eigenvector*l_singular_vector;
+ NumberType product_2 = eigenvector*r_singular_vector;
+ //the tolerance is reduced for the singular vectors
+ AssertThrow((std::abs(product_1)-1)<tol*10,ExcMessage("left singular vectors and eigenvectors do not coincide"));
+ AssertThrow((std::abs(product_2)-1)<tol*10,ExcMessage("right singular vectors and eigenvectors do not coincide"));
+ }
+ pcout << " with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors" << std::endl;
+ pcout << std::endl;
}
int main (int argc,char **argv)
{
- Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
- const std::vector<unsigned int> blocks_i = {{16,32,64}};
- const std::vector<unsigned int> blocks_j = {{16,32,64}};
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
- for (const auto &s : blocks_i)
- for (const auto &b : blocks_j)
- test<float>(s,b);
+ const std::vector<unsigned int> sizes = {{200,400,600}};
+ const std::vector<unsigned int> blocks = {{32,64}};
+ const double tol_double = 1e-10;
- for (const auto &s : blocks_i)
- for (const auto &b : blocks_j)
- test<double>(s,b);
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<double>(s,b,tol_double);
}
-2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
-1D grid matrix: 500x500; blocks=16x16; grid=1x1
-single process matrix: 500x500; blocks=16x16; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
-1D grid matrix: 500x500; blocks=32x32; grid=1x1
-single process matrix: 500x500; blocks=16x32; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
-1D grid matrix: 500x500; blocks=64x64; grid=1x1
-single process matrix: 500x500; blocks=16x64; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
-1D grid matrix: 500x500; blocks=16x16; grid=1x1
-single process matrix: 500x500; blocks=32x16; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
-1D grid matrix: 500x500; blocks=32x32; grid=1x1
-single process matrix: 500x500; blocks=32x32; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
-1D grid matrix: 500x500; blocks=64x64; grid=1x1
-single process matrix: 500x500; blocks=32x64; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
-1D grid matrix: 500x500; blocks=16x16; grid=1x1
-single process matrix: 500x500; blocks=64x16; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
-1D grid matrix: 500x500; blocks=32x32; grid=1x1
-single process matrix: 500x500; blocks=64x32; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
-1D grid matrix: 500x500; blocks=64x64; grid=1x1
-single process matrix: 500x500; blocks=64x64; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
-1D grid matrix: 500x500; blocks=16x16; grid=1x1
-single process matrix: 500x500; blocks=16x16; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
-1D grid matrix: 500x500; blocks=32x32; grid=1x1
-single process matrix: 500x500; blocks=16x32; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=16x16; grid=1x1
-1D grid matrix: 500x500; blocks=64x64; grid=1x1
-single process matrix: 500x500; blocks=16x64; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
-1D grid matrix: 500x500; blocks=16x16; grid=1x1
-single process matrix: 500x500; blocks=32x16; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
-1D grid matrix: 500x500; blocks=32x32; grid=1x1
-single process matrix: 500x500; blocks=32x32; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=32x32; grid=1x1
-1D grid matrix: 500x500; blocks=64x64; grid=1x1
-single process matrix: 500x500; blocks=32x64; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
-1D grid matrix: 500x500; blocks=16x16; grid=1x1
-single process matrix: 500x500; blocks=64x16; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
-1D grid matrix: 500x500; blocks=32x32; grid=1x1
-single process matrix: 500x500; blocks=64x32; grid=1x1
-
-2D grid matrix: dim=500x500; blocks=64x64; grid=1x1
-1D grid matrix: 500x500; blocks=64x64; grid=1x1
-single process matrix: 500x500; blocks=64x64; grid=1x1
+200 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+200 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+400 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+400 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+600 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+600 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
--- /dev/null
+200 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+200 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+400 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+400 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+600 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+600 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
--- /dev/null
+200 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+200 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+400 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+400 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+600 32
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+600 64
+comparing the SVD and Eigendecomposition of a s.p.d matrix
+ with respect to the given tolerance the singular and eigenvalues coincide
+ with respect to the given tolerance the right and left singular vectors coincide with the eigenvectors
+
+++ /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"
-#include "../lapack/create_matrix.h"
-
-// test eigenpairs_symmetric(const bool, const std::pair<int,int>&, const std::pair<NumberType,NumberType>&)
-
-#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/base/process_grid.h>
-
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/lapack_full_matrix.h>
-#include <deal.II/lac/lapack_templates.h>
-#include <deal.II/lac/scalapack.h>
-
-#include <fstream>
-#include <iostream>
-#include <algorithm>
-#include <memory>
-
-
-template <typename NumberType>
-void test(const unsigned int size, const unsigned int block_size, const NumberType tol)
-{
- 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));
-
- std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
-
- pcout << size << " " << block_size << std::endl;
-
- // Create SPD matrices of requested size:
- FullMatrix<NumberType> full_A(size);
- std::vector<NumberType> lapack_A(size*size);
-
- create_spd (full_A);
- for (unsigned int i = 0; i < size; ++i)
- for (unsigned int j = 0; j < size; ++j)
- lapack_A[i*size+j] = full_A(i,j);
-
- std::vector<NumberType> 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<NumberType> work (1);
-
- //by setting lwork to -1 a workspace query for work is done
- //as matrix is symmetric: LDA == size of matrix
- lwork = -1;
- syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
- lwork=work[0];
- work.resize (lwork);
- syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
-
- AssertThrow (info==0, LAPACKSupport::ExcErrorCode("syev", info));
- }
- unsigned int n_eigenvalues = eigenvalues_Lapack.size(), max_n_eigenvalues=5;
-
- std::vector<Vector<NumberType>> s_eigenvectors_ (max_n_eigenvalues,Vector<NumberType>(size));
- for (int i=0; i<max_n_eigenvalues; ++i)
- for (int j=0; j<size; ++j)
- s_eigenvectors_[i][j] = lapack_A[(size-1-i)*size+j];
-
- pcout << "comparing " << max_n_eigenvalues << " eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:" << std::endl;
- std::vector<NumberType> eigenvalues_psyev;
- ScaLAPACKMatrix<NumberType> scalapack_syev (size, grid, block_size);
- scalapack_syev.set_property(LAPACKSupport::Property::symmetric);
- scalapack_syev = full_A;
- eigenvalues_psyev = scalapack_syev.eigenpairs_symmetric(true);
- FullMatrix<NumberType> p_eigenvectors (size,size);
- scalapack_syev.copy_to(p_eigenvectors);
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- {
- AssertThrow ( std::abs(eigenvalues_psyev[n_eigenvalues-i-1]-eigenvalues_Lapack[n_eigenvalues-i-1]) / std::abs(eigenvalues_Lapack[n_eigenvalues-i-1]) < tol,
- dealii::ExcInternalError());
- }
- pcout << " with respect to the given tolerance the eigenvalues coincide" << std::endl;
-
- std::vector<Vector<NumberType>> p_eigenvectors_ (max_n_eigenvalues,Vector<NumberType>(size));
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- for (unsigned int j=0; j<size; ++j)
- p_eigenvectors_[i][j] = p_eigenvectors(j,size-1-i);
-
- //product of eigenvectors computed using Lapack and ScaLapack has to be either 1 or -1
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- {
- NumberType product = p_eigenvectors_[i] * s_eigenvectors_[i];
-
- //the requirement for alignment of the eigenvectors has to be released (primarily for floats)
- AssertThrow (std::abs(std::abs(product)-1) < tol*10, dealii::ExcInternalError());
- }
- pcout << " with respect to the given tolerance also the eigenvectors coincide" << std::endl << std::endl;
-
- pcout << "comparing " << max_n_eigenvalues << " eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyevx:" << std::endl;
- std::vector<NumberType> eigenvalues_psyevx_partial;
- ScaLAPACKMatrix<NumberType> scalapack_syevx_partial (size, grid, block_size);
- scalapack_syevx_partial.set_property(LAPACKSupport::Property::symmetric);
- scalapack_syevx_partial = full_A;
- eigenvalues_psyevx_partial = scalapack_syevx_partial.eigenpairs_symmetric(true, std::make_pair(size-max_n_eigenvalues+1,size));
- scalapack_syevx_partial.copy_to(p_eigenvectors);
- for (unsigned int i=eigenvalues_psyevx_partial.size()-1; i>0; --i)
- {
- AssertThrow ( std::abs(eigenvalues_psyevx_partial[i]-eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i]) /
- std::abs(eigenvalues_Lapack[size-eigenvalues_psyevx_partial.size()+i]) < tol,
- dealii::ExcInternalError());
- }
- pcout << " with respect to the given tolerance the eigenvalues coincide" << std::endl;
-
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- for (unsigned int j=0; j<size; ++j)
- p_eigenvectors_[i][j] = p_eigenvectors(j,max_n_eigenvalues-1-i);
-
- //product of eigenvectors computed using Lapack and ScaLapack has to be either 1 or -1
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- {
- NumberType product = p_eigenvectors_[i] * s_eigenvectors_[i];
-
- //the requirement for alignment of the eigenvectors has to be released (primarily for floats)
- AssertThrow (std::abs(std::abs(product)-1) < tol*10, dealii::ExcInternalError());
- }
- pcout << " with respect to the given tolerance also the eigenvectors coincide" << 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);
-
- const std::vector<unsigned int> sizes = {{200,400,600}};
- const std::vector<unsigned int> blocks = {{32,64}};
-
- const double tol_double = 1e-10;
- const float tol_float = 1e-5;
-
- for (const auto &s : sizes)
- for (const auto &b : blocks)
- if (b <= s)
- test<float>(s,b,tol_float);
-
- for (const auto &s : sizes)
- for (const auto &b : blocks)
- if (b <= s)
- test<double>(s,b,tol_double);
-}
+++ /dev/null
-200 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-411.91<->411.91 396.894<->396.894 395.588<->395.588 392.905<->392.905 391.261<->391.261
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-411.91<->411.91 396.894<->396.894 395.588<->395.588 392.905<->392.905 391.261<->391.261
- with respect to the given tolerance the eigenvalues coincide
-
-
-200 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-412.36<->412.36 397.285<->397.285 395.54<->395.54 392.855<->392.855 391.447<->391.447
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-412.36<->412.36 397.285<->397.285 395.54<->395.54 392.855<->392.855 391.447<->391.447
- with respect to the given tolerance the eigenvalues coincide
-
-
-400 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-825.747<->825.747 797.968<->797.968 795.8<->795.8 793.229<->793.229 790.87<->790.87
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-825.747<->825.747 797.968<->797.968 795.8<->795.8 793.229<->793.229 790.87<->790.87
- with respect to the given tolerance the eigenvalues coincide
-
-
-400 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-825.723<->825.723 798.22<->798.22 795.562<->795.562 793.391<->793.391 791.088<->791.088
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-825.723<->825.723 798.22<->798.22 795.562<->795.562 793.391<->793.391 791.088<->791.088
- with respect to the given tolerance the eigenvalues coincide
-
-
-600 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-1238.43<->1238.43 1197.98<->1197.98 1195.9<->1195.9 1193.53<->1193.53 1191.69<->1191.69
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-1238.43<->1238.43 1197.98<->1197.98 1195.9<->1195.9 1193.53<->1193.53 1191.69<->1191.69
- with respect to the given tolerance the eigenvalues coincide
-
-
-600 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-1238.82<->1238.82 1197.67<->1197.67 1195.39<->1195.39 1192.85<->1192.85 1191.38<->1191.38
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-1238.82<->1238.82 1197.67<->1197.67 1195.39<->1195.39 1192.85<->1192.85 1191.38<->1191.38
- with respect to the given tolerance the eigenvalues coincide
-
-
+++ /dev/null
-200 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-411.91<->411.91 396.894<->396.894 395.588<->395.588 392.905<->392.905 391.261<->391.261
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-411.91<->411.91 396.894<->396.894 395.588<->395.588 392.905<->392.905 391.261<->391.261
- with respect to the given tolerance the eigenvalues coincide
-
-
-200 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-412.36<->412.36 397.285<->397.285 395.54<->395.54 392.855<->392.855 391.447<->391.447
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-412.36<->412.36 397.285<->397.285 395.54<->395.54 392.855<->392.855 391.447<->391.447
- with respect to the given tolerance the eigenvalues coincide
-
-
-400 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-825.747<->825.747 797.968<->797.968 795.8<->795.8 793.229<->793.229 790.87<->790.87
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-825.747<->825.747 797.968<->797.968 795.8<->795.8 793.229<->793.229 790.87<->790.87
- with respect to the given tolerance the eigenvalues coincide
-
-
-400 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-825.723<->825.723 798.22<->798.22 795.562<->795.562 793.391<->793.391 791.088<->791.088
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-825.723<->825.723 798.22<->798.22 795.562<->795.562 793.391<->793.391 791.088<->791.088
- with respect to the given tolerance the eigenvalues coincide
-
-
-600 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-1238.43<->1238.43 1197.98<->1197.98 1195.9<->1195.9 1193.53<->1193.53 1191.69<->1191.69
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-1238.43<->1238.43 1197.98<->1197.98 1195.9<->1195.9 1193.53<->1193.53 1191.69<->1191.69
- with respect to the given tolerance the eigenvalues coincide
-
-
-600 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-1238.82<->1238.82 1197.67<->1197.67 1195.39<->1195.39 1192.85<->1192.85 1191.38<->1191.38
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-1238.82<->1238.82 1197.67<->1197.67 1195.39<->1195.39 1192.85<->1192.85 1191.38<->1191.38
- with respect to the given tolerance the eigenvalues coincide
-
-
+++ /dev/null
-200 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-411.91<->411.91 396.894<->396.894 395.588<->395.588 392.905<->392.905 391.261<->391.261
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-411.91<->411.91 396.894<->396.894 395.588<->395.588 392.905<->392.905 391.261<->391.261
- with respect to the given tolerance the eigenvalues coincide
-
-
-200 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-412.36<->412.36 397.285<->397.285 395.54<->395.54 392.855<->392.855 391.447<->391.447
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-412.36<->412.36 397.285<->397.285 395.54<->395.54 392.855<->392.855 391.447<->391.447
- with respect to the given tolerance the eigenvalues coincide
-
-
-400 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-825.747<->825.747 797.968<->797.968 795.8<->795.8 793.229<->793.229 790.87<->790.87
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-825.747<->825.747 797.968<->797.968 795.8<->795.8 793.229<->793.229 790.87<->790.87
- with respect to the given tolerance the eigenvalues coincide
-
-
-400 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-825.723<->825.723 798.22<->798.22 795.562<->795.562 793.391<->793.391 791.088<->791.088
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-825.723<->825.723 798.22<->798.22 795.562<->795.562 793.391<->793.391 791.088<->791.088
- with respect to the given tolerance the eigenvalues coincide
-
-
-600 32
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-1238.43<->1238.43 1197.98<->1197.98 1195.9<->1195.9 1193.53<->1193.53 1191.69<->1191.69
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-1238.43<->1238.43 1197.98<->1197.98 1195.9<->1195.9 1193.53<->1193.53 1191.69<->1191.69
- with respect to the given tolerance the eigenvalues coincide
-
-
-600 64
-comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK:
-pdsyev
-1238.82<->1238.82 1197.67<->1197.67 1195.39<->1195.39 1192.85<->1192.85 1191.38<->1191.38
- with respect to the given tolerance the eigenvalues coincide
-
-pdsyevx partial
-1238.82<->1238.82 1197.67<->1197.67 1195.39<->1195.39 1192.85<->1192.85 1191.38<->1191.38
- with respect to the given tolerance the eigenvalues coincide
-
-
projects / dealii.git / commitdiff