]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add template wrappers for ScaLAPACK functions and introduce support for float
authorBenjamin Brands <benjamin.brands@fau.de>
Wed, 10 Jan 2018 10:45:04 +0000 (11:45 +0100)
committerBenjamin Brands <benjamin.brands@fau.de>
Wed, 10 Jan 2018 10:45:04 +0000 (11:45 +0100)
28 files changed:
include/deal.II/lac/scalapack.templates.h
source/base/process_grid.cc
source/lac/scalapack.cc
tests/scalapack/scalapack_01.cc
tests/scalapack/scalapack_01.mpirun=1.output
tests/scalapack/scalapack_01.mpirun=11.output
tests/scalapack/scalapack_01.mpirun=4.output
tests/scalapack/scalapack_02.cc
tests/scalapack/scalapack_02.mpirun=1.output
tests/scalapack/scalapack_02.mpirun=10.output
tests/scalapack/scalapack_02.mpirun=4.output
tests/scalapack/scalapack_03.cc
tests/scalapack/scalapack_03.mpirun=1.output
tests/scalapack/scalapack_03.mpirun=10.output
tests/scalapack/scalapack_03.mpirun=4.output
tests/scalapack/scalapack_04.cc
tests/scalapack/scalapack_05.cc
tests/scalapack/scalapack_05.mpirun=1.output
tests/scalapack/scalapack_05.mpirun=11.output
tests/scalapack/scalapack_05.mpirun=4.output
tests/scalapack/scalapack_06.cc
tests/scalapack/scalapack_06.mpirun=1.output
tests/scalapack/scalapack_06.mpirun=11.output
tests/scalapack/scalapack_06.mpirun=4.output
tests/scalapack/scalapack_07.cc
tests/scalapack/scalapack_07.mpirun=1.output
tests/scalapack/scalapack_07.mpirun=11.output
tests/scalapack/scalapack_07.mpirun=4.output

index 7a7d2e4565736f1f5baa46579fc2c8761c637836..fb3ad897d7cfced32d77efe1387dcaa204d49af3 100644 (file)
@@ -47,7 +47,8 @@ extern "C"
    *
    * 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.
@@ -56,7 +57,9 @@ extern "C"
    *
    * 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
@@ -68,19 +71,29 @@ extern "C"
    *
    * 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.
@@ -91,7 +104,8 @@ extern "C"
    * 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.
@@ -103,7 +117,20 @@ extern "C"
    *
    * 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
@@ -111,7 +138,20 @@ extern "C"
    *
    * 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.
@@ -123,7 +163,11 @@ extern "C"
    *
    * 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
@@ -135,8 +179,19 @@ extern "C"
    */
   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
@@ -150,7 +205,17 @@ extern "C"
    */
   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);
 
   /**
@@ -164,10 +229,27 @@ extern "C"
    */
   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);
 
   /**
@@ -179,8 +261,19 @@ extern "C"
   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.
@@ -189,22 +282,30 @@ extern "C"
    *
    * 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
@@ -217,15 +318,38 @@ extern "C"
    * @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);
 
   /**
@@ -235,28 +359,75 @@ extern "C"
    * 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
@@ -267,11 +438,36 @@ extern "C"
    * 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
@@ -280,10 +476,605 @@ extern "C"
    * 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
index 4f05f91c112447a980a3559dbc2c76713b9acca6..8aa2ea9eb4118c3b984b6c6df5c953274f291d2c 100644 (file)
@@ -34,7 +34,8 @@ namespace
    * 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:
@@ -242,6 +243,7 @@ namespace Utilities
 // 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
 
index 7308231001525177ccf6f9f9c668bf14708a389b..87f1d776f0f82e56c9f20999d30156034923d42b 100644 (file)
@@ -215,8 +215,9 @@ void ScaLAPACKMatrix<NumberType>::compute_cholesky_factorization()
     {
       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;
@@ -234,8 +235,8 @@ void ScaLAPACKMatrix<NumberType>::invert()
     {
       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;
 }
@@ -268,16 +269,16 @@ std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenvalues_symmetric()
       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
@@ -326,16 +327,16 @@ std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric()
       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
@@ -374,8 +375,8 @@ NumberType ScaLAPACKMatrix<NumberType>::reciprocal_condition_number(const Number
       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);
@@ -442,7 +443,7 @@ NumberType ScaLAPACKMatrix<NumberType>::norm(const char type) const
                         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;
@@ -450,6 +451,7 @@ NumberType ScaLAPACKMatrix<NumberType>::norm(const char type) const
 
 // instantiations
 template class ScaLAPACKMatrix<double>;
+template class ScaLAPACKMatrix<float>;
 
 
 DEAL_II_NAMESPACE_CLOSE
index f5a99eeec55e9f7b38ff9128ef5bba388dcb3f28..7e75e5f7596a7a02792c387ea1a6ecc38f59f05d 100644 (file)
@@ -61,7 +61,7 @@ void test(const unsigned int size, const unsigned int block_size)
   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)
     {
@@ -85,6 +85,10 @@ int main (int argc,char **argv)
   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);
index f949c6cec348627837c70ed465c3f059b04e85b3..993f2358986665eda795a684642800c5b8556e7a 100644 (file)
@@ -14,3 +14,19 @@ Ok
 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
index 226c7603678388b06ceb079e1cdfac36e9b12801..6a589474fae8fa5de266f51ce8990fb5a77de896 100644 (file)
@@ -14,3 +14,19 @@ 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
index 020bfd4d9a951be0de6b746aa7e5dc754e4d282d..b092cd7941a2b9ae81a458f1ffaaef222df88cae 100644 (file)
@@ -14,3 +14,19 @@ 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
index 311a4120856e845710797fb648f8de568075853f..db6049c4c8edf44a9bcca95b862b0fc6900cc42e 100644 (file)
@@ -33,7 +33,7 @@
 
 
 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));
@@ -64,9 +64,9 @@ void test(const unsigned int size, const unsigned int block_size)
   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;
@@ -90,8 +90,16 @@ int main (int argc,char **argv)
   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);
 }
index 404b40375c086abe3110dbb2a1c1d7c58b593193..21fc377823d49b4197d5010947fee6aaf4131275 100644 (file)
@@ -16,3 +16,21 @@ Ok
 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
index 999b4e2fcb96f89f2c328381a62c96d4219743fd..384c6b7acd4fd0f22664da607ad5d7e8ddd1e037 100644 (file)
@@ -16,3 +16,21 @@ 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
index 639333e12c6a2b335d418afa675640ed3bcd2972..827afa216bef7cdb9d5eaa254e28f107dfccf882 100644 (file)
@@ -16,3 +16,21 @@ 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
index dd55f1f00d58ea17023e1f0d624d7d96e440cfc8..8e6bb32922b239d3eb2cbb918619463c5213f3e8 100644 (file)
@@ -67,10 +67,10 @@ void test(const unsigned int size, const unsigned int block_size)
   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
@@ -79,7 +79,7 @@ void test(const unsigned int size, const unsigned int block_size)
   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)
     {
@@ -111,6 +111,11 @@ int main (int argc,char **argv)
   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)
index 404b40375c086abe3110dbb2a1c1d7c58b593193..21fc377823d49b4197d5010947fee6aaf4131275 100644 (file)
@@ -16,3 +16,21 @@ Ok
 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
index 999b4e2fcb96f89f2c328381a62c96d4219743fd..384c6b7acd4fd0f22664da607ad5d7e8ddd1e037 100644 (file)
@@ -16,3 +16,21 @@ 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
index 639333e12c6a2b335d418afa675640ed3bcd2972..827afa216bef7cdb9d5eaa254e28f107dfccf882 100644 (file)
@@ -16,3 +16,21 @@ 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
index adcfe259f0e125098fbae401d84ad65e5ed2a555..761e64d0754dec51fe7cb4d3460945fb6fbd5500 100644 (file)
@@ -32,7 +32,7 @@
 #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));
@@ -51,28 +51,28 @@ void test(const unsigned int size, const unsigned int block_size)
 
   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());
 }
 
 
@@ -85,9 +85,16 @@ int main (int argc,char **argv)
   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);
 }
index 97b059aef69a6600a195983ebe1f3d1dcff53d58..b9105209eb0350ab21132d5857e6403262c1981d 100644 (file)
@@ -52,14 +52,14 @@ void test(const unsigned int size, const unsigned int block_size)
   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;
@@ -77,6 +77,11 @@ int main (int argc,char **argv)
   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);
 }
index cb7466c9e3ef63d48d226304261b80250c945e0e..199da4b6bb6196881b234b5cb8d881bd4e83b852 100644 (file)
 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
index f9e7a0569f101357c59d1673025f439c12bc1628..41d09fcab0a71c77caad4b1d01a01fb91ea15d33 100644 (file)
 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
index 673831fda67d430dc94e2ab33938464ca826d0db..2cd3be59cebce17e31c1af14f168d4acf9b96ef4 100644 (file)
 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
index 78e7e09128700175cce5a11a699beb642063925a..249aa3eadfb99289c739a4912723f10ef739b35c 100644 (file)
 
 #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));
@@ -71,15 +67,15 @@ void test(const unsigned int size, const unsigned int block_size)
     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));
   }
@@ -88,20 +84,10 @@ void test(const unsigned int size, const unsigned int block_size)
   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;
 }
 
 
@@ -114,9 +100,16 @@ int main (int argc,char **argv)
   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);
 }
index b6faabf3eb6abf85a660b843368fb9b05ae0bc8a..724105f262d476ef7b95ecd35ef6b9d517bc43c7 100644 (file)
@@ -1,54 +1,72 @@
 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
 
index 5addc81fe7e4f13f3ff59bbacef69fd0238e4860..95a7fc9423a77adb0c7d1a0285ebf5816147ae8e 100644 (file)
@@ -1,54 +1,72 @@
 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
 
index fba5ffd4adffd6da2017e748a5db089ebdc2c0e4..28fc898f2ea23fad8fd9ac23b6ae57f1727d2203 100644 (file)
@@ -1,54 +1,72 @@
 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
 
index d639eff0290a1d64d612f7509a5b949970eb9ec2..4dbaf97d035aa9346d88f914f5eed6f81244cf3d 100644 (file)
 #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));
@@ -73,15 +69,15 @@ void test(const unsigned int size, const unsigned int block_size)
     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));
   }
@@ -92,38 +88,36 @@ void test(const unsigned int size, const unsigned int block_size)
   scalapack_A.copy_to(p_eigenvectors);
   unsigned int n_eigenvalues = eigenvalues_ScaLapack.size(), max_n_eigenvalues=5;
 
-  pcout << "First " << max_n_eigenvalues << " ScaLapack eigenvalues" << std::endl;
-
+  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;
 }
 
 
@@ -136,8 +130,16 @@ int main (int argc,char **argv)
   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);
 }
index b6faabf3eb6abf85a660b843368fb9b05ae0bc8a..7159ec7532c450c68db3fb829fbdb4c84c2dda12 100644 (file)
@@ -1,54 +1,90 @@
 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
 
index 5addc81fe7e4f13f3ff59bbacef69fd0238e4860..ef4d2498816d57bdfe3bd0696ab0411c8280ef60 100644 (file)
@@ -1,54 +1,90 @@
 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
 
index fba5ffd4adffd6da2017e748a5db089ebdc2c0e4..6bf2dcd729f62241e578cb68d51a5a76dd82fd67 100644 (file)
@@ -1,54 +1,90 @@
 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
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.