*
* https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbpnf.htm
*/
- void Cblacs_pinfo(int *rank, int *nprocs);
+ void Cblacs_pinfo(int *rank,
+ int *nprocs);
/**
* Return internal BLACS value in @p val based on the input @p what and @p icontxt.
*
* https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbget.htm
*/
- void Cblacs_get(int icontxt, int what, int *val);
+ void Cblacs_get(int icontxt,
+ int what,
+ int *val);
/**
* Map the processes sequentially in row-major or column-major order
*
* https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbint.htm
*/
- void Cblacs_gridinit(int *context, const char *order, int grid_height, int grid_width);
+ void Cblacs_gridinit(int *context,
+ const char *order,
+ int grid_height,
+ int grid_width);
/**
* Return the process row and column index.
*
* https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbinfo.htm
*/
- void Cblacs_gridinfo(int context, int *grid_height, int *grid_width, int *grid_row, int *grid_col);
+ void Cblacs_gridinfo(int context,
+ int *grid_height,
+ int *grid_width,
+ int *grid_row,
+ int *grid_col);
/**
* Given the system process number, return the row and column coordinates in the BLACS' process grid.
*/
- void Cblacs_pcoord(int ictxt, int pnum, int *prow, int *pcol);
+ void Cblacs_pcoord(int ictxt,
+ int pnum,
+ int *prow,
+ int *pcol);
/**
* Release a BLACS context.
* This routines holds up execution of all processes within the indicated
* scope until they have all called the routine.
*/
- void Cblacs_barrier(int, const char *);
+ void Cblacs_barrier(int,
+ const char *);
/**
* Free all BLACS contexts and releases all allocated memory.
*
* https://software.intel.com/en-us/mkl-developer-reference-c-gerv2d
*/
- void Cdgerv2d(int context, int M, int N, double *A, int lda, int rsrc, int csrc);
+ void Cdgerv2d(int context,
+ int M,
+ int N,
+ double *A,
+ int lda,
+ int rsrc,
+ int csrc);
+ void Csgerv2d(int context,
+ int M,
+ int N,
+ float *A,
+ int lda,
+ int rsrc,
+ int csrc);
/**
* Sends the general rectangular matrix A to the destination
*
* https://software.intel.com/en-us/mkl-developer-reference-c-2018-beta-gesd2d
*/
- void Cdgesd2d(int context , int M, int N, double *A, int lda, int rdest, int cdest);
+ void Cdgesd2d(int context,
+ int M,
+ int N,
+ double *A,
+ int lda,
+ int rdest,
+ int cdest);
+ void Csgesd2d(int context,
+ int M,
+ int N,
+ float *A,
+ int lda,
+ int rdest,
+ int cdest);
/**
* Get BLACS context from MPI @p comm.
*
* https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dnumy.htm
*/
- int numroc_ (const int *n, const int *nb, const int *iproc, const int *isproc, const int *nprocs);
+ int numroc_ (const int *n,
+ const int *nb,
+ const int *iproc,
+ const int *isproc,
+ const int *nprocs);
/**
* Compute the Cholesky factorization of an N-by-N real
*/
void pdpotrf_(const char *UPLO,
const int *N,
- double *A, const int *IA, const int *JA, const int *DESCA,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
int *INFO);
+ void pspotrf_(const char *UPLO,
+ const int *N,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO);
+
/**
* Compute the inverse of a real symmetric positive definite
*/
void pdpotri_(const char *UPLO,
const int *N,
- double *A, const int *IA, const int *JA, const int *DESCA,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO);
+ void pspotri_(const char *UPLO,
+ const int *N,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
int *INFO);
/**
*/
void pdpocon_(const char *uplo,
const int *N,
- const double *A, const int *IA, const int *JA, const int *DESCA,
+ const double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
const double *ANORM, double *RCOND,
- double *WORK, const int *LWORK,
- int *IWORK, const int *LIWORK,
+ double *WORK,
+ const int *LWORK,
+ int *IWORK,
+ const int *LIWORK,
+ int *INFO);
+ void pspocon_(const char *uplo,
+ const int *N,
+ const float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const float *ANORM, float *RCOND,
+ float *WORK,
+ const int *LWORK,
+ int *IWORK,
+ const int *LIWORK,
int *INFO);
/**
double pdlansy_(const char *norm,
const char *uplo,
const int *N,
- const double *A, const int *IA, const int *JA, const int *DESCA,
+ const double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
double *work);
+ float pslansy_(const char *norm,
+ const char *uplo,
+ const int *N,
+ const float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ float *work);
/**
* Compute the Least Common Multiple (LCM) of two positive integers @p M and @p N.
*
* http://www.netlib.org/scalapack/explore-html/d0/d9b/ilcm_8f_source.html
*/
- int ilcm_(const int *M, const int *N);
+ int ilcm_(const int *M,
+ const int *N);
/**
* Return the ceiling of the division of two integers.
*
* http://www.netlib.org/scalapack/explore-html/df/d07/iceil_8f_source.html
*/
- int iceil_(const int *i1, const int *i2);
+ int iceil_(const int *i1,
+ const int *i2);
/**
* Initialize the descriptor vector with the 8 input arguments
*/
void descinit_ (int *desc,
- const int *m, const int *n, const int *mb, const int *nb,
- const int *irsrc, const int *icsrc,
- const int *ictxt, const int *lld, int *info);
+ const int *m,
+ const int *n,
+ const int *mb,
+ const int *nb,
+ const int *irsrc,
+ const int *icsrc,
+ const int *ictxt,
+ const int *lld,
+ int *info);
/**
* Compute the global index of a distributed matrix entry
* @param isrcproc The coordinate of the process that possesses the first row/column of the distributed matrix
* @param nprocs The total number processes over which the distributed matrix is distributed
*/
- int indxl2g_ (const int *indxloc, const int *nb, const int *iproc, const int *isrcproc, const int *nprocs);
+ int indxl2g_ (const int *indxloc,
+ const int *nb,
+ const int *iproc,
+ const int *isrcproc,
+ const int *nprocs);
/**
* Compute the solution to a real system of linear equations
*/
- void pdgesv_(const int *n, const int *nrhs,
- double *A, const int *ia, const int *ja, const int *desca,
+ void pdgesv_(const int *n,
+ const int *nrhs,
+ double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ int *ipiv,
+ double *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ int *info);
+ void psgesv_(const int *n,
+ const int *nrhs,
+ float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
int *ipiv,
- double *B, const int *ib, const int *jb, const int *descb,
+ float *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
int *info);
/**
* sub( C ) denotes C(IC:IC+M-1,JC:JC+N-1), and, op( X ) is one of
* op( X ) = X or op( X ) = X'.
*/
- void pdgemm_(const char *transa, const char *transb,
- const int *m, const int *n, const int *k,
+ void pdgemm_(const char *transa,
+ const char *transb,
+ const int *m,
+ const int *n,
+ const int *k,
const double *alpha,
- double *A, const int *IA, const int *JA, const int *DESCA,
- double *B, const int *IB, const int *JB, const int *DESCB,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ double *B,
+ const int *IB,
+ const int *JB,
+ const int *DESCB,
const double *beta,
- double *C, const int *IC, const int *JC, const int *DESCC);
+ double *C,
+ const int *IC,
+ const int *JC,
+ const int *DESCC);
+ void psgemm_(const char *transa,
+ const char *transb,
+ const int *m,
+ const int *n,
+ const int *k,
+ const float *alpha,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ float *B,
+ const int *IB,
+ const int *JB,
+ const int *DESCB,
+ const float *beta,
+ float *C,
+ const int *IC,
+ const int *JC,
+ const int *DESCC);
/**
* Return the value of the one norm, or the Frobenius norm, or the infinity norm,
* or the element of largest absolute value of a distributed matrix
*/
double pdlange_(char const *norm,
- int const &m, int const &n,
- double *A, int const &ia, int const &ja, int *desca,
+ int const &m,
+ int const &n,
+ double *A,
+ int const &ia,
+ int const &ja,
+ int *desca,
double *work);
+ double pslange_(char const *norm,
+ int const &m,
+ int const &n,
+ float *A,
+ int const &ia,
+ int const &ja,
+ int *desca,
+ float *work);
/**
* Compute the process coordinate which possesses the entry of a
* distributed matrix specified by a global index
*/
- int indxg2p_(const int *glob, const int *nb, const int *iproc, const int *isproc, const int *nprocs);
+ int indxg2p_(const int *glob,
+ const int *nb,
+ const int *iproc,
+ const int *isproc,
+ const int *nprocs);
/**
* Compute all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A
* http://www.netlib.org/scalapack/explore-html/d0/d1a/pdsyev_8f.html
* https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lsyev.htm#lsyev
*/
- void pdsyev_(const char *jobz, const char *uplo,
- const int *m, double *A, const int *ia, const int *ja, int *desca,
+ void pdsyev_(const char *jobz,
+ const char *uplo,
+ const int *m,
+ double *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
double *w,
- double *z, const int *iz, const int *jz, int *descz,
- double *work, const int *lwork, int *info);
+ double *z,
+ const int *iz,
+ const int *jz,
+ int *descz,
+ double *work,
+ const int *lwork,
+ int *info);
+ void pssyev_(const char *jobz,
+ const char *uplo,
+ const int *m,
+ float *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ float *w,
+ float *z,
+ const int *iz,
+ const int *jz,
+ int *descz,
+ float *work,
+ const int *lwork,
+ int *info);
/**
* Copy all or a part of a distributed matrix A to another
* A(ia:ia+m-1,ja:ja+n-1) and sub(B) denotes B(ib:ib+m-1,jb:jb+n-1)
*/
void pdlacpy_(const char *uplo,
- const int *m, const int *n, double *A, const int *ia, const int *ja, int *desca,
- double *B, const int *ib, const int *jb, int *descb);
+ const int *m,
+ const int *n,
+ double *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ double *B,
+ const int *ib,
+ const int *jb,
+ int *descb);
+ void pslacpy_(const char *uplo,
+ const int *m,
+ const int *n,
+ float *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ float *B,
+ const int *ib,
+ const int *jb,
+ int *descb);
+
+ /**
+ * Copies the content of a general rectangular distributed matrix @p A to another distributed matrix @p B
+ * It is not required that the matrices A and B have the same process grid or block size, e.g. copying a
+ * matrix from a one-dimensional to a two-dimensional process grid
+ * @p ictxt is a context which is at least a union of all processes in context A and B
+ */
+ void pdgemr2d_(const int *m,
+ const int *n,
+ const double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ double *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ const int *ictxt);
+ void psgemr2d_(const int *m,
+ const int *n,
+ const float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ float *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ const int *ictxt);
+}
+
+
+/*
+ * 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 Cgerv2d(int context,
+ int M,
+ int N,
+ number *A,
+ int lda,
+ int rsrc,
+ int csrc)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void Cgerv2d(int context,
+ int M,
+ int N,
+ double *A,
+ int lda,
+ int rsrc,
+ int csrc)
+{
+ Cdgerv2d (context, M, N, A, lda, rsrc, csrc);
+}
+
+inline void Cgerv2d(int context,
+ int M,
+ int N,
+ float *A,
+ int lda,
+ int rsrc,
+ int csrc)
+{
+ Csgerv2d (context, M, N, A, lda, rsrc, csrc);
+}
+
+
+template <typename number>
+inline void Cgesd2d(int context,
+ int M,
+ int N,
+ number *A,
+ int lda,
+ int rdest,
+ int cdest)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void Cgesd2d (int context,
+ int M,
+ int N,
+ double *A,
+ int lda,
+ int rdest,
+ int cdest)
+{
+ Cdgesd2d (context, M, N, A, lda, rdest, cdest);
+}
+
+inline void Cgesd2d (int context,
+ int M,
+ int N,
+ float *A,
+ int lda,
+ int rdest,
+ int cdest)
+{
+ Csgesd2d (context, M, N, A, lda, rdest, cdest);
+}
+
+
+template <typename number>
+inline void ppotrf(const char *UPLO,
+ const int *N,
+ number *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void ppotrf(const char *UPLO,
+ const int *N,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO)
+{
+ pdpotrf_(UPLO, N, A, IA, JA, DESCA, INFO);
+}
+
+inline void ppotrf(const char *UPLO,
+ const int *N,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO)
+{
+ pspotrf_(UPLO, N, A, IA, JA, DESCA, INFO);
+}
+
+
+template <typename number>
+inline void ppotri(const char *UPLO,
+ const int *N,
+ number *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void ppotri(const char *UPLO,
+ const int *N,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO)
+{
+ pdpotri_(UPLO, N, A, IA, JA, DESCA, INFO);
+}
+
+inline void ppotri(const char *UPLO,
+ const int *N,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ int *INFO)
+{
+ pspotri_(UPLO, N, A, IA, JA, DESCA, INFO);
+}
+
+
+template <typename number>
+inline void ppocon(const char *uplo,
+ const int *N,
+ const number *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const number *ANORM,
+ number *RCOND,
+ number *WORK,
+ const int *LWORK,
+ int *IWORK,
+ const int *LIWORK,
+ int *INFO)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void ppocon(const char *uplo,
+ const int *N,
+ const double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const double *ANORM,
+ double *RCOND,
+ double *WORK,
+ const int *LWORK,
+ int *IWORK,
+ const int *LIWORK,
+ int *INFO)
+{
+ pdpocon_(uplo, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, IWORK, LIWORK, INFO);
+}
+
+inline void ppocon(const char *uplo,
+ const int *N,
+ const float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ const float *ANORM,
+ float *RCOND,
+ float *WORK,
+ const int *LWORK,
+ int *IWORK,
+ const int *LIWORK,
+ int *INFO)
+{
+ pspocon_(uplo, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, IWORK, LIWORK, INFO);
+}
+
+
+template <typename number>
+inline number plansy(const char *norm,
+ const char *uplo,
+ const int *N,
+ const number *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ number *work)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline double plansy(const char *norm,
+ const char *uplo,
+ const int *N,
+ const double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ double *work)
+{
+ return pdlansy_(norm, uplo, N, A, IA, JA, DESCA, work);
+}
+
+inline float plansy(const char *norm,
+ const char *uplo,
+ const int *N,
+ const float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ float *work)
+{
+ return pslansy_(norm, uplo, N, A, IA, JA, DESCA, work);
+}
+
+
+template <typename number>
+inline void pgesv(const int *n,
+ const int *nrhs,
+ number *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ int *ipiv,
+ number *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ int *info)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgesv(const int *n,
+ const int *nrhs,
+ double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ int *ipiv,
+ double *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ int *info)
+{
+ pdgesv_(n, nrhs, A, ia, ja, desca, ipiv, B, ib, jb, descb, info);
+}
+
+inline void pgesv(const int *n,
+ const int *nrhs,
+ float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ int *ipiv,
+ float *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ int *info)
+{
+ psgesv_(n, nrhs, A, ia, ja, desca, ipiv, B, ib, jb, descb, info);
+}
+
+
+template <typename number>
+inline void pgemm(const char *transa,
+ const char *transb,
+ const int *m,
+ const int *n,
+ const int *k,
+ const number *alpha,
+ number *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ number *B,
+ const int *IB,
+ const int *JB,
+ const int *DESCB,
+ const number *beta,
+ number *C,
+ const int *IC,
+ const int *JC,
+ const int *DESCC)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgemm(const char *transa,
+ const char *transb,
+ const int *m,
+ const int *n,
+ const int *k,
+ const double *alpha,
+ double *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ double *B,
+ const int *IB,
+ const int *JB,
+ const int *DESCB,
+ const double *beta,
+ double *C,
+ const int *IC,
+ const int *JC,
+ const int *DESCC)
+{
+ pdgemm_(transa, transb, m, n, k, alpha, A, IA, JA, DESCA, B, IB, JB, DESCB, beta, C, IC, JC, DESCC);
+}
+
+inline void pgemm(const char *transa,
+ const char *transb,
+ const int *m,
+ const int *n,
+ const int *k,
+ const float *alpha,
+ float *A,
+ const int *IA,
+ const int *JA,
+ const int *DESCA,
+ float *B,
+ const int *IB,
+ const int *JB,
+ const int *DESCB,
+ const float *beta,
+ float *C,
+ const int *IC,
+ const int *JC,
+ const int *DESCC)
+{
+ psgemm_(transa, transb, m, n, k, alpha, A, IA, JA, DESCA, B, IB, JB, DESCB, beta, C, IC, JC, DESCC);
+}
+
+
+template <typename number>
+inline number plange(char const *norm,
+ int const &m,
+ int const &n,
+ number *A,
+ int const &ia,
+ int const &ja,
+ int *desca,
+ number *work)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline double plange(char const *norm,
+ int const &m,
+ int const &n,
+ double *A,
+ int const &ia,
+ int const &ja,
+ int *desca,
+ double *work)
+{
+ return pdlange_(norm, m, n, A, ia, ja,desca, work);
+}
+
+inline float plange(char const *norm,
+ int const &m,
+ int const &n,
+ float *A,
+ int const &ia,
+ int const &ja,
+ int *desca,
+ float *work)
+{
+ return pslange_(norm, m, n, A, ia, ja,desca, work);
}
+
+template <typename number>
+inline void psyev(const char *jobz,
+ const char *uplo,
+ const int *m,
+ number *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ number *w,
+ number *z,
+ const int *iz,
+ const int *jz,
+ int *descz,
+ number *work,
+ const int *lwork,
+ int *info)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void psyev(const char *jobz,
+ const char *uplo,
+ const int *m,
+ double *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ double *w,
+ double *z,
+ const int *iz,
+ const int *jz,
+ int *descz,
+ double *work,
+ const int *lwork,
+ int *info)
+{
+ pdsyev_(jobz, uplo, m, A, ia, ja, desca, w, z, iz, jz, descz, work, lwork, info);
+}
+
+inline void psyev(const char *jobz,
+ const char *uplo,
+ const int *m,
+ float *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ float *w,
+ float *z,
+ const int *iz,
+ const int *jz,
+ int *descz,
+ float *work,
+ const int *lwork,
+ int *info)
+{
+ pssyev_(jobz, uplo, m, A, ia, ja, desca, w, z, iz, jz, descz, work, lwork, info);
+}
+
+
+template <typename number>
+inline void placpy(const char *uplo,
+ const int *m,
+ const int *n,
+ number *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ number *B,
+ const int *ib,
+ const int *jb,
+ int *descb)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void placpy(const char *uplo,
+ const int *m,
+ const int *n,
+ double *A,
+ const int *ia,
+ const int *ja,
+ int *desca, double *B,
+ const int *ib,
+ const int *jb,
+ int *descb)
+{
+ pdlacpy_(uplo, m, n, A, ia, ja, desca, B, ib, jb, descb);
+}
+
+inline void placpy(const char *uplo,
+ const int *m,
+ const int *n,
+ float *A,
+ const int *ia,
+ const int *ja,
+ int *desca,
+ float *B,
+ const int *ib,
+ const int *jb,
+ int *descb)
+{
+ pslacpy_(uplo, m, n, A, ia, ja, desca, B, ib, jb, descb);
+}
+
+
+template <typename number>
+inline void pgemr2d(const int *m,
+ const int *n,
+ const number *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ number *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ const int *ictxt)
+{
+ Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgemr2d(const int *m,
+ const int *n,
+ const double *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ double *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ const int *ictxt)
+{
+ pdgemr2d_(m,n,A,ia,ja,desca,B,ib,jb,descb,ictxt);
+}
+
+inline void pgemr2d(const int *m,
+ const int *n,
+ const float *A,
+ const int *ia,
+ const int *ja,
+ const int *desca,
+ float *B,
+ const int *ib,
+ const int *jb,
+ const int *descb,
+ const int *ictxt)
+{
+ psgemr2d_(m,n,A,ia,ja,desca,B,ib,jb,descb,ictxt);
+}
+
+
#endif // DEAL_II_WITH_SCALAPACK
#endif // dealii_scalapack_templates_h
* https://github.com/elemental/Elemental/blob/master/src/core/Grid.cpp#L67-L91
*/
inline
- std::pair<int,int> compute_processor_grid_sizes(MPI_Comm mpi_comm, const unsigned int m, const unsigned int n,
+ std::pair<int,int> compute_processor_grid_sizes(MPI_Comm mpi_comm,
+ const unsigned int m, const unsigned int n,
const unsigned int block_size_m, const unsigned int block_size_n)
{
// Few notes from the ScaLAPACK user guide:
// instantiations
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;
DEAL_II_NAMESPACE_CLOSE
{
int info = 0;
NumberType *A_loc = &this->values[0];
- pdpotrf_(&uplo,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,&info);
- AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pdpotrf", info));
+ //pdpotrf_(&uplo,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,&info);
+ ppotrf(&uplo,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,&info);
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotrf", info));
}
property = (uplo=='L' ? LAPACKSupport::lower_triangular : LAPACKSupport::upper_triangular);
state = LAPACKSupport::cholesky;
{
int info = 0;
NumberType *A_loc = &this->values[0];
- pdpotri_ (&uplo,&n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,&info);
- AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pdpotri", info));
+ ppotri (&uplo,&n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,&info);
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotri", info));
}
state = LAPACKSupport::inverse_matrix;
}
NumberType *Z_loc = &Z.values[0];
work.resize(1);
- pdsyev_(&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);
+ 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);
- pdsyev_(&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);
+ 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("pdsyev", info));
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyev", info));
}
/*
* send the eigenvalues to processors not being part of the process grid
NumberType *eigenvectors_loc = &eigenvectors.values[0];
work.resize(1);
- pdsyev_(&jobz, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, &ev[0],
- eigenvectors_loc, &eigenvectors.submatrix_row, &eigenvectors.submatrix_column, eigenvectors.descriptor, &work[0], &lwork, &info);
+ 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);
lwork=work[0];
work.resize (lwork);
- pdsyev_(&jobz, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, &ev[0],
- eigenvectors_loc, &eigenvectors.submatrix_row, &eigenvectors.submatrix_column, eigenvectors.descriptor, &work[0], &lwork, &info);
+ 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("pdsyev", info));
+ AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyev", info));
// copy eigenvectors to original matrix
// as the temporary matrix eigenvectors has identical dimensions and
iwork.resize(liwork);
int info = 0;
const NumberType *A_loc = &this->values[0];
- pdpocon_(&uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,
- &a_norm, &rcond, &work[0], &lwork, &iwork[0], &liwork, &info);
+ ppocon(&uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,
+ &a_norm, &rcond, &work[0], &lwork, &iwork[0], &liwork, &info);
AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pdpocon", info));
}
grid->send_to_inactive(&rcond);
2*Nq0+Np0+ldw;
work.resize(lwork);
const NumberType *A_loc = &this->values[0];
- res = pdlansy_(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &work[0]);
+ res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &work[0]);
}
grid->send_to_inactive(&res);
return res;
// instantiations
template class ScaLAPACKMatrix<double>;
+template class ScaLAPACKMatrix<float>;
DEAL_II_NAMESPACE_CLOSE
diff.add( 1., full_in);
diff.add(-1., full_out);
- const double error = diff.frobenius_norm();
+ const NumberType error = diff.frobenius_norm();
if (error > 1e-10 && this_mpi_process == 0)
{
const std::vector<unsigned int> sizes = {{64,120,320,640}};
const std::vector<unsigned int> blocks = {{32,64}};
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ test<float>(s,b);
+
for (const auto &s : sizes)
for (const auto &b : blocks)
test<double>(s,b);
Ok
640 64 1 1
Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
Ok
640 64 3 3
Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
Ok
640 64 2 2
Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok
template <typename NumberType>
-void test(const unsigned int size, const unsigned int block_size)
+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));
diff.add( 1., inverse);
diff.add(-1., full_out);
- const double error = diff.frobenius_norm();
+ const NumberType error = diff.frobenius_norm() / inverse.frobenius_norm();
- if (error > 1e-10 && this_mpi_process == 0)
+ if (error > tol && this_mpi_process == 0)
{
std::cout << "Error!" << std::endl
<< "expected to have:" << std::endl;
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;
+
+ 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);
+ test<double>(s,b,tol_double);
}
Ok
640 64 1 1
Ok
+32 32 1 1
+Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
Ok
640 64 3 3
Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
Ok
640 64 2 2
Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok
inverse.invert(full_in);
inverse.mmult(prod1, full_in);
prod1.add(-1., one);
- const double lapack_error = prod1.linfty_norm();
+ const NumberType lapack_error = prod1.linfty_norm();
// estimated condition number from 1-norm:
- const double k = full_in.l1_norm() * inverse.l1_norm();
+ const NumberType k = full_in.l1_norm() * inverse.l1_norm();
const NumberType tol = k * 1000 * std::numeric_limits<NumberType>::epsilon();
// invert via ScaLAPACK
scalapack_matrix.copy_to(full_out);
full_out.mmult(prod2, full_in);
prod2.add(-1., one);
- const double error = prod2.linfty_norm();
+ const NumberType error = prod2.linfty_norm();
if ( error > tol && this_mpi_process == 0)
{
const std::vector<unsigned int> sizes = {{32,64,120,320,640}};
const std::vector<unsigned int> blocks = {{32,64}};
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<float>(s,b);
+
for (const auto &s : sizes)
for (const auto &b : blocks)
if (b <= s)
Ok
640 64 1 1
Ok
+32 32 1 1
+Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
Ok
640 64 3 3
Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
Ok
640 64 2 2
Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok
#include <iostream>
template <typename NumberType>
-void test(const unsigned int size, const unsigned int block_size)
+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));
create_spd (full_A);
- const double l1 = full_A.l1_norm();
- const double linfty = full_A.linfty_norm();
- const double frobenius = full_A.frobenius_norm();
+ const NumberType l1 = full_A.l1_norm();
+ const NumberType linfty = full_A.linfty_norm();
+ const NumberType frobenius = full_A.frobenius_norm();
scalapack_A = full_A;
// local result on this core:
- const double s_l1 = scalapack_A.l1_norm();
- const double s_linfty = scalapack_A.linfty_norm();
- const double s_frobenius = scalapack_A.frobenius_norm();
+ const NumberType s_l1 = scalapack_A.l1_norm();
+ const NumberType s_linfty = scalapack_A.linfty_norm();
+ const NumberType s_frobenius = scalapack_A.frobenius_norm();
// make sure we have the same result on all cores, do average:
- const double as_l1 = dealii::Utilities::MPI::sum(s_l1, mpi_communicator)/n_mpi_processes;
- const double as_linfty = dealii::Utilities::MPI::sum(s_linfty, mpi_communicator)/n_mpi_processes;
- const double as_frobenius = dealii::Utilities::MPI::sum(s_frobenius, mpi_communicator)/n_mpi_processes;
+ const NumberType as_l1 = dealii::Utilities::MPI::sum(s_l1, mpi_communicator)/n_mpi_processes;
+ const NumberType as_linfty = dealii::Utilities::MPI::sum(s_linfty, mpi_communicator)/n_mpi_processes;
+ const NumberType as_frobenius = dealii::Utilities::MPI::sum(s_frobenius, mpi_communicator)/n_mpi_processes;
pcout << l1 << " " << s_l1 << " " << as_l1 << std::endl
<< linfty << " " << s_linfty << " " << as_linfty << std::endl
<< frobenius << " " << s_frobenius << " " << as_frobenius << std::endl;
- AssertThrow (std::abs(l1 -as_l1) < 1e-10*std::abs(l1), dealii::ExcInternalError());
- AssertThrow (std::abs(linfty -as_linfty) < 1e-10*std::abs(linfty), dealii::ExcInternalError());
- AssertThrow (std::abs(frobenius -as_frobenius) < 1e-10*std::abs(frobenius), dealii::ExcInternalError());
+ AssertThrow (std::abs(l1 -as_l1) < tol*std::abs(l1), dealii::ExcInternalError());
+ AssertThrow (std::abs(linfty -as_linfty) < tol*std::abs(linfty), dealii::ExcInternalError());
+ AssertThrow (std::abs(frobenius -as_frobenius) < tol*std::abs(frobenius), dealii::ExcInternalError());
}
const std::vector<unsigned int> sizes = {{32,64,120,320,640}};
const std::vector<unsigned int> blocks = {{32,64}};
- for (const auto &s : sizes)
+ 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<double>(s,b);
+ 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);
}
create_spd (full_A);
inv_A.invert(full_A);
- const double l1 = full_A.l1_norm();
- const double inv_l1 = inv_A.l1_norm();
+ const NumberType l1 = full_A.l1_norm();
+ const NumberType inv_l1 = inv_A.l1_norm();
// Scalapack:
scalapack_A = full_A;
- const double scalapack_l1 = scalapack_A.l1_norm();
+ const NumberType scalapack_l1 = scalapack_A.l1_norm();
scalapack_A.compute_cholesky_factorization();
- const double rcond = scalapack_A.reciprocal_condition_number(scalapack_l1);
+ const NumberType rcond = scalapack_A.reciprocal_condition_number(scalapack_l1);
pcout << 1./(l1*inv_l1) << " " << rcond << std::endl;
for (const auto &s : sizes)
for (const auto &b : blocks)
if (b <= s)
- test<double>(s,b);
+ test<float>(s,b);
+
+ for (const auto &s : sizes)
+ for (const auto &b : blocks)
+ if (b <= s)
+ test<double>(s,b);
}
0.472941 0.487731
640 64 1 1
0.471001 0.48693
+32 32 1 1
+0.443368 0.511715
+64 32 1 1
+0.441083 0.474693
+64 64 1 1
+0.446963 0.485188
+120 32 1 1
+0.445718 0.476956
+120 64 1 1
+0.46023 0.484744
+320 32 1 1
+0.46245 0.484004
+320 64 1 1
+0.460946 0.48346
+640 32 1 1
+0.472221 0.492506
+640 64 1 1
+0.470315 0.486349
0.472941 0.487731
640 64 3 3
0.471001 0.48693
+32 32 1 1
+0.443368 0.511715
+64 32 2 2
+0.441083 0.474693
+64 64 1 1
+0.446963 0.485188
+120 32 3 3
+0.445718 0.476956
+120 64 2 2
+0.46023 0.484744
+320 32 3 3
+0.46245 0.484004
+320 64 3 3
+0.460946 0.48346
+640 32 3 3
+0.472221 0.492506
+640 64 3 3
+0.470315 0.486349
0.472941 0.487731
640 64 2 2
0.471001 0.48693
+32 32 1 1
+0.443368 0.511715
+64 32 2 2
+0.441083 0.474693
+64 64 1 1
+0.446963 0.485188
+120 32 2 2
+0.445718 0.476956
+120 64 2 2
+0.46023 0.484744
+320 32 2 2
+0.46245 0.484004
+320 64 2 2
+0.460946 0.48346
+640 32 2 2
+0.472221 0.492506
+640 64 2 2
+0.470315 0.486349
#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>
-extern "C" //Some Lapack routines
-{
- void dsyev_(char *jobz, char *uplo, int *n, double *A, int *lda, double *w, double *work, int *lwork, int *info);
- void ssyev_(char *jobz, char *uplo, int *n, float *A, int *lda, float *w, float *work, int *lwork, int *info);
-}
+
template <typename NumberType>
-void test(const unsigned int size, const unsigned int block_size)
+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));
char uplo = 'U'; //storage format of the matrix A; not so important as matrix is symmetric
int LDA = size; //leading dimension of the matrix A
int lwork; //length of vector/array work
- std::vector<double> work (1);
+ 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;
- dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
+ syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
lwork=work[0];
work.resize (lwork);
- dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
+ syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
AssertThrow (info==0, LAPACKSupport::ExcErrorCode("syev", info));
}
eigenvalues_ScaLapack = scalapack_A.eigenvalues_symmetric();
unsigned int n_eigenvalues = eigenvalues_ScaLapack.size(), max_n_eigenvalues=5;
- pcout << "First " << max_n_eigenvalues << " ScaLapack eigenvalues" << std::endl;
-
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- pcout << eigenvalues_ScaLapack[i] << " ";
-
- pcout << std::endl << "First " << max_n_eigenvalues << " Lapack eigenvalues" << std::endl;
-
+ pcout << "comparing " << max_n_eigenvalues << " eigenvalues computed using LAPACK and ScaLAPACK:" << std::endl;
for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- pcout << eigenvalues_Lapack[i] << " ";
-
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- AssertThrow ( std::abs(eigenvalues_ScaLapack[i]-eigenvalues_Lapack[i]) < std::abs(eigenvalues_Lapack[i])*1e-10, dealii::ExcInternalError());
-
- pcout << std::endl << std::endl;
+ 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;
}
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;
+
for (const auto &s : sizes)
for (const auto &b : blocks)
if (b <= s)
- test<double>(s,b);
+ 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);
}
32 32 1 1
-First 5 ScaLapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
-First 5 Lapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
64 32 1 1
-First 5 ScaLapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
-First 5 Lapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
64 64 1 1
-First 5 ScaLapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
-First 5 Lapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
120 32 1 1
-First 5 ScaLapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
-First 5 Lapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
120 64 1 1
-First 5 ScaLapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
-First 5 Lapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
320 32 1 1
-First 5 ScaLapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
-First 5 Lapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
320 64 1 1
-First 5 ScaLapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
-First 5 Lapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
640 32 1 1
-First 5 ScaLapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
-First 5 Lapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
640 64 1 1
-First 5 ScaLapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
-First 5 Lapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
+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:
+ with respect to the given tolerance the eigenvalues coincide
+
+64 32 1 1
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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:
+ with respect to the given tolerance the eigenvalues coincide
+
+120 64 1 1
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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:
+ with respect to the given tolerance the eigenvalues coincide
+
+640 32 1 1
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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
-First 5 ScaLapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
-First 5 Lapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
64 32 2 2
-First 5 ScaLapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
-First 5 Lapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
64 64 1 1
-First 5 ScaLapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
-First 5 Lapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
120 32 3 3
-First 5 ScaLapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
-First 5 Lapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
120 64 2 2
-First 5 ScaLapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
-First 5 Lapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
320 32 3 3
-First 5 ScaLapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
-First 5 Lapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
320 64 3 3
-First 5 ScaLapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
-First 5 Lapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
640 32 3 3
-First 5 ScaLapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
-First 5 Lapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
640 64 3 3
-First 5 ScaLapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
-First 5 Lapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
+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:
+ with respect to the given tolerance the eigenvalues coincide
+
+64 32 2 2
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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:
+ with respect to the given tolerance the eigenvalues coincide
+
+120 64 2 2
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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:
+ with respect to the given tolerance the eigenvalues coincide
+
+640 32 3 3
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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
-First 5 ScaLapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
-First 5 Lapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
64 32 2 2
-First 5 ScaLapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
-First 5 Lapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
64 64 1 1
-First 5 ScaLapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
-First 5 Lapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
120 32 2 2
-First 5 ScaLapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
-First 5 Lapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
120 64 2 2
-First 5 ScaLapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
-First 5 Lapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
320 32 2 2
-First 5 ScaLapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
-First 5 Lapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
320 64 2 2
-First 5 ScaLapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
-First 5 Lapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
640 32 2 2
-First 5 ScaLapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
-First 5 Lapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
640 64 2 2
-First 5 ScaLapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
-First 5 Lapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
+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:
+ with respect to the given tolerance the eigenvalues coincide
+
+64 32 2 2
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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:
+ with respect to the given tolerance the eigenvalues coincide
+
+120 64 2 2
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues 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:
+ with respect to the given tolerance the eigenvalues coincide
+
+640 32 2 2
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
+
+640 64 2 2
+comparing 5 eigenvalues computed using LAPACK and ScaLAPACK:
+ with respect to the given tolerance the eigenvalues coincide
#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 <algorithm>
#include <memory>
-extern "C" //Some Lapack routines
-{
- void dsyev_(char *jobz, char *uplo, int *n, double *A, int *lda, double *w, double *work, int *lwork, int *info);
- void ssyev_(char *jobz, char *uplo, int *n, float *A, int *lda, float *w, float *work, int *lwork, int *info);
-}
template <typename NumberType>
-void test(const unsigned int size, const unsigned int block_size)
+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));
char uplo = 'U'; //storage format of the matrix A; not so important as matrix is symmetric
int LDA = size; //leading dimension of the matrix A
int lwork; //length of vector/array work
- std::vector<double> work (1);
+ 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;
- dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
+ syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
lwork=work[0];
work.resize (lwork);
- dsyev_(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
+ syev(&jobz, &uplo, &LDA, & *lapack_A.begin(), &LDA, & *eigenvalues_Lapack.begin(), & *work.begin(), &lwork, &info);
AssertThrow (info==0, LAPACKSupport::ExcErrorCode("syev", info));
}
scalapack_A.copy_to(p_eigenvectors);
unsigned int n_eigenvalues = eigenvalues_ScaLapack.size(), max_n_eigenvalues=5;
- pcout << "First " << max_n_eigenvalues << " ScaLapack eigenvalues" << std::endl;
-
+ 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)
- pcout << eigenvalues_ScaLapack[i] << " ";
-
- pcout << std::endl << "First " << max_n_eigenvalues << " Lapack eigenvalues" << std::endl;
-
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- pcout << eigenvalues_Lapack[i] << " ";
- pcout << std::endl;
- pcout << std::endl;
-
- for (unsigned int i=0; i<max_n_eigenvalues; ++i)
- AssertThrow ( std::fabs(eigenvalues_ScaLapack[i]-eigenvalues_Lapack[i]) < std::fabs(eigenvalues_Lapack[i])*1e-10, dealii::ExcInternalError());
+ 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<size; ++i)
+ for (unsigned int i=0; i<max_n_eigenvalues; ++i)
{
- Vector<double> p_ev(size), s_ev(size);
- for (unsigned int j=0; j<size; ++j)
- {
- p_ev[j] = p_eigenvectors(j,i);
- s_ev[j] = s_eigenvectors(i,j);
- }
- double product = p_ev * s_ev;
- AssertThrow (std::fabs(std::fabs(product)-1) < 1e-6, dealii::ExcInternalError());
+ 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 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;
+
+ 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);
+ test<double>(s,b,tol_double);
}
32 32 1 1
-First 5 ScaLapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
-First 5 Lapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
+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
-First 5 ScaLapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
-First 5 Lapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
+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
-First 5 ScaLapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
-First 5 Lapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
+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
-First 5 ScaLapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
-First 5 Lapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
+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
-First 5 ScaLapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
-First 5 Lapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
+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
-First 5 ScaLapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
-First 5 Lapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
+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
-First 5 ScaLapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
-First 5 Lapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
+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
-First 5 ScaLapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
-First 5 Lapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
+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
-First 5 ScaLapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
-First 5 Lapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
+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
32 32 1 1
-First 5 ScaLapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
-First 5 Lapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
+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
-First 5 ScaLapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
-First 5 Lapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
+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
-First 5 ScaLapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
-First 5 Lapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
+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
-First 5 ScaLapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
-First 5 Lapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
+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
-First 5 ScaLapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
-First 5 Lapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
+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 3 3
-First 5 ScaLapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
-First 5 Lapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
+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
-First 5 ScaLapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
-First 5 Lapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
+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 3 3
-First 5 ScaLapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
-First 5 Lapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
+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 3 3
-First 5 ScaLapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
-First 5 Lapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
+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 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
+
+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 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:
+ 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:
+ 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
-First 5 ScaLapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
-First 5 Lapack eigenvalues
-29.4327 29.6255 29.7881 30.092 30.3062
+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
-First 5 ScaLapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
-First 5 Lapack eigenvalues
-59.5716 59.9842 60.1172 60.3363 60.4051
+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
-First 5 ScaLapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
-First 5 Lapack eigenvalues
-59.7161 59.8977 60.3065 60.5131 60.6934
+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
-First 5 ScaLapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
-First 5 Lapack eigenvalues
-114.098 114.163 114.52 114.662 114.877
+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
-First 5 ScaLapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
-First 5 Lapack eigenvalues
-113.953 114.207 114.496 114.525 114.8
+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
-First 5 ScaLapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
-First 5 Lapack eigenvalues
-309.886 310.093 310.385 310.648 310.719
+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
-First 5 ScaLapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
-First 5 Lapack eigenvalues
-309.771 310.096 310.171 310.37 310.471
+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
-First 5 ScaLapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
-First 5 Lapack eigenvalues
-625.62 625.745 625.974 625.996 626.13
+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
-First 5 ScaLapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
-First 5 Lapack eigenvalues
-625.381 625.696 625.855 625.983 626.074
+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