]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement the solution of complex linear systems using UMFPACK.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 30 Dec 2019 19:07:39 +0000 (12:07 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 6 Jan 2020 19:29:50 +0000 (12:29 -0700)
include/deal.II/lac/sparse_direct.h
source/lac/sparse_direct.cc

index a3546365311580578c112a3b9758728fca685369..d4cd647a068ee52bd4017ea34389dc4a79eb93ea 100644 (file)
@@ -244,6 +244,21 @@ public:
   void
   solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
 
+  /**
+   * Like the previous function, but for a complex-valued right hand side
+   * and solution vector.
+   *
+   * If the matrix that was previously factorized had complex-valued entries,
+   * then the `rhs_and_solution` vector will, upon return from this function,
+   * simply contain the solution of the linear system $Ax=b$. If the matrix
+   * was real-valued, then this is also true, but the solution will simply
+   * be computed by applying the factorized $A^{-1}$ to both the
+   * real and imaginary parts of the right hand side vector.
+   */
+  void
+  solve(Vector<std::complex<double>> &rhs_and_solution,
+        const bool                    transpose = false) const;
+
   /**
    * Same as before, but for block vectors.
    */
@@ -251,6 +266,13 @@ public:
   solve(BlockVector<double> &rhs_and_solution,
         const bool           transpose = false) const;
 
+  /**
+   * Same as before, but for complex-valued block vectors.
+   */
+  void
+  solve(BlockVector<std::complex<double>> &rhs_and_solution,
+        const bool                         transpose = false) const;
+
   /**
    * Call the two functions factorize() and solve() in that order, i.e.
    * perform the whole solution process for the given right hand side vector.
@@ -263,6 +285,15 @@ public:
         Vector<double> &rhs_and_solution,
         const bool      transpose = false);
 
+  /**
+   * Same as before, but for complex-valued solution vectors.
+   */
+  template <class Matrix>
+  void
+  solve(const Matrix &                matrix,
+        Vector<std::complex<double>> &rhs_and_solution,
+        const bool                    transpose = false);
+
   /**
    * Same as before, but for block vectors.
    */
@@ -272,6 +303,15 @@ public:
         BlockVector<double> &rhs_and_solution,
         const bool           transpose = false);
 
+  /**
+   * Same as before, but for complex-valued block vectors.
+   */
+  template <class Matrix>
+  void
+  solve(const Matrix &                     matrix,
+        BlockVector<std::complex<double>> &rhs_and_solution,
+        const bool                         transpose = false);
+
   /**
    * @}
    */
@@ -362,11 +402,20 @@ private:
   sort_arrays(const BlockSparseMatrix<number> &);
 
   /**
-   * The arrays in which we store the data for the solver.
+   * The arrays in which we store the data for the solver. These are documented
+   * in the descriptions of the umfpack_*_symbolic() and umfpack_*_numeric()
+   * functions, but in short:
+   * - `Ap` is the array saying which row starts where in `Ai`
+   * - `Ai` is the array that stores the column indices of nonzero entries
+   * - `Ax` is the array that stores the values of nonzero entries; if the
+   *   matrix is complex-valued, then it stores the real parts
+   * - `Az` is the array that stores the imaginary parts of nonzero entries,
+   *   and is used only if the matrix is complex-valued.
    */
   std::vector<types::suitesparse_index> Ap;
   std::vector<types::suitesparse_index> Ai;
   std::vector<double>                   Ax;
+  std::vector<double>                   Az;
 
   /**
    * Control and work arrays for the solver routines.
index 9a225192ab2a7cb5ee5cc6d0de7a47f0e34e44f7..b158a656164e26c6cf9baa5924815db18b573404 100644 (file)
@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/numbers.h>
 
 #include <deal.II/lac/block_sparse_matrix.h>
 #include <deal.II/lac/sparse_direct.h>
@@ -21,6 +22,7 @@
 #include <deal.II/lac/vector.h>
 
 #include <cerrno>
+#include <complex>
 #include <iostream>
 #include <list>
 #include <typeinfo>
@@ -93,6 +95,11 @@ SparseDirectUMFPACK::clear()
     tmp.swap(Ax);
   }
 
+  {
+    std::vector<double> tmp;
+    tmp.swap(Az);
+  }
+
   umfpack_dl_defaults(control.data());
 }
 
@@ -126,7 +133,11 @@ SparseDirectUMFPACK::sort_arrays(const SparseMatrix<number> &matrix)
       while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
         {
           std::swap(Ai[cursor], Ai[cursor + 1]);
+
           std::swap(Ax[cursor], Ax[cursor + 1]);
+          if (numbers::NumberTraits<number>::is_complex == true)
+            std::swap(Az[cursor], Az[cursor + 1]);
+
           ++cursor;
         }
     }
@@ -145,7 +156,11 @@ SparseDirectUMFPACK::sort_arrays(const SparseMatrixEZ<number> &matrix)
       while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
         {
           std::swap(Ai[cursor], Ai[cursor + 1]);
+
           std::swap(Ax[cursor], Ax[cursor + 1]);
+          if (numbers::NumberTraits<number>::is_complex == true)
+            std::swap(Az[cursor], Az[cursor + 1]);
+
           ++cursor;
         }
     }
@@ -181,7 +196,11 @@ SparseDirectUMFPACK::sort_arrays(const BlockSparseMatrix<number> &matrix)
           while ((element < Ap[row + 1] - 1) && (Ai[element] > Ai[element + 1]))
             {
               std::swap(Ai[element], Ai[element + 1]);
+
               std::swap(Ax[element], Ax[element + 1]);
+              if (numbers::NumberTraits<number>::is_complex == true)
+                std::swap(Az[cursor], Az[cursor + 1]);
+
               ++element;
             }
         }
@@ -194,9 +213,11 @@ template <class Matrix>
 void
 SparseDirectUMFPACK::factorize(const Matrix &matrix)
 {
-  Assert(matrix.m() == matrix.n(), ExcNotQuadratic())
+  Assert(matrix.m() == matrix.n(), ExcNotQuadratic());
 
-    clear();
+  clear();
+
+  using number = typename Matrix::value_type;
 
   n_rows = matrix.m();
   n_cols = matrix.n();
@@ -224,6 +245,8 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix)
   Ap.resize(N + 1);
   Ai.resize(matrix.n_nonzero_elements());
   Ax.resize(matrix.n_nonzero_elements());
+  if (numbers::NumberTraits<number>::is_complex == true)
+    Az.resize(matrix.n_nonzero_elements());
 
   // first fill row lengths array
   Ap[0] = 0;
@@ -250,7 +273,9 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix)
           {
             // write entry into the first free one for this row
             Ai[row_pointers[row]] = p->column();
-            Ax[row_pointers[row]] = p->value();
+            Ax[row_pointers[row]] = std::real(p->value());
+            if (numbers::NumberTraits<number>::is_complex == true)
+              Az[row_pointers[row]] = std::imag(p->value());
 
             // then move pointer ahead
             ++row_pointers[row];
@@ -268,24 +293,45 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix)
   sort_arrays(matrix);
 
   int status;
-  status = umfpack_dl_symbolic(N,
-                               N,
-                               Ap.data(),
-                               Ai.data(),
-                               Ax.data(),
-                               &symbolic_decomposition,
-                               control.data(),
-                               nullptr);
+  if (numbers::NumberTraits<number>::is_complex == false)
+    status = umfpack_dl_symbolic(N,
+                                 N,
+                                 Ap.data(),
+                                 Ai.data(),
+                                 Ax.data(),
+                                 &symbolic_decomposition,
+                                 control.data(),
+                                 nullptr);
+  else
+    status = umfpack_zl_symbolic(N,
+                                 N,
+                                 Ap.data(),
+                                 Ai.data(),
+                                 Ax.data(),
+                                 Az.data(),
+                                 &symbolic_decomposition,
+                                 control.data(),
+                                 nullptr);
   AssertThrow(status == UMFPACK_OK,
               ExcUMFPACKError("umfpack_dl_symbolic", status));
 
-  status = umfpack_dl_numeric(Ap.data(),
-                              Ai.data(),
-                              Ax.data(),
-                              symbolic_decomposition,
-                              &numeric_decomposition,
-                              control.data(),
-                              nullptr);
+  if (numbers::NumberTraits<number>::is_complex == false)
+    status = umfpack_dl_numeric(Ap.data(),
+                                Ai.data(),
+                                Ax.data(),
+                                symbolic_decomposition,
+                                &numeric_decomposition,
+                                control.data(),
+                                nullptr);
+  else
+    status = umfpack_zl_numeric(Ap.data(),
+                                Ai.data(),
+                                Ax.data(),
+                                Az.data(),
+                                symbolic_decomposition,
+                                &numeric_decomposition,
+                                control.data(),
+                                nullptr);
   AssertThrow(status == UMFPACK_OK,
               ExcUMFPACKError("umfpack_dl_numeric", status));
 
@@ -303,6 +349,13 @@ SparseDirectUMFPACK::solve(Vector<double> &rhs_and_solution,
   Assert(Ai.size() != 0, ExcNotInitialized());
   Assert(Ai.size() == Ax.size(), ExcNotInitialized());
 
+  Assert(Az.size() == 0,
+         ExcMessage("You have previously factored a matrix using this class "
+                    "that had complex-valued entries. This then requires "
+                    "applying the factored matrix to a complex-valued "
+                    "vector, but you are only providing a real-valued vector "
+                    "here."));
+
   Vector<double> rhs(rhs_and_solution.size());
   rhs = rhs_and_solution;
 
@@ -326,6 +379,111 @@ SparseDirectUMFPACK::solve(Vector<double> &rhs_and_solution,
 }
 
 
+
+void
+SparseDirectUMFPACK::solve(Vector<std::complex<double>> &rhs_and_solution,
+                           const bool transpose /*=false*/) const
+{
+  // make sure that some kind of factorize() call has happened before
+  Assert(Ap.size() != 0, ExcNotInitialized());
+  Assert(Ai.size() != 0, ExcNotInitialized());
+  Assert(Ai.size() == Ax.size(), ExcNotInitialized());
+
+  // First see whether the matrix that was factorized was complex-valued.
+  // If so, just apply the complex factorization to the vector.
+  if (Az.size() != 0)
+    {
+      Assert(Ax.size() == Az.size(), ExcInternalError());
+
+      // It would be nice if we could just present a pointer to the
+      // first element of the complex-valued solution vector and let
+      // UMFPACK fill both the real and imaginary parts of the solution
+      // vector at that address. UMFPACK calls this the 'packed' format,
+      // and in those cases it only takes one pointer to the entire
+      // vector, rather than a pointer to the real and one pointer to
+      // the imaginary parts of the vector. The problem is that if we
+      // want to pack, then we also need to pack the matrix, and the
+      // functions above have already decided that we don't want to pack
+      // the matrix but instead deal with split format for the matrix,
+      // and then UMFPACK decides that it can't deal with a split matrix
+      // and a packed vector. We have to choose one or the other, not
+      // mix.
+      //
+      // So create four vectors, one each for the real and imaginary parts
+      // of the right hand side and solution.
+      Vector<double> rhs_re(rhs_and_solution.size());
+      Vector<double> rhs_im(rhs_and_solution.size());
+      for (unsigned int i = 0; i < rhs_and_solution.size(); ++i)
+        {
+          rhs_re(i) = std::real(rhs_and_solution(i));
+          rhs_im(i) = std::imag(rhs_and_solution(i));
+        }
+
+      Vector<double> solution_re(rhs_and_solution.size());
+      Vector<double> solution_im(rhs_and_solution.size());
+
+      // Solve the system. note that since UMFPACK wants compressed column
+      // storage instead of the compressed row storage format we use in
+      // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
+      //
+      // Conversely, if we solve for the transpose, we have to use UMFPACK_A
+      // instead.
+      //
+      // Note that for the complex case here, the transpose is selected using
+      // UMFPACK_Aat, not UMFPACK_At.
+      const int status = umfpack_zl_solve(transpose ? UMFPACK_A : UMFPACK_Aat,
+                                          Ap.data(),
+                                          Ai.data(),
+                                          Ax.data(),
+                                          Az.data(),
+                                          solution_re.data(),
+                                          solution_im.data(),
+                                          rhs_re.data(),
+                                          rhs_im.data(),
+                                          numeric_decomposition,
+                                          control.data(),
+                                          nullptr);
+      AssertThrow(status == UMFPACK_OK,
+                  ExcUMFPACKError("umfpack_dl_solve", status));
+
+      // Now put things back together into the output vector
+      for (unsigned int i = 0; i < rhs_and_solution.size(); ++i)
+        rhs_and_solution(i) = {solution_re(i), solution_im(i)};
+    }
+  else
+    {
+      // We have factorized a real-valued matrix, but the rhs and solution
+      // vectors are complex-valued. UMFPACK does not natively support this
+      // case, but we can just apply the factorization to real and imaginary
+      // parts of the right hand side separately
+      const Vector<std::complex<double>> rhs = rhs_and_solution;
+
+      // Get the real part of the right hand side, solve with it, and copy the
+      // results into the result vector by just copying the real output
+      // into the complex-valued result vector (which implies setting the
+      // imaginary parts to zero):
+      Vector<double> rhs_real_or_imag(rhs_and_solution.size());
+      for (unsigned int i = 0; i < rhs.size(); ++i)
+        rhs_real_or_imag(i) = std::real(rhs(i));
+
+      solve(rhs_real_or_imag, transpose);
+
+      rhs_and_solution = rhs_real_or_imag;
+
+      // Then repeat the whole thing with the imaginary part. The copying step
+      // is more complicated now because we can only touch the imaginary
+      // component of the output vector:
+      for (unsigned int i = 0; i < rhs.size(); ++i)
+        rhs_real_or_imag(i) = std::imag(rhs(i));
+
+      solve(rhs_real_or_imag, transpose);
+
+      for (unsigned int i = 0; i < rhs.size(); ++i)
+        rhs_and_solution(i).imag(rhs_real_or_imag(i));
+    }
+}
+
+
 void
 SparseDirectUMFPACK::solve(BlockVector<double> &rhs_and_solution,
                            const bool           transpose /*=false*/) const
@@ -341,6 +499,21 @@ SparseDirectUMFPACK::solve(BlockVector<double> &rhs_and_solution,
 
 
 
+void
+SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &rhs_and_solution,
+                           const bool transpose /*=false*/) const
+{
+  // the UMFPACK functions want a contiguous array of elements, so
+  // there is no way around copying data around. thus, just copy the
+  // data into a regular vector and back
+  Vector<std::complex<double>> tmp(rhs_and_solution.size());
+  tmp = rhs_and_solution;
+  solve(tmp, transpose);
+  rhs_and_solution = tmp;
+}
+
+
+
 template <class Matrix>
 void
 SparseDirectUMFPACK::solve(const Matrix &  matrix,
@@ -352,6 +525,19 @@ SparseDirectUMFPACK::solve(const Matrix &  matrix,
 }
 
 
+
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix &                matrix,
+                           Vector<std::complex<double>> &rhs_and_solution,
+                           const bool                    transpose /*=false*/)
+{
+  factorize(matrix);
+  solve(rhs_and_solution, transpose);
+}
+
+
+
 template <class Matrix>
 void
 SparseDirectUMFPACK::solve(const Matrix &       matrix,
@@ -363,6 +549,18 @@ SparseDirectUMFPACK::solve(const Matrix &       matrix,
 }
 
 
+
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix &                     matrix,
+                           BlockVector<std::complex<double>> &rhs_and_solution,
+                           const bool transpose /*=false*/)
+{
+  factorize(matrix);
+  solve(rhs_and_solution, transpose);
+}
+
+
 #else
 
 
@@ -402,6 +600,17 @@ SparseDirectUMFPACK::solve(Vector<double> &, const bool) const
 
 
 
+void
+SparseDirectUMFPACK::solve(Vector<std::complex<double>> &, const bool) const
+{
+  AssertThrow(
+    false,
+    ExcMessage(
+      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
+}
+
+
+
 void
 SparseDirectUMFPACK::solve(BlockVector<double> &, const bool) const
 {
@@ -412,6 +621,19 @@ SparseDirectUMFPACK::solve(BlockVector<double> &, const bool) const
 }
 
 
+
+void
+SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &,
+                           const bool) const
+{
+  AssertThrow(
+    false,
+    ExcMessage(
+      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
+}
+
+
+
 template <class Matrix>
 void
 SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, const bool)
@@ -424,6 +646,20 @@ SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, const bool)
 
 
 
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix &,
+                           Vector<std::complex<double>> &,
+                           const bool)
+{
+  AssertThrow(
+    false,
+    ExcMessage(
+      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
+}
+
+
+
 template <class Matrix>
 void
 SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, const bool)
@@ -434,6 +670,20 @@ SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, const bool)
       "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
 }
 
+
+
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix &,
+                           BlockVector<std::complex<double>> &,
+                           const bool)
+{
+  AssertThrow(
+    false,
+    ExcMessage(
+      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
+}
+
 #endif
 
 
@@ -497,17 +747,23 @@ SparseDirectUMFPACK::n() const
 
 
 // explicit instantiations for SparseMatrixUMFPACK
-#define InstantiateUMFPACK(MatrixType)                              \
-  template void SparseDirectUMFPACK::factorize(const MatrixType &); \
-  template void SparseDirectUMFPACK::solve(const MatrixType &,      \
-                                           Vector<double> &,        \
-                                           bool);                   \
-  template void SparseDirectUMFPACK::solve(const MatrixType &,      \
-                                           BlockVector<double> &,   \
-                                           bool);                   \
-  template void SparseDirectUMFPACK::initialize(const MatrixType &, \
+#define InstantiateUMFPACK(MatrixType)                                     \
+  template void SparseDirectUMFPACK::factorize(const MatrixType &);        \
+  template void SparseDirectUMFPACK::solve(const MatrixType &,             \
+                                           Vector<double> &,               \
+                                           const bool);                    \
+  template void SparseDirectUMFPACK::solve(const MatrixType &,             \
+                                           Vector<std::complex<double>> &, \
+                                           const bool);                    \
+  template void SparseDirectUMFPACK::solve(const MatrixType &,             \
+                                           BlockVector<double> &,          \
+                                           const bool);                    \
+  template void SparseDirectUMFPACK::solve(                                \
+    const MatrixType &, BlockVector<std::complex<double>> &, const bool);  \
+  template void SparseDirectUMFPACK::initialize(const MatrixType &,        \
                                                 const AdditionalData)
 
+// Instantiate everything for real-valued matrices
 InstantiateUMFPACK(SparseMatrix<double>);
 InstantiateUMFPACK(SparseMatrix<float>);
 InstantiateUMFPACK(SparseMatrixEZ<double>);
@@ -515,5 +771,11 @@ InstantiateUMFPACK(SparseMatrixEZ<float>);
 InstantiateUMFPACK(BlockSparseMatrix<double>);
 InstantiateUMFPACK(BlockSparseMatrix<float>);
 
+// Now also for complex-valued matrices
+InstantiateUMFPACK(SparseMatrix<std::complex<double>>);
+InstantiateUMFPACK(SparseMatrix<std::complex<float>>);
+InstantiateUMFPACK(BlockSparseMatrix<std::complex<double>>);
+InstantiateUMFPACK(BlockSparseMatrix<std::complex<float>>);
+
 
 DEAL_II_NAMESPACE_CLOSE

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.