]> https://gitweb.dealii.org/ - dealii.git/commitdiff
reviewer corrections 5734/head
authorBenjamin Brands <benjamin.brands@fau.de>
Wed, 17 Jan 2018 12:09:03 +0000 (13:09 +0100)
committerbenjamin <benjamin.brands@fau.de>
Fri, 19 Jan 2018 19:38:35 +0000 (20:38 +0100)
include/deal.II/lac/lapack_support.h
include/deal.II/lac/scalapack.h
include/deal.II/lac/scalapack.templates.h
source/lac/scalapack.cc
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_06.mpirun=9.output
tests/scalapack/scalapack_08.cc

index 08c4e3e111295f8e2fdfe8b97c5ba73fa411e864..df40f7f5a87e4b8d9d851c1f80f8ad9fcd3869d3 100644 (file)
@@ -214,7 +214,6 @@ namespace LAPACKSupport
                  << "file for information on how to ensure that deal.II "
                  << "picks up an existing BLAS and LAPACK installation at "
                  << "configuration time.");
-
 }
 
 
index 405b1248fd2480ee3e956b677806ba281088245c..4155ab4e14c4f45e5a9bb3248e5d8b10c8dc3c94 100644 (file)
@@ -158,7 +158,6 @@ public:
    */
   void copy_to (FullMatrix<NumberType> &matrix) const;
 
-
   /**
    * Copy the contents of the distributed matrix into a differently distributed matrix @p dest.
    * The function also works for matrices with different process grids
@@ -166,7 +165,6 @@ public:
    */
   void copy_to (ScaLAPACKMatrix<NumberType> &dest) const;
 
-
   /**
    * Compute the Cholesky factorization of the matrix using ScaLAPACK
    * function <code>pXpotrf</code>. The result of the factorization is stored in this object.
@@ -180,85 +178,97 @@ public:
    */
   void invert();
 
-
-
   /**
-   * Function to compute selected eigenvalues and, optionally, the eigenvectors.
-   * If the function is called with the default arguments all eigenvalues are computed but no eigenvectors.
-   * The eigenvalues/eigenvectors are selected by either prescribing a range of indices @p index_limits
-   * or a range of values @p value_limits for the eigenvalues. The funtion will throw an exception
-   * if both ranges are prescribed (meaning that both ranges differ from the default value)
-   * as this ambiguity is prohibited.
+   * Computing selected eigenvalues and, optionally, the eigenvectors of the real symmetric
+   * matrix $A \in \mathbb{R}^{M \times M}$.
+   *
+   * The eigenvalues/eigenvectors are selected by prescribing a range of indices @p index_limits.
+   *
    * If successful, the computed eigenvalues are arranged in ascending order.
    * The eigenvectors are stored in the columns of the matrix, thereby
    * overwriting the original content of the matrix.
+   *
+   * If all eigenvalues/eigenvectors have to be computed, pass the closed interval $ \left[ 0, M-1 \right] $ in @p index_limits.
+   *
+   * Pass the closed interval $ \left[ M-r, M-1 \right] $ if the $r$ largest eigenvalues/eigenvectors are desired.
    */
-  std::vector<NumberType> eigenpairs_symmetric(const bool compute_eigenvectors=false,
-                                               const std::pair<int,int> &index_limits = std::make_pair(-1,-1),
-                                               const std::pair<NumberType,NumberType> &value_limits = std::make_pair(-1,-1));
-
+  std::vector<NumberType> eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &index_limits,
+                                                        const bool compute_eigenvectors);
 
+  /**
+   * Computing selected eigenvalues and, optionally, the eigenvectors.
+   * The eigenvalues/eigenvectors are selected by prescribing a range of values @p value_limits for the eigenvalues.
+   *
+   * If successful, the computed eigenvalues are arranged in ascending order.
+   * The eigenvectors are stored in the columns of the matrix, thereby
+   * overwriting the original content of the matrix.
+   */
+  std::vector<NumberType> eigenpairs_symmetric_by_value(const std::pair<NumberType,NumberType> &value_limits,
+                                                        const bool compute_eigenvectors);
 
   /**
-  * Funcion to compute the singular value decomposition (SVD) of an
-  * M-by-N matrix A, optionally computing the left and/or right
-  * singular vectors. The SVD is written as A = U * SIGMA * transpose(V)
-  * where SIGMA is an M-by-N diagonal matrix, @p U is an M-by-M orthogonal matrix,
-  * and @p V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
-  * are the singular values of A and the columns of U and V are the
+  * Computing the singular value decomposition (SVD) of a
+  * matrix $A \in \mathbb{R}^{M \times N}$, optionally computing the left and/or right
+  * singular vectors. The SVD is written as $A = U * \Sigma * V^T$
+  * with $\Sigma \in \mathbb{R}^{M \times N}$ as a diagonal matrix,
+  * $U \in \mathbb{R}^{M \times M}$ and $U \in \mathbb{R}^{M \times M}$
+  * as orthogonal matrices. The diagonal elements of $\Sigma$
+  * are the singular values of $A$ and the columns of $U$ and $V$ are the
   * corresponding left and right singular vectors, respectively. The
-  * singular values are returned in decreasing order and only the first min(M,N)
-  * columns of U and rows of VT = transpose(V) are computed.
+  * singular values are returned in decreasing order and only the first $\min(M,N)$
+  * columns of $U$ and rows of VT = $V^T$ are computed.
+  *
   * Upon return the content of the matrix is unusable.
-  * The matrix A must have identical block cyclic distribution for the rows and column
-  * If left singular vectors are required matrices A and U
+  * The matrix A must have identical block cyclic distribution for the rows and column.
+  *
+  * If left singular vectors are required matrices $A$ and $U$
   * have to be constructed with the same process grid and block cyclic distribution.
-  * If right singular vectors are required matrices A and VT
+  * If right singular vectors are required matrices $A$ and $V^T$
   * have to be constructed with the same process grid  and block cyclic distribution.
+  *
+  * To avoid computing the left and/or right singular vectors the function accepts <code>nullptr</code>
+  * for @p U and/or @p VT.
   */
-  std::vector<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> &U,
-                                      ScaLAPACKMatrix<NumberType> &VT,
-                                      const bool left_singluar_vectors=false,
-                                      const bool right_singluar_vectors=false);
-
+  std::vector<NumberType> compute_SVD(ScaLAPACKMatrix<NumberType> *U = nullptr,
+                                      ScaLAPACKMatrix<NumberType> *VT = nullptr);
 
   /**
-  * Function solves overdetermined or underdetermined real linear
-  * systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A.
+  * Solving overdetermined or underdetermined real linear
+  * systems involving matrix $A \in \mathbb{R}^{M \times N}$, or its transpose $A^T$,
+  * using a QR or LQ factorization of $A$ for $N_{\rm RHS}$ RHS vectors in the columns of matrix $B$
   *
-  * It is assumed that A has full rank: \f$rank(A) = \min(M,N)\f$.
-  * Upon exit the columns of B contain the solutions and
-  * the following options are supported:
-  * - 1. If transpose==false and \f$M \geq N\f$: least squares solution of overdetermined system
-  *      \f$\min \Vert B - A X\Vert\f$.
+  * It is assumed that $A$ has full rank: $rank(A) = \min(M,N)$.
   *
-  *      Upon exit the rows 0 to N-1 contain the least square solution vectors. The residual sum of squares
+  * The following options are supported:
+  * - 1. If <code>transpose==false</code> and $M \geq N$: least squares solution of overdetermined system
+  *      $\min \Vert B - A*X\Vert$.
+  *
+  *      Upon exit the rows 0 to N-1 of $B$ contain the least square solution vectors. The residual sum of squares
   *      for each column is given by the sum of squares of elements N to M-1 in that column
   *
-  * - 2. If transpose==false and \f$M < N\f$: find minimum norm solutions of underdetermined systems
-  *      \f$A X = B\f$.
+  * - 2. If <code>transpose==false</code> and $M < N$: find minimum norm solutions of underdetermined systems
+  *      $A * X = B$.
   *
-  *      Upon exit the columns of B contain the minimum norm solution vectors
+  *      Upon exit the columns of $B$ contain the minimum norm solution vectors
   *
-  * - 3. If transpose==true and and \f$M \geq N\f$: find minimum norm solutions of underdetermined system
-  *      \f$ A^\top X = B\f$
+  * - 3. If <code>transpose==true</code> and $M \geq N$: find minimum norm solutions of underdetermined system
+  *      $ A^\top X = B$.
   *
-  *      Upon exit the columns of B contain the minimum norm solution vectors
+  *      Upon exit the columns of $B$ contain the minimum norm solution vectors.
   *
-  * - 4. If transpose==true and \f$M < N\f$: least squares solution of overdetermined system
-  *      \f$\min \Vert B - A^\top X\Vert\f$.
+  * - 4. If <code>transpose==true</code> and $M < N$: least squares solution of overdetermined system
+  *      $\min \Vert B - A^\top X\Vert$.
   *
   *      Upon exit the rows 0 to M-1 contain the least square solution vectors. The residual sum of squares
   *      for each column is given by the sum of squares of elements M to N-1 in that column
   * .
-  * If transpose==false B is M x NRHS matrix, otherwise it is NxNRHS.
-  * The matrices A and B must have an identical block cyclic distribution for rows and columns
+  * If <code>transpose==false</code> then $B \in \mathbb{R}^{M \times N_{\rm RHS}}$,
+  * otherwise $B \in \mathbb{R}^{N \times N_{\rm RHS}}}$.
+  * The matrices $A$ and $B$ must have an identical block cyclic distribution for rows and columns.
   */
   void least_squares(ScaLAPACKMatrix<NumberType> &B,
                      const bool transpose=false);
 
-
-
   /**
    * Estimate the the condition number of a SPD matrix in the $l_1$-norm.
    * The matrix has to be in the Cholesky state (see compute_cholesky_factorization()).
@@ -337,6 +347,22 @@ private:
    */
   NumberType norm(const char type) const;
 
+  /**
+   * Computing selected eigenvalues and, optionally, the eigenvectors.
+   * The eigenvalues/eigenvectors are selected by either prescribing a range of indices @p index_limits
+   * or a range of values @p value_limits for the eigenvalues. The funtion will throw an exception
+   * if both ranges are prescribed (meaning that both ranges differ from the default value)
+   * as this ambiguity is prohibited.
+   * If successful, the computed eigenvalues are arranged in ascending order.
+   * The eigenvectors are stored in the columns of the matrix, thereby
+   * overwriting the original content of the matrix.
+   */
+  std::vector<NumberType> eigenpairs_symmetric(const bool compute_eigenvectors,
+                                               const std::pair<unsigned int,unsigned int> &index_limits=
+                                                 std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int),
+                                               const std::pair<NumberType,NumberType> &value_limits=
+                                                 std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),std::numeric_limits<NumberType>::quiet_NaN()));
+
   /**
    * Since ScaLAPACK operations notoriously change the meaning of the matrix
    * entries, we record the current state after the last operation here.
index 7ef24d5dc557417c0d8c7bb0f8c0acd9db7c1b6f..e1ec89df91f69e0a1c1532950ac54025013b2e59 100644 (file)
@@ -606,7 +606,6 @@ extern "C"
    * M-by-N matrix A, optionally computing the left and/or right
    * singular vectors
    */
-  //FIXME
   void pdgesvd_(const char *jobu,
                 const char *jobvt,
                 const int *m,
index 832199c2bb87bc91fa2fae75c11fe94cde064935..15c85da517178d3925ccc63d09d13956c13fadf1 100644 (file)
@@ -280,6 +280,9 @@ ScaLAPACKMatrix<NumberType>::copy_to (ScaLAPACKMatrix<NumberType> &dest) const
     //process is active in the process grid
     if (this->grid->mpi_process_is_active)
       dest.values = this->values;
+
+  dest.state = state;
+  dest.property = property;
 }
 
 
@@ -322,10 +325,44 @@ void ScaLAPACKMatrix<NumberType>::invert()
 
 
 
+template <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &index_limits,
+    const bool compute_eigenvectors)
+{
+  // check validity of index limits
+  Assert (index_limits.first < (unsigned int)n_rows,ExcIndexRange(index_limits.first,0,n_rows));
+  Assert (index_limits.second < (unsigned int)n_rows,ExcIndexRange(index_limits.second,0,n_rows));
+
+  std::pair<unsigned int,unsigned int> idx = std::make_pair(std::min(index_limits.first,index_limits.second),
+                                                            std::max(index_limits.first,index_limits.second));
+
+  // compute all eigenvalues/eigenvectors
+  if (idx.first==0 && idx.second==(unsigned int)n_rows-1)
+    return eigenpairs_symmetric(compute_eigenvectors);
+  else
+    return eigenpairs_symmetric(compute_eigenvectors,idx);
+}
+
+
+
+template <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_by_value(const std::pair<NumberType,NumberType> &value_limits,
+    const bool compute_eigenvectors)
+{
+  Assert (!std::isnan(value_limits.first),ExcMessage("value_limits.first is NaN"));
+  Assert (!std::isnan(value_limits.second),ExcMessage("value_limits.second is NaN"));
+
+  std::pair<unsigned int,unsigned int> indices = std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int);
+
+  return eigenpairs_symmetric(compute_eigenvectors,indices,value_limits);
+}
+
+
+
 template <typename NumberType>
 std::vector<NumberType>
 ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvectors,
-                                                  const std::pair<int,int> &eigenvalue_idx,
+                                                  const std::pair<unsigned int, unsigned int> &eigenvalue_idx,
                                                   const std::pair<NumberType,NumberType> &eigenvalue_limits)
 {
   Assert (state == LAPACKSupport::matrix,
@@ -335,14 +372,16 @@ ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvector
 
   Threads::Mutex::ScopedLock lock (mutex);
 
+  const bool use_values = (std::isnan(eigenvalue_limits.first) || std::isnan(eigenvalue_limits.second)) ? false : true;
+  const bool use_indices = ((eigenvalue_idx.first==numbers::invalid_unsigned_int) || (eigenvalue_idx.second==numbers::invalid_unsigned_int)) ? false : true;
+
+  Assert(!(use_values && use_indices),ExcMessage("Prescribing both the index and value range for the eigenvalues is ambiguous"));
+
   // if computation of eigenvectors is not required use a sufficiently small distributed matrix
   std::unique_ptr<ScaLAPACKMatrix<NumberType>> eigenvectors = compute_eigenvectors ?
-                                                              std::make_unique<ScaLAPACKMatrix<NumberType>>(n_rows, grid, row_block_size)
-                                                              :
-                                                              std::make_unique<ScaLAPACKMatrix<NumberType>>(grid->n_process_rows, grid->n_process_columns,
-                                                                  grid,1,1);
+                                                              std::make_unique<ScaLAPACKMatrix<NumberType>>(n_rows,grid,row_block_size) :
+                                                              std::make_unique<ScaLAPACKMatrix<NumberType>>(grid->n_process_rows,grid->n_process_columns,grid,1,1);
 
-  //ScaLAPACKMatrix<NumberType> eigenvectors (n_rows, grid, row_block_size);
   eigenvectors->property = property;
   // number of eigenvalues to be returned; upon successful exit ev contains the m seclected eigenvalues in ascending order
   int m = n_rows;
@@ -383,10 +422,10 @@ ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvector
       std::vector<NumberType> gap(n_local_rows * n_local_columns);
 
       // index range for eigenvalues is not specified
-      if (eigenvalue_idx.first==-1 && eigenvalue_idx.second==-1)
+      if (!use_indices)
         {
           // interval for eigenvalues is not specified and consequently all eigenvalues/eigenpairs will be computed
-          if (std::abs(eigenvalue_limits.first-eigenvalue_limits.second)<1e-12 && std::abs(eigenvalue_limits.first+1)<1e-12)
+          if (!use_values)
             {
               range = 'A';
               all_eigenpairs = true;
@@ -401,13 +440,11 @@ ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvector
         }
       else
         {
-          Assert(std::abs(eigenvalue_limits.first-eigenvalue_limits.second)<1e-12 && std::abs(eigenvalue_limits.first+1)<1e-12,
-                 ExcMessage("Prescribing both the index and value range for the eigenvalues is ambiguous"));
-
           range = 'I';
           all_eigenpairs = false;
-          il = std::min(eigenvalue_idx.first,eigenvalue_idx.second);
-          iu = std::max(eigenvalue_idx.first,eigenvalue_idx.second);
+          //as Fortran starts counting/indexing from 1 unlike C/C++, where it starts from 0
+          il = std::min(eigenvalue_idx.first,eigenvalue_idx.second) + 1;
+          iu = std::max(eigenvalue_idx.first,eigenvalue_idx.second) + 1;
         }
       NumberType *A_loc = &this->values[0];
       /*
@@ -415,7 +452,7 @@ ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvector
        */
       int lwork=-1;
       int liwork=-1;
-      NumberType *eigenvectors_loc = (compute_eigenvectors ? &eigenvectors->values[0] : NULL);
+      NumberType *eigenvectors_loc = (compute_eigenvectors ? &eigenvectors->values[0] : nullptr);
       work.resize(1);
       iwork.resize (1);
 
@@ -491,38 +528,45 @@ ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvector
    * if only eigenvalues are queried the content of the matrix will be destroyed
    * if the eigenpairs are queried matrix A on exit stores the eigenvectors in the columns
    */
-  property = LAPACKSupport::Property::general;
-  state = LAPACKSupport::eigenvalues;
+  if (compute_eigenvectors)
+    {
+      property = LAPACKSupport::Property::general;
+      state = LAPACKSupport::eigenvalues;
+    }
+  else
+    state = LAPACKSupport::unusable;
 
   return ev;
 }
 
 
+
 template <typename NumberType>
-std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix<NumberType> &U,
-    ScaLAPACKMatrix<NumberType> &VT,
-    const bool left_singluar_vectors,
-    const bool right_singluar_vectors)
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix<NumberType> *U,
+    ScaLAPACKMatrix<NumberType> *VT)
 {
   Assert (state == LAPACKSupport::matrix,
           ExcMessage("Matrix has to be in Matrix state before calling this function."));
-  Assert(row_block_size==column_block_size,ExcDimensionMismatch(n_rows,U.n_rows));
+  Assert(row_block_size==column_block_size,ExcDimensionMismatch(row_block_size,column_block_size));
+
+  const bool left_singluar_vectors = (U != nullptr) ? true : false;
+  const bool right_singluar_vectors = (VT != nullptr) ? true : false;
 
   if (left_singluar_vectors)
     {
-      Assert(n_rows==U.n_rows,ExcDimensionMismatch(n_rows,U.n_rows));
-      Assert(U.n_rows==U.n_columns,ExcDimensionMismatch(U.n_rows,U.n_columns));
-      Assert(row_block_size==U.row_block_size,ExcDimensionMismatch(row_block_size,U.row_block_size));
-      Assert(column_block_size==U.column_block_size,ExcDimensionMismatch(column_block_size,U.column_block_size));
-      Assert(grid->blacs_context==U.grid->blacs_context,ExcDimensionMismatch(grid->blacs_context,U.grid->blacs_context));
+      Assert(n_rows==U->n_rows,ExcDimensionMismatch(n_rows,U->n_rows));
+      Assert(U->n_rows==U->n_columns,ExcDimensionMismatch(U->n_rows,U->n_columns));
+      Assert(row_block_size==U->row_block_size,ExcDimensionMismatch(row_block_size,U->row_block_size));
+      Assert(column_block_size==U->column_block_size,ExcDimensionMismatch(column_block_size,U->column_block_size));
+      Assert(grid->blacs_context==U->grid->blacs_context,ExcDimensionMismatch(grid->blacs_context,U->grid->blacs_context));
     }
   if (right_singluar_vectors)
     {
-      Assert(n_columns==VT.n_rows,ExcDimensionMismatch(n_columns,VT.n_rows));
-      Assert(VT.n_rows==VT.n_columns,ExcDimensionMismatch(VT.n_rows,VT.n_columns));
-      Assert(row_block_size==VT.row_block_size,ExcDimensionMismatch(row_block_size,VT.row_block_size));
-      Assert(column_block_size==VT.column_block_size,ExcDimensionMismatch(column_block_size,VT.column_block_size));
-      Assert(grid->blacs_context==VT.grid->blacs_context,ExcDimensionMismatch(grid->blacs_context,VT.grid->blacs_context));
+      Assert(n_columns==VT->n_rows,ExcDimensionMismatch(n_columns,VT->n_rows));
+      Assert(VT->n_rows==VT->n_columns,ExcDimensionMismatch(VT->n_rows,VT->n_columns));
+      Assert(row_block_size==VT->row_block_size,ExcDimensionMismatch(row_block_size,VT->row_block_size));
+      Assert(column_block_size==VT->column_block_size,ExcDimensionMismatch(column_block_size,VT->column_block_size));
+      Assert(grid->blacs_context==VT->grid->blacs_context,ExcDimensionMismatch(grid->blacs_context,VT->grid->blacs_context));
     }
   Threads::Mutex::ScopedLock lock (mutex);
 
@@ -533,8 +577,8 @@ std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix
       char jobu = left_singluar_vectors ? 'V' : 'N';
       char jobvt = right_singluar_vectors ? 'V' : 'N';
       NumberType *A_loc = &this->values[0];
-      NumberType *U_loc = left_singluar_vectors ? &(U.values[0]) : NULL;
-      NumberType *VT_loc = right_singluar_vectors ? &(VT.values[0]) : NULL;
+      NumberType *U_loc = left_singluar_vectors ? &(U->values[0]) : nullptr;
+      NumberType *VT_loc = right_singluar_vectors ? &(VT->values[0]) : nullptr;
       int info = 0;
       /*
        * by setting lwork to -1 a workspace query for optimal length of work is performed
@@ -543,8 +587,8 @@ std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix
       work.resize(1);
 
       pgesvd(&jobu,&jobvt,&n_rows,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,
-             & *sv.begin(),U_loc,&U.submatrix_row,&U.submatrix_column,U.descriptor,
-             VT_loc,&VT.submatrix_row,&VT.submatrix_column,VT.descriptor,
+             & *sv.begin(),U_loc,&U->submatrix_row,&U->submatrix_column,U->descriptor,
+             VT_loc,&VT->submatrix_row,&VT->submatrix_column,VT->descriptor,
              &work[0],&lwork,&info);
       AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgesvd", info));
 
@@ -552,8 +596,8 @@ std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix
       work.resize(lwork);
 
       pgesvd(&jobu,&jobvt,&n_rows,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,
-             & *sv.begin(),U_loc,&U.submatrix_row,&U.submatrix_column,U.descriptor,
-             VT_loc,&VT.submatrix_row,&VT.submatrix_column,VT.descriptor,
+             & *sv.begin(),U_loc,&U->submatrix_row,&U->submatrix_column,U->descriptor,
+             VT_loc,&VT->submatrix_row,&VT->submatrix_column,VT->descriptor,
              &work[0],&lwork,&info);
       AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgesvd", info));
     }
@@ -615,6 +659,7 @@ void ScaLAPACKMatrix<NumberType>::least_squares(ScaLAPACKMatrix<NumberType> &B,
             B_loc,&B.submatrix_row,&B.submatrix_column,B.descriptor,&work[0],&lwork,&info);
       AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgels", info));
     }
+  state = LAPACKSupport::State::unusable;
 }
 
 
index 3ede3474ae28ccdba23f0a47e65f73d42d7e5d99..977ff53f5db5a42d8583e0608aac1f660e55b321 100644 (file)
@@ -16,7 +16,7 @@
 #include "../tests.h"
 #include "../lapack/create_matrix.h"
 
-// test eigenpairs_symmetric(const bool, const std::pair<int,int>&, const std::pair<NumberType,NumberType>&)
+// test eigenpairs_symmetric_by_index(const std::pair<unsigned int,unsigned int> &, const bool)
 
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/utilities.h>
@@ -47,7 +47,7 @@ void test(const unsigned int size, const unsigned int block_size, const NumberTy
 
   std::shared_ptr<Utilities::MPI::ProcessGrid> grid = std::make_shared<Utilities::MPI::ProcessGrid>(mpi_communicator,size,size,block_size,block_size);
 
-  pcout << size << " " << block_size << std::endl;
+  pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
 
   // Create SPD matrices of requested size:
   FullMatrix<NumberType> full_A(size);
@@ -90,7 +90,7 @@ void test(const unsigned int size, const unsigned int block_size, const NumberTy
   ScaLAPACKMatrix<NumberType> scalapack_syev (size, grid, block_size);
   scalapack_syev.set_property(LAPACKSupport::Property::symmetric);
   scalapack_syev = full_A;
-  eigenvalues_psyev = scalapack_syev.eigenpairs_symmetric(true);
+  eigenvalues_psyev = scalapack_syev.eigenpairs_symmetric_by_index(std::make_pair(0,size-1),true);
   FullMatrix<NumberType> p_eigenvectors (size,size);
   scalapack_syev.copy_to(p_eigenvectors);
   for (unsigned int i=0; i<max_n_eigenvalues; ++i)
@@ -120,7 +120,7 @@ void test(const unsigned int size, const unsigned int block_size, const NumberTy
   ScaLAPACKMatrix<NumberType> scalapack_syevx_partial (size, grid, block_size);
   scalapack_syevx_partial.set_property(LAPACKSupport::Property::symmetric);
   scalapack_syevx_partial = full_A;
-  eigenvalues_psyevx_partial = scalapack_syevx_partial.eigenpairs_symmetric(true, std::make_pair(size-max_n_eigenvalues+1,size));
+  eigenvalues_psyevx_partial = scalapack_syevx_partial.eigenpairs_symmetric_by_index(std::make_pair(size-max_n_eigenvalues,size-1),true);
   scalapack_syevx_partial.copy_to(p_eigenvectors);
   for (unsigned int i=eigenvalues_psyevx_partial.size()-1; i>0; --i)
     {
index 77e6a56658902b73e1d1aa015f9e4c88de582872..c4d61c891c474279e49332a5d7df3fb461810efd 100644 (file)
@@ -1,4 +1,4 @@
-200 32
+200 32 1 1
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -8,7 +8,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-200 64
+200 64 1 1
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -18,7 +18,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 32
+400 32 1 1
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -28,7 +28,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 64
+400 64 1 1
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -38,7 +38,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 32
+600 32 1 1
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -48,7 +48,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 64
+600 64 1 1
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
index 77e6a56658902b73e1d1aa015f9e4c88de582872..da927e1d1936da824961d5ac08a57865e65b2e1a 100644 (file)
@@ -1,4 +1,4 @@
-200 32
+200 32 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -8,7 +8,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-200 64
+200 64 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -18,7 +18,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 32
+400 32 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -28,7 +28,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 64
+400 64 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -38,7 +38,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 32
+600 32 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -48,7 +48,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 64
+600 64 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
index 77e6a56658902b73e1d1aa015f9e4c88de582872..98129f132b318a8d4cf0698f50e8fe26913fb2cc 100644 (file)
@@ -1,4 +1,4 @@
-200 32
+200 32 2 2
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -8,7 +8,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-200 64
+200 64 2 2
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -18,7 +18,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 32
+400 32 2 2
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -28,7 +28,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 64
+400 64 2 2
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -38,7 +38,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 32
+600 32 2 2
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -48,7 +48,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 64
+600 64 2 2
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
index 77e6a56658902b73e1d1aa015f9e4c88de582872..da927e1d1936da824961d5ac08a57865e65b2e1a 100644 (file)
@@ -1,4 +1,4 @@
-200 32
+200 32 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -8,7 +8,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-200 64
+200 64 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -18,7 +18,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 32
+400 32 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -28,7 +28,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-400 64
+400 64 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -38,7 +38,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 32
+600 32 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
@@ -48,7 +48,7 @@ comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pds
    with respect to the given tolerance also the eigenvectors coincide
 
 
-600 64
+600 64 3 3
 comparing 5 eigenvalues and eigenvectors computed using LAPACK and ScaLAPACK pdsyev:
    with respect to the given tolerance the eigenvalues coincide
    with respect to the given tolerance also the eigenvectors coincide
index 399b235cff7f1359dcf2335ffa156a8176fd14c9..18430dd8d2511437bd66b6297afee77b594a4b8a 100644 (file)
@@ -16,7 +16,7 @@
 #include "../tests.h"
 #include "../lapack/create_matrix.h"
 
-// test compute_SVD(ScaLAPACKMatrix<NumberType>&,ScaLAPACKMatrix<NumberType>&,const bool,const bool)
+// test compute_SVD(ScaLAPACKMatrix<NumberType>*,ScaLAPACKMatrix<NumberType>*)
 
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/utilities.h>
@@ -57,7 +57,7 @@ void test(const unsigned int size, const unsigned int block_size, const NumberTy
   ScaLAPACKMatrix<NumberType> scalapack_A_ev (size, grid_2d, block_size);
   scalapack_A_ev.set_property(LAPACKSupport::Property::symmetric);
   scalapack_A_ev = full_A;
-  std::vector<NumberType> eigenvalues = scalapack_A_ev.eigenpairs_symmetric(true);
+  std::vector<NumberType> eigenvalues = scalapack_A_ev.eigenpairs_symmetric_by_index(std::make_pair(0,size-1),true);
   FullMatrix<NumberType> eigenvectors (size,size);
   scalapack_A_ev.copy_to(eigenvectors);
 
@@ -67,7 +67,7 @@ void test(const unsigned int size, const unsigned int block_size, const NumberTy
   ScaLAPACKMatrix<NumberType> scalapack_VT (size, grid_2d, block_size);
   scalapack_A_sv.set_property(LAPACKSupport::Property::symmetric);
   scalapack_A_sv = full_A;
-  std::vector<NumberType> singular_values = scalapack_A_sv.compute_SVD(scalapack_U,scalapack_VT,true,true);
+  std::vector<NumberType> singular_values = scalapack_A_sv.compute_SVD(&scalapack_U,&scalapack_VT);
   FullMatrix<NumberType> l_singular_vectors (size,size);
   FullMatrix<NumberType> r_singular_vectors (size,size);
   scalapack_U.copy_to(l_singular_vectors);

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.