<< "file for information on how to ensure that deal.II "
<< "picks up an existing BLAS and LAPACK installation at "
<< "configuration time.");
-
}
*/
void copy_to (FullMatrix<NumberType> &matrix) const;
-
/**
* Copy the contents of the distributed matrix into a differently distributed matrix @p dest.
* The function also works for matrices with different process grids
*/
void copy_to (ScaLAPACKMatrix<NumberType> &dest) const;
-
/**
* Compute the Cholesky factorization of the matrix using ScaLAPACK
* function <code>pXpotrf</code>. The result of the factorization is stored in this object.
*/
void invert();
-
-
/**
- * 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.
+ * Computing selected eigenvalues and, optionally, the eigenvectors of the real symmetric
+ * matrix $A \in \mathbb{R}^{M \times M}$.
+ *
+ * The eigenvalues/eigenvectors are selected by prescribing a range of indices @p index_limits.
+ *
* 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.
+ *
+ * If all eigenvalues/eigenvectors have to be computed, pass the closed interval $ \left[ 0, M-1 \right] $ in @p index_limits.
+ *
+ * Pass the closed interval $ \left[ M-r, M-1 \right] $ if the $r$ largest eigenvalues/eigenvectors are desired.
*/
- 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));
-
+ std::vector<NumberType> eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &index_limits,
+ const bool compute_eigenvectors);
+ /**
+ * Computing selected eigenvalues and, optionally, the eigenvectors.
+ * The eigenvalues/eigenvectors are selected by prescribing a range of values @p value_limits for the eigenvalues.
+ *
+ * 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_by_value(const std::pair<NumberType,NumberType> &value_limits,
+ const bool compute_eigenvectors);
/**
- * 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
+ * Computing the singular value decomposition (SVD) of a
+ * matrix $A \in \mathbb{R}^{M \times N}$, optionally computing the left and/or right
+ * singular vectors. The SVD is written as $A = U * \Sigma * V^T$
+ * with $\Sigma \in \mathbb{R}^{M \times N}$ as a diagonal matrix,
+ * $U \in \mathbb{R}^{M \times M}$ and $U \in \mathbb{R}^{M \times M}$
+ * as orthogonal matrices. 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.
+ * singular values are returned in decreasing order and only the first $\min(M,N)$
+ * columns of $U$ and rows of VT = $V^T$ 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
+ * 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
+ * If right singular vectors are required matrices $A$ and $V^T$
* have to be constructed with the same process grid and block cyclic distribution.
+ *
+ * To avoid computing the left and/or right singular vectors the function accepts <code>nullptr</code>
+ * for @p U and/or @p VT.
*/
- std::vector<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> &U,
- ScaLAPACKMatrix<NumberType> &VT,
- const bool left_singluar_vectors=false,
- const bool right_singluar_vectors=false);
-
+ std::vector<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> *U = nullptr,
+ ScaLAPACKMatrix<NumberType> *VT = nullptr);
/**
- * Function solves overdetermined or underdetermined real linear
- * systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A.
+ * Solving overdetermined or underdetermined real linear
+ * systems involving matrix $A \in \mathbb{R}^{M \times N}$, or its transpose $A^T$,
+ * using a QR or LQ factorization of $A$ for $N_{\rm RHS}$ RHS vectors in the columns of matrix $B$
*
- * It is assumed that A has full rank: \f$rank(A) = \min(M,N)\f$.
- * Upon exit the columns of B contain the solutions and
- * the following options are supported:
- * - 1. If transpose==false and \f$M \geq N\f$: least squares solution of overdetermined system
- * \f$\min \Vert B - A X\Vert\f$.
+ * It is assumed that $A$ has full rank: $rank(A) = \min(M,N)$.
*
- * Upon exit the rows 0 to N-1 contain the least square solution vectors. The residual sum of squares
+ * The following options are supported:
+ * - 1. If <code>transpose==false</code> and $M \geq N$: least squares solution of overdetermined system
+ * $\min \Vert B - A*X\Vert$.
+ *
+ * Upon exit the rows 0 to N-1 of $B$ contain the least square solution vectors. The residual sum of squares
* for each column is given by the sum of squares of elements N to M-1 in that column
*
- * - 2. If transpose==false and \f$M < N\f$: find minimum norm solutions of underdetermined systems
- * \f$A X = B\f$.
+ * - 2. If <code>transpose==false</code> and $M < N$: find minimum norm solutions of underdetermined systems
+ * $A * X = B$.
*
- * Upon exit the columns of B contain the minimum norm solution vectors
+ * Upon exit the columns of $B$ contain the minimum norm solution vectors
*
- * - 3. If transpose==true and and \f$M \geq N\f$: find minimum norm solutions of underdetermined system
- * \f$ A^\top X = B\f$
+ * - 3. If <code>transpose==true</code> and $M \geq N$: find minimum norm solutions of underdetermined system
+ * $ A^\top X = B$.
*
- * Upon exit the columns of B contain the minimum norm solution vectors
+ * Upon exit the columns of $B$ contain the minimum norm solution vectors.
*
- * - 4. If transpose==true and \f$M < N\f$: least squares solution of overdetermined system
- * \f$\min \Vert B - A^\top X\Vert\f$.
+ * - 4. If <code>transpose==true</code> and $M < N$: least squares solution of overdetermined system
+ * $\min \Vert B - A^\top X\Vert$.
*
* Upon exit the rows 0 to M-1 contain the least square solution vectors. The residual sum of squares
* for each column is given by the sum of squares of elements M to N-1 in that column
* .
- * If transpose==false B is M x NRHS matrix, otherwise it is NxNRHS.
- * The matrices A and B must have an identical block cyclic distribution for rows and columns
+ * If <code>transpose==false</code> then $B \in \mathbb{R}^{M \times N_{\rm RHS}}$,
+ * otherwise $B \in \mathbb{R}^{N \times N_{\rm RHS}}}$.
+ * The matrices $A$ and $B$ must have an identical block cyclic distribution for rows and columns.
*/
void least_squares(ScaLAPACKMatrix<NumberType> &B,
const bool transpose=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()).
*/
NumberType norm(const char type) const;
+ /**
+ * Computing selected eigenvalues and, optionally, the 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(const bool compute_eigenvectors,
+ const std::pair<unsigned int,unsigned int> &index_limits=
+ std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int),
+ const std::pair<NumberType,NumberType> &value_limits=
+ std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),std::numeric_limits<NumberType>::quiet_NaN()));
+
/**
* Since ScaLAPACK operations notoriously change the meaning of the matrix
* entries, we record the current state after the last operation here.
* 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,
//process is active in the process grid
if (this->grid->mpi_process_is_active)
dest.values = this->values;
+
+ dest.state = state;
+ dest.property = property;
}
+template <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &index_limits,
+ const bool compute_eigenvectors)
+{
+ // check validity of index limits
+ Assert (index_limits.first < (unsigned int)n_rows,ExcIndexRange(index_limits.first,0,n_rows));
+ Assert (index_limits.second < (unsigned int)n_rows,ExcIndexRange(index_limits.second,0,n_rows));
+
+ std::pair<unsigned int,unsigned int> idx = std::make_pair(std::min(index_limits.first,index_limits.second),
+ std::max(index_limits.first,index_limits.second));
+
+ // compute all eigenvalues/eigenvectors
+ if (idx.first==0 && idx.second==(unsigned int)n_rows-1)
+ return eigenpairs_symmetric(compute_eigenvectors);
+ else
+ return eigenpairs_symmetric(compute_eigenvectors,idx);
+}
+
+
+
+template <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_by_value(const std::pair<NumberType,NumberType> &value_limits,
+ const bool compute_eigenvectors)
+{
+ Assert (!std::isnan(value_limits.first),ExcMessage("value_limits.first is NaN"));
+ Assert (!std::isnan(value_limits.second),ExcMessage("value_limits.second is NaN"));
+
+ std::pair<unsigned int,unsigned int> indices = std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int);
+
+ return eigenpairs_symmetric(compute_eigenvectors,indices,value_limits);
+}
+
+
+
template <typename NumberType>
std::vector<NumberType>
ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvectors,
- const std::pair<int,int> &eigenvalue_idx,
+ const std::pair<unsigned int, unsigned int> &eigenvalue_idx,
const std::pair<NumberType,NumberType> &eigenvalue_limits)
{
Assert (state == LAPACKSupport::matrix,
Threads::Mutex::ScopedLock lock (mutex);
+ const bool use_values = (std::isnan(eigenvalue_limits.first) || std::isnan(eigenvalue_limits.second)) ? false : true;
+ const bool use_indices = ((eigenvalue_idx.first==numbers::invalid_unsigned_int) || (eigenvalue_idx.second==numbers::invalid_unsigned_int)) ? false : true;
+
+ Assert(!(use_values && use_indices),ExcMessage("Prescribing both the index and value range for the eigenvalues is ambiguous"));
+
// 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);
+ 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;
// 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> gap(n_local_rows * n_local_columns);
// index range for eigenvalues is not specified
- if (eigenvalue_idx.first==-1 && eigenvalue_idx.second==-1)
+ if (!use_indices)
{
// 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)
+ if (!use_values)
{
range = 'A';
all_eigenpairs = true;
}
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);
+ //as Fortran starts counting/indexing from 1 unlike C/C++, where it starts from 0
+ il = std::min(eigenvalue_idx.first,eigenvalue_idx.second) + 1;
+ iu = std::max(eigenvalue_idx.first,eigenvalue_idx.second) + 1;
}
NumberType *A_loc = &this->values[0];
/*
*/
int lwork=-1;
int liwork=-1;
- NumberType *eigenvectors_loc = (compute_eigenvectors ? &eigenvectors->values[0] : NULL);
+ NumberType *eigenvectors_loc = (compute_eigenvectors ? &eigenvectors->values[0] : nullptr);
work.resize(1);
iwork.resize (1);
* 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 (compute_eigenvectors)
+ {
+ property = LAPACKSupport::Property::general;
+ state = LAPACKSupport::eigenvalues;
+ }
+ else
+ state = LAPACKSupport::unusable;
return ev;
}
+
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)
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix<NumberType> *U,
+ ScaLAPACKMatrix<NumberType> *VT)
{
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));
+ Assert(row_block_size==column_block_size,ExcDimensionMismatch(row_block_size,column_block_size));
+
+ const bool left_singluar_vectors = (U != nullptr) ? true : false;
+ const bool right_singluar_vectors = (VT != nullptr) ? true : false;
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));
+ 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));
+ 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);
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;
+ NumberType *U_loc = left_singluar_vectors ? &(U->values[0]) : nullptr;
+ NumberType *VT_loc = right_singluar_vectors ? &(VT->values[0]) : nullptr;
int info = 0;
/*
* by setting lwork to -1 a workspace query for optimal length of work is performed
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,
+ & *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));
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,
+ & *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));
}
B_loc,&B.submatrix_row,&B.submatrix_column,B.descriptor,&work[0],&lwork,&info);
AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgels", info));
}
+ state = LAPACKSupport::State::unusable;
}
#include "../tests.h"
#include "../lapack/create_matrix.h"
-// test eigenpairs_symmetric(const bool, const std::pair<int,int>&, const std::pair<NumberType,NumberType>&)
+// test eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &, const bool)
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
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;
+ pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
// Create SPD matrices of requested size:
FullMatrix<NumberType> full_A(size);
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);
+ eigenvalues_psyev = scalapack_syev.eigenpairs_symmetric_by_index(std::make_pair(0,size-1),true);
FullMatrix<NumberType> p_eigenvectors (size,size);
scalapack_syev.copy_to(p_eigenvectors);
for (unsigned int i=0; i<max_n_eigenvalues; ++i)
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));
+ eigenvalues_psyevx_partial = scalapack_syevx_partial.eigenpairs_symmetric_by_index(std::make_pair(size-max_n_eigenvalues,size-1),true);
scalapack_syevx_partial.copy_to(p_eigenvectors);
for (unsigned int i=eigenvalues_psyevx_partial.size()-1; i>0; --i)
{
-200 32
+200 32 1 1
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
with respect to the given tolerance also the eigenvectors coincide
-200 64
+200 64 1 1
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
with respect to the given tolerance also the eigenvectors coincide
-400 32
+400 32 1 1
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
with respect to the given tolerance also the eigenvectors coincide
-400 64
+400 64 1 1
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
with respect to the given tolerance also the eigenvectors coincide
-600 32
+600 32 1 1
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
with respect to the given tolerance also the eigenvectors coincide
-600 64
+600 64 1 1
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
-200 32
+200 32 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-200 64
+200 64 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-400 32
+400 32 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-400 64
+400 64 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-600 32
+600 32 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-600 64
+600 64 3 3
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
-200 32
+200 32 2 2
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
with respect to the given tolerance also the eigenvectors coincide
-200 64
+200 64 2 2
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
with respect to the given tolerance also the eigenvectors coincide
-400 32
+400 32 2 2
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
with respect to the given tolerance also the eigenvectors coincide
-400 64
+400 64 2 2
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
with respect to the given tolerance also the eigenvectors coincide
-600 32
+600 32 2 2
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
with respect to the given tolerance also the eigenvectors coincide
-600 64
+600 64 2 2
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
-200 32
+200 32 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-200 64
+200 64 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-400 32
+400 32 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-400 64
+400 64 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-600 32
+600 32 3 3
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
with respect to the given tolerance also the eigenvectors coincide
-600 64
+600 64 3 3
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
#include "../tests.h"
#include "../lapack/create_matrix.h"
-// test compute_SVD(ScaLAPACKMatrix<NumberType>&,ScaLAPACKMatrix<NumberType>&,const bool,const bool)
+// test compute_SVD(ScaLAPACKMatrix<NumberType>*,ScaLAPACKMatrix<NumberType>*)
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
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);
+ std::vector<NumberType> eigenvalues = scalapack_A_ev.eigenpairs_symmetric_by_index(std::make_pair(0,size-1),true);
FullMatrix<NumberType> eigenvectors (size,size);
scalapack_A_ev.copy_to(eigenvectors);
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);
+ std::vector<NumberType> singular_values = scalapack_A_sv.compute_SVD(&scalapack_U,&scalapack_VT);
FullMatrix<NumberType> l_singular_vectors (size,size);
FullMatrix<NumberType> r_singular_vectors (size,size);
scalapack_U.copy_to(l_singular_vectors);
projects / dealii.git / commitdiff