<< "file for information on how to ensure that deal.II "
<< "picks up an existing BLAS and LAPACK installation at "
<< "configuration time.");
+
}
*/
void invert();
- /**
- * Compute all eigenvalues of a real symmetric matrix using <code>pXsyev</code>.
- * If successful, the computed eigenvalues are arranged in ascending order.
- * After this function is called, the content of the matrix is overwritten
- * making it unusable.
- */
- std::vector<NumberType> eigenvalues_symmetric();
+
/**
- * Compute all eigenpairs of a real symmetric matrix using <code>pXsyev</code>.
+ * Function to compute selected eigenvalues and, optionally, the eigenvectors.
+ * If the function is called with the default arguments all eigenvalues are computed but no eigenvectors.
+ * The eigenvalues/eigenvectors are selected by either prescribing a range of indices @p index_limits
+ * or a range of values @p value_limits for the eigenvalues. The funtion will throw an exception
+ * if both ranges are prescribed (meaning that both ranges differ from the default value)
+ * as this ambiguity is prohibited.
* If successful, the computed eigenvalues are arranged in ascending order.
* The eigenvectors are stored in the columns of the matrix, thereby
* overwriting the original content of the matrix.
*/
- std::vector<NumberType> eigenpairs_symmetric ();
+ std::vector<NumberType> eigenpairs_symmetric(const bool compute_eigenvectors=false,
+ const std::pair<int,int> &index_limits = std::make_pair(-1,-1),
+ const std::pair<NumberType,NumberType> &value_limits = std::make_pair(-1,-1));
/**
* Estimate the the condition number of a SPD matrix in the $l_1$-norm.
const int *jb,
const int *descb,
const int *ictxt);
-}
+
+ /**
+ * helper routines determining machine precision
+ */
+ double pdlamch_(const int *ictxt,
+ const char *cmach);
+ float pslamch_(const int *ictxt,
+ const char *cmach);
+ /**
+ * psyevx computes selected eigenvalues and, optionally, eigenvectors
+ * of a real symmetric matrix A. Eigenvalues/vectors can be selected by
+ * specifying a range of values or a range of indices for the desired
+ * eigenvalues.
+ */
+ void pdsyevx_(const char *jobz,
+ const char *range,
+ const char *uplo,
+ const int *n,
+ double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ const double *VL,
+ const double *VU,
+ const int *il,
+ const int *iu,
+ const double *abstol,
+ const int *m,
+ const int *nz,
+ double *w,
+ double *orfac,
+ double *Z,
+ const int *iz,
+ const int *jz,
+ const int *descz,
+ double *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *ifail,
+ int *iclustr,
+ double *gap,
+ int *info);
+ void pssyevx_(const char *jobz,
+ const char *range,
+ const char *uplo,
+ const int *n,
+ float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ const float *VL,
+ const float *VU,
+ const int *il,
+ const int *iu,
+ const float *abstol,
+ const int *m,
+ const int *nz,
+ float *w,
+ float *orfac,
+ float *Z,
+ const int *iz,
+ const int *jz,
+ const int *descz,
+ float *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *ifail,
+ int *iclustr,
+ float *gap,
+ 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 plamch(const int *ictxt,
+ const char *cmach,
+ number &val)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void plamch(const int *ictxt,
+ const char *cmach,
+ double &val)
+{
+ val = pdlamch_(ictxt,cmach);
+}
+
+inline void plamch(const int *ictxt,
+ const char *cmach,
+ float &val)
+{
+ val = pslamch_(ictxt,cmach);
+}
+
+
+template <typename number>
+inline void psyevx(const char *jobz,
+ const char *range,
+ const char *uplo,
+ const int *n,
+ number *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ number *VL,
+ number *VU,
+ const int *il,
+ const int *iu,
+ number *abstol,
+ const int *m,
+ const int *nz,
+ number *w,
+ number *orfac,
+ number *Z,
+ const int *iz,
+ const int *jz,
+ const int *descz,
+ number *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *ifail,
+ int *iclustr,
+ number *gap,
+ int *info)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void psyevx(const char *jobz,
+ const char *range,
+ const char *uplo,
+ const int *n,
+ double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ double *VL,
+ double *VU,
+ const int *il,
+ const int *iu,
+ double *abstol,
+ const int *m,
+ const int *nz,
+ double *w,
+ double *orfac,
+ double *Z,
+ const int *iz,
+ const int *jz,
+ const int *descz,
+ double *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *ifail,
+ int *iclustr,
+ double *gap,
+ int *info)
+{
+ pdsyevx_(jobz,range,uplo,n,A,ia,ja,desca,VL,VU,il,iu,abstol,m,nz,w,orfac,Z,iz,jz,descz,work,lwork,iwork,liwork,ifail,iclustr,gap,info);
+}
+
+inline void psyevx(const char *jobz,
+ const char *range,
+ const char *uplo,
+ const int *n,
+ float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ float *VL,
+ float *VU,
+ const int *il,
+ const int *iu,
+ float *abstol,
+ const int *m,
+ const int *nz,
+ float *w,
+ float *orfac,
+ float *Z,
+ const int *iz,
+ const int *jz,
+ const int *descz,
+ float *work,
+ int *lwork,
+ int *iwork,
+ int *liwork,
+ int *ifail,
+ int *iclustr,
+ float *gap,
+ int *info)
+{
+ pssyevx_(jobz,range,uplo,n,A,ia,ja,desca,VL,VU,il,iu,abstol,m,nz,w,orfac,Z,iz,jz,descz,work,lwork,iwork,liwork,ifail,iclustr,gap,info);
+}
+
+
#endif // DEAL_II_WITH_SCALAPACK
#endif // dealii_scalapack_templates_h
template <typename NumberType>
-std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenvalues_symmetric()
+std::vector<NumberType>
+ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvectors,
+ const std::pair<int,int> &eigenvalue_idx,
+ const std::pair<NumberType,NumberType> &eigenvalue_limits)
{
- Assert (state == LAPACKSupport::matrix,
- ExcMessage("Matrix has to be in Matrix state before calling this function."));
- Assert (property == LAPACKSupport::symmetric,
- ExcMessage("Matrix has to be symmetric for this operation."));
- Threads::Mutex::ScopedLock lock (mutex);
-
- ScaLAPACKMatrix<NumberType> Z (grid->n_mpi_processes, grid, 1);
- std::vector<NumberType> ev (n_rows);
-
- if (grid->mpi_process_is_active)
- {
- int info = 0;
-
- char jobz = 'N';
- NumberType *A_loc = &this->values[0];
-
- /*
- * by setting lwork to -1 a workspace query for optimal length of work is performed
- */
- int lwork=-1;
- NumberType *Z_loc = &Z.values[0];
- work.resize(1);
-
- psyev(&jobz, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, &ev[0],
- Z_loc, &Z.submatrix_row, &Z.submatrix_column, Z.descriptor, &work[0], &lwork, &info);
-
- lwork=work[0];
- work.resize (lwork);
-
- psyev(&jobz, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, &ev[0],
- Z_loc, &Z.submatrix_row, &Z.submatrix_column, Z.descriptor, &work[0], &lwork, &info);
-
- AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyev", info));
- }
- /*
- * send the eigenvalues to processors not being part of the process grid
- */
- grid->send_to_inactive(ev.data(), ev.size());
-
- /*
- * 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
- */
- state = LAPACKSupport::unusable;
+ //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"));
- return ev;
-}
-
-
-
-template <typename NumberType>
-std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric()
-{
Assert (state == LAPACKSupport::matrix,
ExcMessage("Matrix has to be in Matrix state before calling this function."));
Assert (property == LAPACKSupport::symmetric,
Threads::Mutex::ScopedLock lock (mutex);
- ScaLAPACKMatrix<NumberType> eigenvectors (n_rows, grid, row_block_size);
- eigenvectors.property = property;
+ // if computation of eigenvectors is not required use a sufficiently small distributed matrix
+ std::unique_ptr<ScaLAPACKMatrix<NumberType>> eigenvectors = compute_eigenvectors ?
+ std::make_unique<ScaLAPACKMatrix<NumberType>>(n_rows, grid, row_block_size)
+ :
+ std::make_unique<ScaLAPACKMatrix<NumberType>>(grid->n_process_rows, grid->n_process_columns,
+ grid,1,1);
+
+ //ScaLAPACKMatrix<NumberType> eigenvectors (n_rows, grid, row_block_size);
+ eigenvectors->property = property;
std::vector<NumberType> ev(n_rows);
if (grid->mpi_process_is_active)
{
int info = 0;
-
/*
- * for jobz = 'V' all eigenpairs of the matrix are computed
+ * for jobz==N only eigenvalues are computed, for jobz='V' also the eigenvectors of the matrix are computed
*/
- char jobz = 'V';
- NumberType *A_loc = &this->values[0];
+ char jobz = compute_eigenvectors ? 'V' : 'N';
+ char range;
+ // default value is to compute all eigenvalues and optionally eigenvectors
+ 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;
+ NumberType abstol;
+ char cmach = compute_eigenvectors ? 'U' : 'S';
+
+ // orfac decides which eigenvectors should be reorthogonalized
+ // see http://www.netlib.org/scalapack/explore-html/df/d1a/pdsyevx_8f_source.html for explanation
+ // to keeps simple no reorthogonalized will be done by setting orfac to 0
+ NumberType orfac = 0;
+ //contains the indices of eigenvectors that failed to converge
+ std::vector<int> ifail;
+ // This array contains indices of eigenvectors corresponding to
+ // a cluster of eigenvalues that could not be reorthogonalized
+ // due to insufficient workspace
+ // see http://www.netlib.org/scalapack/explore-html/df/d1a/pdsyevx_8f_source.html for explanation
+ std::vector<int> iclustr;
+ // This array contains the gap between eigenvalues whose
+ // eigenvectors could not be reorthogonalized.
+ // see http://www.netlib.org/scalapack/explore-html/df/d1a/pdsyevx_8f_source.html for explanation
+ std::vector<NumberType> gap(n_local_rows * n_local_columns);
+
+ // index range for eigenvalues is not specified
+ if (eigenvalue_idx.first==-1 && eigenvalue_idx.second==-1)
+ {
+ // interval for eigenvalues is not specified and consequently all eigenvalues/eigenpairs will be computed
+ if (std::abs(eigenvalue_limits.first-eigenvalue_limits.second)<1e-12 && std::abs(eigenvalue_limits.first+1)<1e-12)
+ {
+ range = 'A';
+ all_eigenpairs = true;
+ }
+ else
+ {
+ range = 'V';
+ all_eigenpairs = false;
+ vl = std::min(eigenvalue_limits.first,eigenvalue_limits.second);
+ vu = std::max(eigenvalue_limits.first,eigenvalue_limits.second);
+ }
+ }
+ else
+ {
+ Assert(std::abs(eigenvalue_limits.first-eigenvalue_limits.second)<1e-12 && std::abs(eigenvalue_limits.first+1)<1e-12,
+ ExcMessage("Prescribing both the index and value range for the eigenvalues is ambiguous"));
+ range = 'I';
+ all_eigenpairs = false;
+ il = std::min(eigenvalue_idx.first,eigenvalue_idx.second);
+ iu = std::max(eigenvalue_idx.first,eigenvalue_idx.second);
+ }
+ NumberType *A_loc = &this->values[0];
/*
* by setting lwork to -1 a workspace query for optimal length of work is performed
*/
int lwork=-1;
- NumberType *eigenvectors_loc = &eigenvectors.values[0];
+ int liwork=-1;
+ NumberType *eigenvectors_loc = (compute_eigenvectors ? &eigenvectors->values[0] : NULL);
work.resize(1);
+ iwork.resize (1);
- psyev(&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);
-
+ if (all_eigenpairs)
+ {
+ psyev(&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);
+ }
+ else
+ {
+ plamch( &(this->grid->blacs_context), &cmach, abstol);
+ abstol *= 2;
+ ifail.resize(n_rows);
+ iclustr.resize(n_local_rows * n_local_columns);
+ gap.resize(n_local_rows * n_local_columns);
+
+ psyevx(&jobz, &range, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor,
+ &vl, &vu, &il, &iu, &abstol, &m, &nz, &ev[0], &orfac,
+ eigenvectors_loc, &eigenvectors->submatrix_row, &eigenvectors->submatrix_column, eigenvectors->descriptor,
+ &work[0], &lwork, &iwork[0], &liwork, &ifail[0], &iclustr[0], &gap[0], &info);
+ }
lwork=work[0];
work.resize (lwork);
- psyev(&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);
+ if (all_eigenpairs)
+ {
+ psyev(&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);
- AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyev", info));
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyev", info));
+ }
+ else
+ {
+ liwork = iwork[0];
+ AssertThrow(liwork>0,ExcInternalError());
+ iwork.resize(liwork);
- // copy eigenvectors to original matrix
+ psyevx(&jobz, &range, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor,
+ &vl, &vu, &il, &iu, &abstol, &m, &nz, &ev[0], &orfac,
+ eigenvectors_loc, &eigenvectors->submatrix_row, &eigenvectors->submatrix_column, eigenvectors->descriptor,
+ &work[0], &lwork, &iwork[0], &liwork, &ifail[0], &iclustr[0], &gap[0], &info);
+
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyevx", info));
+ }
+ // if eigenvectors are queried 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);
+ if (compute_eigenvectors)
+ this->values.swap(eigenvectors->values);
+
+ //adapt the size of ev to fit m upon return
+ while ((int)ev.size() > m)
+ ev.pop_back();
}
/*
* send the eigenvalues to processors not being part of the process grid
grid->send_to_inactive(ev.data(), ev.size());
/*
- * On exit matrix A stores the eigenvectors in the columns
+ * if only eigenvalues are queried the content of the matrix will be destroyed
+ * if the eigenpairs are queried matrix A on exit stores the eigenvectors in the columns
*/
property = LAPACKSupport::Property::general;
state = LAPACKSupport::eigenvalues;
if (i==j)
// since A(i,j) < 1 and
// a symmetric diagonally dominant matrix is SPD
- A(i,j) = val + size;
+ A(i,j) = val + size + i*i/size;
else
{
A(i,j) = val;
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;
--- /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