]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Return eigenvectors as full matrix. Fix complex case
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 13 Mar 2023 17:42:03 +0000 (18:42 +0100)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Tue, 14 Mar 2023 13:17:16 +0000 (14:17 +0100)
include/deal.II/lac/lapack_full_matrix.h
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_38.cc
tests/lapack/full_matrix_38.output

index 6d915861777ae5a18ebae91a332422d8917bf450..78b1ad807ef7fdc3a6e12b6856c2da46a0ab3e89 100644 (file)
@@ -75,14 +75,12 @@ public:
    */
   explicit LAPACKFullMatrix(const size_type size = 0);
 
-
   /**
    * Constructor. Initialize the matrix as a rectangular matrix $\rm{rows}
    * \times \rm{cols}$.
    */
   LAPACKFullMatrix(const size_type rows, const size_type cols);
 
-
   /**
    * Copy constructor. This constructor does a deep copy of the matrix.
    * Therefore, it poses a possible efficiency problem, if for example,
@@ -481,8 +479,9 @@ public:
    * @note It is assumed that @p A, @p B and @p V have compatible sizes and that
    * @p C already has the right size.
    *
-   * @note This function is not provided by LAPACK. The function first forms $\rm{diag}(\mathbf V) \cdot \mathbf B$ product and
-   * then uses Xgemm function.
+   * @note This function is not provided by LAPACK. The function first forms
+   * $\rm{diag}(\mathbf V) \cdot \mathbf B$ product and then uses the Xgemm
+   * function.
    */
   void
   Tmmult(LAPACKFullMatrix<number> &      C,
@@ -830,30 +829,32 @@ public:
   compute_inverse_svd_with_kernel(const unsigned int kernel_size);
 
   /**
-   * Retrieve eigenvalue after compute_eigenvalues() was called.
+   * Retrieve eigenvalue after compute_eigenvalues() was called. The return
+   * type is always a complex number: In case the underlying matrix is
+   * real-valued, the type `std::complex<number>` is returned, whereas
+   * `number` is returned for complex numbers.
    */
-  std::complex<number>
+  std::complex<typename numbers::NumberTraits<number>::real_type>
   eigenvalue(const size_type i) const;
 
   /**
-   * After a call to compute_eigenvalues(), this function returns all (right)
-   * eigenvectors as returned by LAPACK. This means that eigenvectors are
-   * contained in column-major ordering in a matrix associated to the given
-   * flat vector. Note that for real-valued matrices, there might appear
-   * complex eigenvalues with complex-conjugate values and eigenvectors. For
-   * those cases, LAPACK places a column of n() entries with the real part and
-   * the next column (corresponding to the complex conjugate eigenvalue) for
-   * the imaginary part of the eigenvector.
+   * After a call to compute_eigenvalues(), this function returns the $n\times
+   * n$ matrix of (right) eigenvectors in a decomposition of the form $A V = V
+   * \Lambda$. Note that this function constructs the associated matrix on the
+   * fly, since LAPACK packs complex-conjugate eigenvalue/eigenvector pairs of
+   * real-valued matrices into a real-valued return matrix. This call only
+   * succeeds in case the respective flag `right_eigenvectors` in
+   * compute_eigenvalues() has been set to true.
    */
-  std::vector<number>
+  FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
   get_right_eigenvectors() const;
 
   /**
    * Return the matrix of left eigenvectors after a call to
-   * compute_eigenvalues(), following the same convention as
+   * compute_eigenvalues(), following the same principles as
    * get_right_eigenvectors().
    */
-  std::vector<number>
+  FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
   get_left_eigenvectors() const;
 
   /**
@@ -1037,6 +1038,7 @@ LAPACKFullMatrix<number>::set(const size_type i,
 }
 
 
+
 template <typename number>
 inline typename LAPACKFullMatrix<number>::size_type
 LAPACKFullMatrix<number>::m() const
@@ -1044,6 +1046,8 @@ LAPACKFullMatrix<number>::m() const
   return static_cast<size_type>(this->n_rows());
 }
 
+
+
 template <typename number>
 inline typename LAPACKFullMatrix<number>::size_type
 LAPACKFullMatrix<number>::n() const
@@ -1051,6 +1055,8 @@ LAPACKFullMatrix<number>::n() const
   return static_cast<size_type>(this->n_cols());
 }
 
+
+
 template <typename number>
 template <typename MatrixType>
 inline void
@@ -1109,6 +1115,7 @@ LAPACKFullMatrix<number>::fill(const MatrixType &M,
 }
 
 
+
 template <typename number>
 template <typename number2>
 void
@@ -1122,6 +1129,7 @@ LAPACKFullMatrix<number>::vmult(Vector<number2> &,
 }
 
 
+
 template <typename number>
 template <typename number2>
 void
@@ -1134,6 +1142,7 @@ LAPACKFullMatrix<number>::vmult_add(Vector<number2> &,
 }
 
 
+
 template <typename number>
 template <typename number2>
 void
@@ -1147,6 +1156,7 @@ LAPACKFullMatrix<number>::Tvmult(Vector<number2> &,
 }
 
 
+
 template <typename number>
 template <typename number2>
 void
@@ -1154,54 +1164,50 @@ LAPACKFullMatrix<number>::Tvmult_add(Vector<number2> &,
                                      const Vector<number2> &) const
 {
   Assert(false,
-         ExcMessage(
-           "LAPACKFullMatrix<number>::Tvmult_add must be called with a "
-           "matching Vector<double> vector type."));
+         ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called "
+                    "with a matching Vector<double> vector type."));
 }
 
 
-template <typename number>
-inline std::complex<number>
-LAPACKFullMatrix<number>::eigenvalue(const size_type i) const
-{
-  Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
-  Assert(wr.size() == this->n_rows(), ExcInternalError());
-  Assert(wi.size() == this->n_rows(), ExcInternalError());
-  AssertIndexRange(i, this->n_rows());
-
-  if (numbers::NumberTraits<number>::is_complex)
-    return std::complex<number>(wi[i]);
-  else
-    return std::complex<number>(wr[i], wi[i]);
-}
-
 
-template <typename number>
-inline std::vector<number>
-LAPACKFullMatrix<number>::get_right_eigenvectors() const
+namespace internal
 {
-  Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
-  Assert(vr.size() == this->n_rows() * this->n_cols(),
-         ExcMessage("Right eigenvectors are not available! Did you "
-                    "set the associated flag in compute_eigenvalues()"));
+  namespace LAPACKFullMatrixImplementation
+  {
+    template <typename RealNumber>
+    std::complex<RealNumber>
+    pack_complex(const RealNumber &real_part, const RealNumber &imaginary_part)
+    {
+      return std::complex<RealNumber>(real_part, imaginary_part);
+    }
+
+    // The eigenvalues in LAPACKFullMatrix with complex-valued matrices are
+    // contained in the 'wi' array, ignoring the 'wr' array.
+    template <typename Number>
+    std::complex<Number>
+    pack_complex(const Number &, const std::complex<Number> &complex_number)
+    {
+      return complex_number;
+    }
+  } // namespace LAPACKFullMatrixImplementation
+} // namespace internal
 
-  return vr;
-}
 
 
 template <typename number>
-inline std::vector<number>
-LAPACKFullMatrix<number>::get_left_eigenvectors() const
+inline std::complex<typename numbers::NumberTraits<number>::real_type>
+LAPACKFullMatrix<number>::eigenvalue(const size_type i) const
 {
   Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
-  Assert(vl.size() == this->n_rows() * this->n_cols(),
-         ExcMessage("Left eigenvectors are not available! Did you "
-                    "set the associated flag in compute_eigenvalues()"));
+  Assert(wr.size() == this->n_rows(), ExcInternalError());
+  Assert(wi.size() == this->n_rows(), ExcInternalError());
+  AssertIndexRange(i, this->n_rows());
 
-  return vl;
+  return internal::LAPACKFullMatrixImplementation::pack_complex(wr[i], wi[i]);
 }
 
 
+
 template <typename number>
 inline number
 LAPACKFullMatrix<number>::singular_value(const size_type i) const
@@ -1214,6 +1220,7 @@ LAPACKFullMatrix<number>::singular_value(const size_type i) const
 }
 
 
+
 template <typename number>
 inline const LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::get_svd_u() const
@@ -1225,6 +1232,7 @@ LAPACKFullMatrix<number>::get_svd_u() const
 }
 
 
+
 template <typename number>
 inline const LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::get_svd_vt() const
index 765aea45297ea83a76e46bbaf9a218a4fe6bbded..1d61dea4cbc0542efd19bca6298736767c947d06 100644 (file)
@@ -95,6 +95,8 @@ namespace internal
            &info);
     }
 
+
+
     template <typename T>
     void
     geev_helper(const char                      vl,
@@ -125,10 +127,10 @@ namespace internal
         Assert(static_cast<std::size_t>(n_rows * n_rows) <=
                  right_eigenvectors.size(),
                ExcInternalError());
-      Assert(std::max<std::size_t>(1, work_flag) <= real_work.size(),
-             ExcInternalError());
       Assert(work_flag == -1 ||
-               std::max<long int>(1, 2 * n_rows) <= (work_flag),
+               std::max<std::size_t>(1, work_flag) <= real_work.size(),
+             ExcInternalError());
+      Assert(work_flag == -1 || std::max<long int>(1, 2 * n_rows) <= work_flag,
              ExcInternalError());
 
       geev(&vl,
@@ -240,6 +242,8 @@ namespace internal
   } // namespace LAPACKFullMatrixImplementation
 } // namespace internal
 
+
+
 template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type n)
   : TransposeTable<number>(n, n)
@@ -248,6 +252,7 @@ LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type n)
 {}
 
 
+
 template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type m, const size_type n)
   : TransposeTable<number>(m, n)
@@ -256,6 +261,7 @@ LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type m, const size_type n)
 {}
 
 
+
 template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
   : TransposeTable<number>(M)
@@ -264,6 +270,7 @@ LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
 {}
 
 
+
 template <typename number>
 LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::operator=(const LAPACKFullMatrix<number> &M)
@@ -370,6 +377,7 @@ LAPACKFullMatrix<number>::reinit(const size_type m, const size_type n)
 }
 
 
+
 template <typename number>
 template <typename number2>
 LAPACKFullMatrix<number> &
@@ -387,6 +395,7 @@ LAPACKFullMatrix<number>::operator=(const FullMatrix<number2> &M)
 }
 
 
+
 template <typename number>
 template <typename number2>
 LAPACKFullMatrix<number> &
@@ -404,6 +413,7 @@ LAPACKFullMatrix<number>::operator=(const SparseMatrix<number2> &M)
 }
 
 
+
 template <typename number>
 LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::operator=(const number d)
@@ -419,6 +429,7 @@ LAPACKFullMatrix<number>::operator=(const number d)
 }
 
 
+
 template <typename number>
 LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::operator*=(const number factor)
@@ -447,6 +458,7 @@ LAPACKFullMatrix<number>::operator*=(const number factor)
 }
 
 
+
 template <typename number>
 LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::operator/=(const number factor)
@@ -780,6 +792,7 @@ LAPACKFullMatrix<number>::vmult(Vector<number> &      w,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::Tvmult(Vector<number> &      w,
@@ -918,6 +931,7 @@ LAPACKFullMatrix<number>::Tvmult(Vector<number> &      w,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::vmult_add(Vector<number> &      w,
@@ -927,6 +941,7 @@ LAPACKFullMatrix<number>::vmult_add(Vector<number> &      w,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::Tvmult_add(Vector<number> &      w,
@@ -936,6 +951,7 @@ LAPACKFullMatrix<number>::Tvmult_add(Vector<number> &      w,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::mmult(LAPACKFullMatrix<number> &      C,
@@ -970,6 +986,7 @@ LAPACKFullMatrix<number>::mmult(LAPACKFullMatrix<number> &      C,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::mmult(FullMatrix<number> &            C,
@@ -1088,6 +1105,7 @@ LAPACKFullMatrix<number>::transpose(LAPACKFullMatrix<number> &B) const
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::scale_rows(const Vector<number> &V)
@@ -1162,6 +1180,7 @@ LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &      C,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::Tmmult(FullMatrix<number> &            C,
@@ -1197,6 +1216,7 @@ LAPACKFullMatrix<number>::Tmmult(FullMatrix<number> &            C,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::mTmult(LAPACKFullMatrix<number> &      C,
@@ -1290,6 +1310,7 @@ LAPACKFullMatrix<number>::mTmult(FullMatrix<number> &            C,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::TmTmult(LAPACKFullMatrix<number> &      C,
@@ -1324,6 +1345,7 @@ LAPACKFullMatrix<number>::TmTmult(LAPACKFullMatrix<number> &      C,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::TmTmult(FullMatrix<number> &            C,
@@ -1359,6 +1381,7 @@ LAPACKFullMatrix<number>::TmTmult(FullMatrix<number> &            C,
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_lu_factorization()
@@ -1647,6 +1670,7 @@ LAPACKFullMatrix<number>::compute_svd()
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
@@ -1876,6 +1900,7 @@ LAPACKFullMatrix<number>::determinant() const
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_eigenvalues(const bool right, const bool left)
@@ -1932,6 +1957,7 @@ LAPACKFullMatrix<number>::compute_eigenvalues(const bool right, const bool left)
 
   // resize workspace array
   work.resize(lwork);
+  real_work.resize(lwork);
 
   // Finally compute the eigenvalues.
   internal::LAPACKFullMatrixImplementation::geev_helper(jobvl,
@@ -1956,6 +1982,111 @@ LAPACKFullMatrix<number>::compute_eigenvalues(const bool right, const bool left)
 }
 
 
+
+namespace
+{
+  // This function extracts complex eigenvectors from the underlying 'number'
+  // array 'vr' of the LAPACK eigenvalue routine. For real-valued matrices
+  // addressed by this function specialization, we might get complex
+  // eigenvalues, which come in complex-conjugate pairs. In LAPACK, a compact
+  // storage scheme is applied that stores the real and imaginary part of
+  // eigenvectors only once, putting the real parts in one column and the
+  // imaginary part in the next of a real-valued array. Here, we do the
+  // unpacking into the usual complex values.
+  template <typename RealNumber>
+  void
+  unpack_lapack_eigenvector_and_increment_index(
+    const std::vector<RealNumber> &       vr,
+    const std::complex<RealNumber> &      eigenvalue,
+    FullMatrix<std::complex<RealNumber>> &result,
+    unsigned int &                        index)
+  {
+    const std::size_t n = result.n();
+    if (eigenvalue.imag() != 0.)
+      {
+        for (std::size_t j = 0; j < n; ++j)
+          {
+            result(j, index).real(vr[index * n + j]);
+            result(j, index + 1).real(vr[index * n + j]);
+            result(j, index).imag(vr[(index + 1) * n + j]);
+            result(j, index + 1).imag(-vr[(index + 1) * n + j]);
+          }
+
+        // we filled two columns with the complex-conjugate pair, so increment
+        // returned index by 2
+        index += 2;
+      }
+    else
+      {
+        for (unsigned int j = 0; j < n; ++j)
+          result(j, index).real(vr[index * n + j]);
+
+        // real-valued case, we only filled one column
+        ++index;
+      }
+  }
+
+  // This specialization fills the eigenvectors for complex-valued matrices,
+  // in which case we simply read off the entry in the 'vr' array.
+  template <typename ComplexNumber>
+  void
+  unpack_lapack_eigenvector_and_increment_index(
+    const std::vector<ComplexNumber> &vr,
+    const ComplexNumber &,
+    FullMatrix<ComplexNumber> &result,
+    unsigned int &             index)
+  {
+    const std::size_t n = result.n();
+    for (unsigned int j = 0; j < n; ++j)
+      result(j, index) = vr[index * n + j];
+
+    // complex-valued case always only fills a single column
+    ++index;
+  }
+} // namespace
+
+
+
+template <typename number>
+FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+LAPACKFullMatrix<number>::get_right_eigenvectors() const
+{
+  Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
+  Assert(vr.size() == this->n_rows() * this->n_cols(),
+         ExcMessage("Right eigenvectors are not available! Did you "
+                    "set the associated flag in compute_eigenvalues()?"));
+
+  FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+    result(n(), n());
+
+  for (unsigned int i = 0; i < n();)
+    unpack_lapack_eigenvector_and_increment_index(vr, eigenvalue(i), result, i);
+
+  return result;
+}
+
+
+
+template <typename number>
+FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+LAPACKFullMatrix<number>::get_left_eigenvectors() const
+{
+  Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
+  Assert(vl.size() == this->n_rows() * this->n_cols(),
+         ExcMessage("Left eigenvectors are not available! Did you "
+                    "set the associated flag in compute_eigenvalues()?"));
+
+  FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
+    result(n(), n());
+
+  for (unsigned int i = 0; i < n();)
+    unpack_lapack_eigenvector_and_increment_index(vl, eigenvalue(i), result, i);
+
+  return result;
+}
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(
@@ -2067,6 +2198,7 @@ LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
@@ -2192,6 +2324,7 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
@@ -2284,6 +2417,7 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::print_formatted(std::ostream &     out,
index 2475522972235dfeb66a9faa44270cc2a889ac5b..f84e30fb55a58dd4f8023befb9d525d7187c22d1 100644 (file)
@@ -57,13 +57,10 @@ check_matrix(const double *matrix_pointer)
       deallog << lambda << " ";
     }
   deallog << std::endl;
-  deallog << "Right eigenvectors ";
-  for (const auto d : LA.get_right_eigenvectors())
-    deallog << d << " ";
-  deallog << std::endl;
-  deallog << "Left eigenvectors ";
-  for (const auto d : LA.get_left_eigenvectors())
-    deallog << d << " ";
+  deallog << "Right eigenvectors" << std::endl;
+  LA.get_right_eigenvectors().print_formatted(deallog.get_file_stream());
+  deallog << "Left eigenvectors" << std::endl;
+  LA.get_left_eigenvectors().print_formatted(deallog.get_file_stream());
   deallog << std::endl;
 }
 
index 1dd51a200a44f25c73f34b015057f77fc5388483..28e5a4ca21713ba07fb79e970d86702b14fe27ca 100644 (file)
@@ -3,11 +3,21 @@ DEAL::Checking matrix
 4.000e+00  -1.000e+00 
 -1.000e+00 4.000e+00  
 DEAL::Eigenvalues: (5.00000,0.00000) (3.00000,0.00000) 
-DEAL::Right eigenvectors 0.707107 -0.707107 0.707107 0.707107 
-DEAL::Left eigenvectors 0.707107 -0.707107 0.707107 0.707107 
+DEAL::Right eigenvectors
+(7.071e-01,0.000e+00) (7.071e-01,0.000e+00) 
+(-7.071e-01,0.000e+00) (7.071e-01,0.000e+00) 
+DEAL::Left eigenvectors
+(7.071e-01,0.000e+00) (7.071e-01,0.000e+00) 
+(-7.071e-01,0.000e+00) (7.071e-01,0.000e+00) 
+DEAL::
 DEAL::Checking matrix 
 1.000e+00  -1.000e+00 
 1.000e+00  1.000e+00  
 DEAL::Eigenvalues: (1.00000,1.00000) (1.00000,-1.00000) 
-DEAL::Right eigenvectors 0.707107 0.00000 0.00000 -0.707107 
-DEAL::Left eigenvectors -0.707107 0.00000 0.00000 0.707107 
+DEAL::Right eigenvectors
+(7.071e-01,0.000e+00) (7.071e-01,0.000e+00) 
+(0.000e+00,-7.071e-01) (0.000e+00,7.071e-01) 
+DEAL::Left eigenvectors
+(-7.071e-01,0.000e+00) (-7.071e-01,0.000e+00) 
+(0.000e+00,7.071e-01) (0.000e+00,-7.071e-01) 
+DEAL::

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.