]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add additional data and tests final commit 18466/head
authorDavide Polverino <davide.polverino@phd.unipi.it>
Tue, 20 May 2025 13:30:26 +0000 (15:30 +0200)
committerDavide Polverino <davide.polverino@phd.unipi.it>
Thu, 22 May 2025 08:28:57 +0000 (10:28 +0200)
Add additional data and tests

Added comments, sym/posdef matrix, refactored test and basic interface ready for pull request

Minor fix on comments

Removed symmetric feature for the time being

Fixed username and some comments

Add changelog

Added again the symmetric functionality

Apply review

Follow up fixes

Additional fixes

13 files changed:
doc/news/changes/minor/20250520Polverino [new file with mode: 0644]
include/deal.II/lac/sparse_direct.h
source/lac/sparse_direct.cc
tests/mumps/mumps_01.cc
tests/mumps/mumps_01.output
tests/mumps/mumps_02.cc
tests/mumps/mumps_02.output
tests/mumps/mumps_03.cc [new file with mode: 0644]
tests/mumps/mumps_03.output [new file with mode: 0644]
tests/mumps/mumps_04.cc [new file with mode: 0644]
tests/mumps/mumps_04.output [new file with mode: 0644]
tests/mumps/mumps_05.cc [new file with mode: 0644]
tests/mumps/mumps_05.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20250520Polverino b/doc/news/changes/minor/20250520Polverino
new file mode 100644 (file)
index 0000000..f21ef76
--- /dev/null
@@ -0,0 +1,4 @@
+Improved: The SparseDirectMUMPS class now takes an AdditionalData to
+control the MUMPS execution.
+<br>
+(Davide Polverino, 2025/05/20)
index c3f0c7df7658de3426940a494895a1ce11e7be62..af369422c8406e484b8fa32a783d90e80d45fa84 100644 (file)
@@ -452,97 +452,202 @@ private:
  */
 class SparseDirectMUMPS : public EnableObserverPointer
 {
-private:
-#ifdef DEAL_II_WITH_MUMPS
-  DMUMPS_STRUC_C id;
-#endif // DEAL_II_WITH_MUMPS
+public:
+  /**
+   * Declare type for container size.
+   */
+  using size_type = types::global_dof_index;
 
-  double                 *a;
-  std::vector<double>     rhs;
-  int                    *irn;
-  int                    *jcn;
-  types::global_dof_index n;
-  types::global_dof_index nz;
+  /**
+   * The AdditionalData struct contains data for controlling the MUMPS
+   * execution.
+   */
+  struct AdditionalData
+  {
+    /**
+     * Struct that contains the data for the Block Low-Rank approximation.
+     */
+    struct BlockLowRank
+    {
+      BlockLowRank(const bool   blr_ucfs          = false,
+                   const double lowrank_threshold = 1e-8)
+        : blr_ucfs(blr_ucfs)
+        , lowrank_threshold(lowrank_threshold)
+      {}
+
+      /**
+       * If true, the Block Low-Rank approximation is used with the UCFS
+       * algorithm, an algorithm with higher compression than the standard one.
+       */
+      bool blr_ucfs;
+
+      /**
+       * Threshold for the low-rank truncation of the blocks.
+       */
+      double lowrank_threshold;
+    };
+
+    AdditionalData(const bool          output_details    = false,
+                   const bool          error_statistics  = false,
+                   const bool          symmetric         = false,
+                   const bool          posdef            = false,
+                   const bool          blr_factorization = false,
+                   const BlockLowRank &blr               = BlockLowRank())
+      : output_details(output_details)
+      , error_statistics(error_statistics)
+      , symmetric(symmetric)
+      , posdef(posdef)
+      , blr_factorization(blr_factorization)
+      , blr(blr)
+    {}
+
+    /*
+     * If true, the MUMPS solver will print out details of the execution.
+     */
+    bool output_details;
+
+    /*
+     * If true, the MUMPS solver will print out error statistics.
+     */
+    bool error_statistics;
+
+    /*
+     * If true, the MUMPS solver will use the symmetric factorization. This is
+     * only possible if the matrix is symmetric.
+     */
+    bool symmetric;
+
+    /*
+     * If true, the MUMPS solver will use the positive definite factorization.
+     * This is only possible if the matrix is symmetric and positive definite.
+     */
+    bool posdef;
+
+    /*
+     * If true, the MUMPS solver will use the Block Low-Rank factorization.
+     */
+    bool blr_factorization;
+
+    /*
+     * Stores Block Low-Rank approximation settings to be used in MUMPS
+     * factorization, if enabled.
+     */
+    BlockLowRank blr;
+  };
+
+  /**
+   * Constructor, takes <tt>AdditionalData</tt> to control MUMPS execution.
+   */
+  SparseDirectMUMPS(const AdditionalData &additional_data = AdditionalData());
+
+  /**
+   * Destructor.
+   */
+  ~SparseDirectMUMPS();
 
   /**
-   * This function initializes a MUMPS instance and hands over the system's
-   * matrix <tt>matrix</tt>.
+   * Exception.
+   */
+  DeclException0(ExcInitializeAlreadyCalled);
+
+  /**
+   * This function computes the LU factorization
+   * of the system's matrix <tt>matrix</tt> with the options given in the
+   * constructor's <tt>AdditionalData</tt>.
    */
   template <class Matrix>
   void
-  initialize_matrix(const Matrix &matrix);
+  initialize(const Matrix &matrix);
 
   /**
-   * Copy the computed solution into the solution vector.
+   * A function in which the inverse of the matrix is applied to the input
+   * vector <tt>src</tt> and the solution is written into the output vector
+   * <tt>dst</tt>.
    */
   void
-  copy_solution(Vector<double> &vector);
+  vmult(Vector<double> &dst, const Vector<double> &src) const;
 
   /**
-   *
+   * A function in which the inverse of the transposed matrix is applied to the
+   * input vector <tt>src</tt> and the solution is written into the output
+   * vector <tt>dst</tt>.
    */
   void
-  copy_rhs_to_mumps(const Vector<double> &rhs);
+  Tvmult(Vector<double> &dst, const Vector<double> &src) const;
 
   /**
-   * Flags storing whether the function <tt>initialize ()</tt> has already
-   * been called.
+   * A function that returns the ICNTL integer array from MUMPS.
+   * The ICNTL array contains integer control parameters for the MUMPS solver.
+   * Keep in mind that MUMPS is a fortran library and the documentation refers
+   * to indices into this array starting from one rather than from zero. To
+   * select the correct index one can use a macro like this:
+   * `#define ICNTL(I) icntl[(I)-1]`. In the MUMPS documentation there is the
+   * description of each parameter of the array. Be aware that ownership of the
+   * array remains in the current class rather than with the caller of this
+   * function.
    */
-  bool initialize_called;
+  int *
+  get_icntl();
+
+private:
+#ifdef DEAL_II_WITH_MUMPS
+  mutable DMUMPS_STRUC_C id;
+#endif // DEAL_II_WITH_MUMPS
 
-public:
   /**
-   * Declare type for container size.
+   * a contains the actual values of the matrix. a[k] is the value of the matrix
+   * entry (i,j) if irn[k] == i and jcn[k] == j.
    */
-  using size_type = types::global_dof_index;
+  double *a;
 
   /**
-   * Constructor
+   * The right-hand side vector. This is the right hand side of the system.
    */
-  SparseDirectMUMPS();
+  mutable std::vector<double> rhs;
 
   /**
-   * Destructor
+   * irn contains the row indices of the non-zero entries of the matrix.
    */
-  ~SparseDirectMUMPS();
+  int *irn;
 
   /**
-   * Exception
+   * jcn contains the column indices of the non-zero entries of the matrix.
    */
-  DeclException0(ExcInitializeAlreadyCalled);
+  int *jcn;
 
   /**
-   * This function initializes a MUMPS instance and hands over the system's
-   * matrix <tt>matrix</tt> and right-hand side <tt>rhs_vector</tt> to the
-   * solver.
+   * The number of rows of the matrix. The matrix is square.
    */
-  template <class Matrix>
-  void
-  initialize(const Matrix &matrix, const Vector<double> &rhs_vector);
+  types::global_dof_index n;
+
+  /**
+   * The number of non-zero entries in the matrix.
+   */
+  std::size_t nnz;
 
   /**
-   * This function initializes a MUMPS instance and computes the factorization
-   * of the system's matrix <tt>matrix</tt>.
+   * This function hands over to MUMPS the system's <tt>matrix</tt>.
    */
   template <class Matrix>
   void
-  initialize(const Matrix &matrix);
+  initialize_matrix(const Matrix &matrix);
 
   /**
-   * A function in which the linear system is solved and the solution vector
-   * is copied into the given <tt>vector</tt>. The right-hand side need to be
-   * supplied in initialize(matrix, vector);
+   * Copy the computed solution into the solution vector.
    */
   void
-  solve(Vector<double> &vector);
+  copy_solution(Vector<double> &vector) const;
 
   /**
-   * A function in which the inverse of the matrix is applied to the input
-   * vector <tt>src</tt> and the solution is written into the output vector
-   * <tt>dst</tt>.
+   * Copy the right-hand side vector into the MUMPS instance.
    */
   void
-  vmult(Vector<double> &dst, const Vector<double> &src);
+  copy_rhs_to_mumps(const Vector<double> &rhs) const;
+
+  /**
+   * Struct that holds the additional data for the MUMPS solver.
+   */
+  AdditionalData additional_data;
 };
 
 DEAL_II_NAMESPACE_CLOSE
index f59d6a93166abe0d1bc9b3dd14beb7a097363368..48f9d4763356c72f2ca9851a04867ec60152a274 100644 (file)
@@ -861,9 +861,56 @@ SparseDirectUMFPACK::n() const
 
 #ifdef DEAL_II_WITH_MUMPS
 
-SparseDirectMUMPS::SparseDirectMUMPS()
-  : initialize_called(false)
-{}
+SparseDirectMUMPS::SparseDirectMUMPS(const AdditionalData &data)
+{
+  this->additional_data = data;
+  // Initialize MUMPS instance:
+  id.job = -1;
+  id.par = 1;
+
+  Assert(!(additional_data.symmetric == false &&
+           additional_data.posdef == true),
+         ExcMessage(
+           "You can't have a positive definite matrix that is not symmetric."));
+
+  if (additional_data.symmetric == true)
+    {
+      if (additional_data.posdef == true)
+        id.sym = 1;
+      else
+        id.sym = 2;
+    }
+  else
+    id.sym = 0;
+
+  // Use MPI_COMM_WORLD as communicator
+  id.comm_fortran = -987654;
+  dmumps_c(&id);
+
+  if (additional_data.output_details == false)
+    {
+      // No outputs
+      id.icntl[0] = -1;
+      id.icntl[1] = -1;
+      id.icntl[2] = -1;
+      id.icntl[3] = 0;
+    }
+
+  if (additional_data.error_statistics == true)
+    id.icntl[10] = 2;
+
+  if (additional_data.blr_factorization == true)
+    {
+      id.icntl[34] = 2;
+      if (additional_data.blr.blr_ucfs == true)
+        id.icntl[35] = 1;
+      Assert(additional_data.blr.lowrank_threshold > 0,
+             ExcMessage("Lowrank threshold must be positive."));
+      id.cntl[6] = additional_data.blr.lowrank_threshold;
+    }
+}
+
+
 
 SparseDirectMUMPS::~SparseDirectMUMPS()
 {
@@ -879,93 +926,81 @@ SparseDirectMUMPS::~SparseDirectMUMPS()
     }
 }
 
+
+
 template <class Matrix>
 void
 SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
 {
   Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square."));
 
-
-  // Check we haven't been here before:
-  Assert(initialize_called == false, ExcInitializeAlreadyCalled());
-
-  // Initialize MUMPS instance:
-  id.job = -1;
-  id.par = 1;
-  id.sym = 0;
-
-  // Use MPI_COMM_WORLD as communicator
-  id.comm_fortran = -987654;
-  dmumps_c(&id);
-
-  // Hand over matrix and right-hand side
+  // Hand over matrix to MUMPS as centralized assembled matrix
   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
     {
-      // Objects denoting a MUMPS data structure:
-      //
       // Set number of unknowns
       n = matrix.n();
 
       // number of nonzero elements in matrix
-      nz = matrix.n_actually_nonzero_elements();
+      nnz = matrix.n_actually_nonzero_elements();
 
       // representation of the matrix
-      a = new double[nz];
+      a = new double[nnz];
 
       // matrix indices pointing to the row and column dimensions
       // respectively of the matrix representation above (a): ie. a[k] is
       // the matrix element (irn[k], jcn[k])
-      irn = new int[nz];
-      jcn = new int[nz];
+      irn = new int[nnz];
+      jcn = new int[nnz];
 
-      size_type index = 0;
+      size_type n_non_zero_elements = 0;
 
       // loop over the elements of the matrix row by row, as suggested in
       // the documentation of the sparse matrix iterator class
-      for (size_type row = 0; row < matrix.m(); ++row)
+      if (additional_data.symmetric == true)
         {
-          for (typename Matrix::const_iterator ptr = matrix.begin(row);
-               ptr != matrix.end(row);
-               ++ptr)
-            if (std::abs(ptr->value()) > 0.0)
-              {
-                a[index]   = ptr->value();
-                irn[index] = row + 1;
-                jcn[index] = ptr->column() + 1;
-                ++index;
-              }
+          for (size_type row = 0; row < matrix.m(); ++row)
+            {
+              for (typename Matrix::const_iterator ptr = matrix.begin(row);
+                   ptr != matrix.end(row);
+                   ++ptr)
+                if (std::abs(ptr->value()) > 0.0 && ptr->column() >= row)
+                  {
+                    a[n_non_zero_elements]   = ptr->value();
+                    irn[n_non_zero_elements] = row + 1;
+                    jcn[n_non_zero_elements] = ptr->column() + 1;
+
+                    ++n_non_zero_elements;
+                  }
+            }
+        }
+      else
+        {
+          for (size_type row = 0; row < matrix.m(); ++row)
+            {
+              for (typename Matrix::const_iterator ptr = matrix.begin(row);
+                   ptr != matrix.end(row);
+                   ++ptr)
+                if (std::abs(ptr->value()) > 0.0)
+                  {
+                    a[n_non_zero_elements]   = ptr->value();
+                    irn[n_non_zero_elements] = row + 1;
+                    jcn[n_non_zero_elements] = ptr->column() + 1;
+                    ++n_non_zero_elements;
+                  }
+            }
         }
-
       id.n   = n;
-      id.nz  = nz;
+      id.nnz = n_non_zero_elements;
       id.irn = irn;
       id.jcn = jcn;
       id.a   = a;
     }
-
-  // No outputs
-  id.icntl[0] = -1;
-  id.icntl[1] = -1;
-  id.icntl[2] = -1;
-  id.icntl[3] = 0;
-
-  // Exit by setting this flag:
-  initialize_called = true;
 }
 
-template <class Matrix>
-void
-SparseDirectMUMPS::initialize(const Matrix         &matrix,
-                              const Vector<double> &vector)
-{
-  // Hand over matrix and right-hand side
-  initialize_matrix(matrix);
 
-  copy_rhs_to_mumps(vector);
-}
 
 void
-SparseDirectMUMPS::copy_rhs_to_mumps(const Vector<double> &new_rhs)
+SparseDirectMUMPS::copy_rhs_to_mumps(const Vector<double> &new_rhs) const
 {
   Assert(n == new_rhs.size(),
          ExcMessage("Matrix size and rhs length must be equal."));
@@ -980,8 +1015,10 @@ SparseDirectMUMPS::copy_rhs_to_mumps(const Vector<double> &new_rhs)
     }
 }
 
+
+
 void
-SparseDirectMUMPS::copy_solution(Vector<double> &vector)
+SparseDirectMUMPS::copy_solution(Vector<double> &vector) const
 {
   Assert(n == vector.size(),
          ExcMessage("Matrix size and solution vector length must be equal."));
@@ -998,48 +1035,52 @@ SparseDirectMUMPS::copy_solution(Vector<double> &vector)
     }
 }
 
+
+
 template <class Matrix>
 void
 SparseDirectMUMPS::initialize(const Matrix &matrix)
 {
   // Initialize MUMPS instance:
   initialize_matrix(matrix);
-  // Start factorization
+
+  // Start analysis + factorization
   id.job = 4;
   dmumps_c(&id);
 }
 
+
+
 void
-SparseDirectMUMPS::solve(Vector<double> &vector)
+SparseDirectMUMPS::vmult(Vector<double> &dst, const Vector<double> &src) const
 {
-  // TODO: this could be implemented similar to SparseDirectUMFPACK where
-  //  the given vector will be used as the RHS. Sadly, there is no easy
-  //  way to do this without breaking the interface.
+  // and that the matrix has at least one nonzero element:
+  Assert(nnz != 0, ExcNotInitialized());
 
-  // Check that the solver has been initialized by the routine above:
-  Assert(initialize_called == true, ExcNotInitialized());
+  Assert(n == dst.size(), ExcMessage("Destination vector has the wrong size."));
+  Assert(n == src.size(), ExcMessage("Source vector has the wrong size."));
 
-  // and that the matrix has at least one nonzero element:
-  Assert(nz != 0, ExcNotInitialized());
+  // Hand over right-hand side
+  copy_rhs_to_mumps(src);
 
   // Start solver
-  id.job = 6; // 6 = analysis, factorization, and solve
+  id.job = 3;
   dmumps_c(&id);
-  copy_solution(vector);
+  copy_solution(dst);
 }
 
+
+
 void
-SparseDirectMUMPS::vmult(Vector<double> &dst, const Vector<double> &src)
+SparseDirectMUMPS::Tvmult(Vector<double> &dst, const Vector<double> &src) const
 {
-  // Check that the solver has been initialized by the routine above:
-  Assert(initialize_called == true, ExcNotInitialized());
-
-  // and that the matrix has at least one nonzero element:
-  Assert(nz != 0, ExcNotInitialized());
+  // The matrix has at least one nonzero element:
+  Assert(nnz != 0, ExcNotInitialized());
 
   Assert(n == dst.size(), ExcMessage("Destination vector has the wrong size."));
   Assert(n == src.size(), ExcMessage("Source vector has the wrong size."));
 
+  id.icntl[8] = 2; // transpose
   // Hand over right-hand side
   copy_rhs_to_mumps(src);
 
@@ -1047,11 +1088,21 @@ SparseDirectMUMPS::vmult(Vector<double> &dst, const Vector<double> &src)
   id.job = 3;
   dmumps_c(&id);
   copy_solution(dst);
+  id.icntl[8] = 1; // reset to default
+}
+
+
+
+int *
+SparseDirectMUMPS::get_icntl()
+{
+  return id.icntl;
 }
 
 #endif // DEAL_II_WITH_MUMPS
 
 
+
 // explicit instantiations for SparseMatrixUMFPACK
 #define InstantiateUMFPACK(MatrixType)                                     \
   template void SparseDirectUMFPACK::factorize(const MatrixType &);        \
@@ -1087,9 +1138,7 @@ InstantiateUMFPACK(BlockSparseMatrix<std::complex<float>>);
 
 // explicit instantiations for SparseDirectMUMPS
 #ifdef DEAL_II_WITH_MUMPS
-#  define InstantiateMUMPS(MATRIX)                                       \
-    template void SparseDirectMUMPS::initialize(const MATRIX &,          \
-                                                const Vector<double> &); \
+#  define InstantiateMUMPS(MATRIX) \
     template void SparseDirectMUMPS::initialize(const MATRIX &);
 
 InstantiateMUMPS(SparseMatrix<double>) InstantiateMUMPS(SparseMatrix<float>)
index 380ef42a74dc96cac7cf76f58ad6520dc8befab4..5420fbf24a50273915896971775eb87798c8912d 100644 (file)
@@ -1,21 +1,19 @@
-// ---------------------------------------------------------------------
+// ------------------------------------------------------------------------
 //
-// Copyright (C) 2002 - 2013 by the deal.II authors
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2001 - 2024 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
 //
-// ---------------------------------------------------------------------
-
+// ------------------------------------------------------------------------
 
 // test the mumps sparse direct solver on a mass matrix
-// - check several ways to solve
+// - check ways to solve
 
 #include <deal.II/base/function.h>
 #include <deal.II/base/quadrature_lib.h>
@@ -56,18 +54,8 @@ solve_and_check(const SparseMatrix<double> &M,
     dst -= solution;
     Assert(dst.l2_norm() < 1e-9, ExcInternalError());
   }
-  {
-    SparseDirectMUMPS solver;
-    solver.initialize(M, rhs);
-    Vector<double> dst(rhs.size());
-    solver.solve(dst);
-    dst -= solution;
-    Assert(dst.l2_norm() < 1e-9, ExcInternalError());
-  }
 }
 
-
-
 template <int dim>
 void
 test()
index bb98ce4852676269a4564d5309dad623802b7005..9a03e0c9694cd6aea2a3b7d0ae2204ac6b5b2778 100644 (file)
@@ -1,25 +1,24 @@
-
 DEAL::1d
 DEAL::Number of dofs = 193
-DEAL::relative norm distance = 2.06357e-16
-DEAL::absolute norms = 6.36406e-13 3084.00
-DEAL::relative norm distance = 2.28670e-16
-DEAL::absolute norms = 1.76304e-12 7709.99
-DEAL::relative norm distance = 2.28670e-16
-DEAL::absolute norms = 3.52609e-12 15420.0
+DEAL::relative norm distance = 2.47865e-16
+DEAL::absolute norms = 7.64416e-13 3084.00
+DEAL::relative norm distance = 2.34248e-16
+DEAL::absolute norms = 1.80605e-12 7709.99
+DEAL::relative norm distance = 2.34248e-16
+DEAL::absolute norms = 3.61211e-12 15420.0
 DEAL::2d
 DEAL::Number of dofs = 1889
-DEAL::relative norm distance = 5.07512e-16
-DEAL::absolute norms = 4.80941e-11 94764.3
-DEAL::relative norm distance = 5.30372e-16
-DEAL::absolute norms = 1.25651e-10 236911.
-DEAL::relative norm distance = 5.30372e-16
-DEAL::absolute norms = 2.51302e-10 473822.
+DEAL::relative norm distance = 4.92318e-16
+DEAL::absolute norms = 4.66541e-11 94764.3
+DEAL::relative norm distance = 4.64003e-16
+DEAL::absolute norms = 1.09927e-10 236911.
+DEAL::relative norm distance = 4.65511e-16
+DEAL::absolute norms = 2.20569e-10 473822.
 DEAL::3d
 DEAL::Number of dofs = 1333
-DEAL::relative norm distance = 1.33240e-15
-DEAL::absolute norms = 7.48349e-11 56165.6
-DEAL::relative norm distance = 1.27330e-15
-DEAL::absolute norms = 1.78789e-10 140414.
-DEAL::relative norm distance = 1.27330e-15
-DEAL::absolute norms = 3.57579e-10 280828.
+DEAL::relative norm distance = 1.21148e-15
+DEAL::absolute norms = 6.80433e-11 56165.6
+DEAL::relative norm distance = 1.43450e-15
+DEAL::absolute norms = 2.01424e-10 140414.
+DEAL::relative norm distance = 1.43535e-15
+DEAL::absolute norms = 4.03086e-10 280828.
\ No newline at end of file
index cda7a2d9d4dc5d99df8c180626a925c2b0cf4377..c6e38f366f4bb8b3142303cca433e1b009b116c9 100644 (file)
@@ -1,17 +1,16 @@
-// ---------------------------------------------------------------------
+// ------------------------------------------------------------------------
 //
-// Copyright (C) 2002 - 2013 by the deal.II authors
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2001 - 2024 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
 //
-// ---------------------------------------------------------------------
+// ------------------------------------------------------------------------
 
 
 // test the mumps sparse direct solver on a mass matrix
@@ -92,8 +91,8 @@ test()
       B.vmult(b, solution);
 
       SparseDirectMUMPS solver;
-      solver.initialize(B, b);
-      solver.solve(x);
+      solver.initialize(B);
+      solver.vmult(x, b);
 
       x -= solution;
       deallog << "relative norm distance = " << x.l2_norm() / solution.l2_norm()
index 8df3f884c1722765a5f173d586e9e6ce9c4f9d8f..bc8c8288bece20e157bd95434cff83b6313a5f39 100644 (file)
@@ -1,25 +1,24 @@
-
 DEAL::1d
 DEAL::Number of dofs = 193
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 3084.00
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 7709.99
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 15420.0
+DEAL::relative norm distance = 2.47865e-16
+DEAL::absolute norms = 7.64416e-13 3084.00
+DEAL::relative norm distance = 2.34248e-16
+DEAL::absolute norms = 1.80605e-12 7709.99
+DEAL::relative norm distance = 2.34248e-16
+DEAL::absolute norms = 3.61211e-12 15420.0
 DEAL::2d
 DEAL::Number of dofs = 1889
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 94764.3
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 236911.
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 473822.
+DEAL::relative norm distance = 4.96258e-16
+DEAL::absolute norms = 4.70275e-11 94764.3
+DEAL::relative norm distance = 4.64464e-16
+DEAL::absolute norms = 1.10036e-10 236911.
+DEAL::relative norm distance = 4.65627e-16
+DEAL::absolute norms = 2.20624e-10 473822.
 DEAL::3d
 DEAL::Number of dofs = 1333
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 56165.6
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 140414.
-DEAL::relative norm distance = 0
-DEAL::absolute norms = 0 280828.
+DEAL::relative norm distance = 1.21940e-15
+DEAL::absolute norms = 6.84883e-11 56165.6
+DEAL::relative norm distance = 1.43458e-15
+DEAL::absolute norms = 2.01434e-10 140414.
+DEAL::relative norm distance = 1.43902e-15
+DEAL::absolute norms = 4.04117e-10 280828.
diff --git a/tests/mumps/mumps_03.cc b/tests/mumps/mumps_03.cc
new file mode 100644 (file)
index 0000000..e85c602
--- /dev/null
@@ -0,0 +1,147 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2001 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+// test the mumps sparse direct solver on a mass matrix
+// test of the transpose option
+// test of other options
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+solve_and_check(const SparseMatrix<double> &M,
+                const Vector<double>       &rhs,
+                const Vector<double>       &solution)
+{
+  {
+    SparseDirectMUMPS::AdditionalData data;
+    SparseDirectMUMPS                 solver(data);
+    data.output_details = false;
+    solver.initialize(M);
+    Vector<double> dst(rhs.size());
+    solver.Tvmult(dst, rhs);
+    dst -= solution;
+    Assert(dst.l2_norm() < 1e-9, ExcInternalError());
+  }
+}
+
+template <int dim>
+void
+test()
+{
+  deallog << dim << 'd' << std::endl;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1);
+  tria.refine_global(1);
+
+  // destroy the uniformity of the matrix by
+  // refining one cell
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(8 - 2 * dim);
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  deallog << "Number of dofs = " << dof_handler.n_dofs() << std::endl;
+
+  SparsityPattern sparsity_pattern;
+  sparsity_pattern.reinit(dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  SparseMatrix<double> B;
+  B.reinit(sparsity_pattern);
+
+  QGauss<dim> qr(2);
+  MatrixTools::create_mass_matrix(dof_handler, qr, B);
+
+  // compute a decomposition of the matrix
+
+  SparseDirectMUMPS::AdditionalData data;
+  data.output_details   = true;
+  data.error_statistics = true;
+  SparseDirectMUMPS Binv(data);
+  Binv.initialize(B);
+
+  // for a number of different solution
+  // vectors, make up a matching rhs vector
+  // and check what the solver finds
+  for (unsigned int i = 0; i < 3; ++i)
+    {
+      Vector<double> solution(dof_handler.n_dofs());
+      Vector<double> x(dof_handler.n_dofs());
+      Vector<double> b(dof_handler.n_dofs());
+
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        solution(j) = j + j * (i + 1) * (i + 1);
+
+      B.Tvmult(b, solution);
+
+      Binv.Tvmult(x, b);
+
+      x -= solution;
+      deallog << "relative norm distance = " << x.l2_norm() / solution.l2_norm()
+              << std::endl;
+      deallog << "absolute norms = " << x.l2_norm() << ' ' << solution.l2_norm()
+              << std::endl;
+      Assert(x.l2_norm() / solution.l2_norm() < 1e-8, ExcInternalError());
+
+
+      // also check solving in different ways:
+      solve_and_check(B, b, solution);
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, numbers::invalid_unsigned_int);
+
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/mumps/mumps_03.output b/tests/mumps/mumps_03.output
new file mode 100644 (file)
index 0000000..09220a8
--- /dev/null
@@ -0,0 +1,24 @@
+DEAL::1d
+DEAL::Number of dofs = 193
+DEAL::relative norm distance = 2.50020e-16
+DEAL::absolute norms = 7.71061e-13 3084.00
+DEAL::relative norm distance = 2.26119e-16
+DEAL::absolute norms = 1.74338e-12 7709.99
+DEAL::relative norm distance = 2.26119e-16
+DEAL::absolute norms = 3.48675e-12 15420.0
+DEAL::2d
+DEAL::Number of dofs = 1889
+DEAL::relative norm distance = 4.62992e-16
+DEAL::absolute norms = 4.38751e-11 94764.3
+DEAL::relative norm distance = 4.44464e-16
+DEAL::absolute norms = 1.05298e-10 236911.
+DEAL::relative norm distance = 4.46851e-16
+DEAL::absolute norms = 2.11727e-10 473822.
+DEAL::3d
+DEAL::Number of dofs = 1333
+DEAL::relative norm distance = 1.27339e-15
+DEAL::absolute norms = 7.15208e-11 56165.6
+DEAL::relative norm distance = 1.13118e-15
+DEAL::absolute norms = 1.58833e-10 140414.
+DEAL::relative norm distance = 1.12371e-15
+DEAL::absolute norms = 3.15569e-10 280828.
\ No newline at end of file
diff --git a/tests/mumps/mumps_04.cc b/tests/mumps/mumps_04.cc
new file mode 100644 (file)
index 0000000..06c2b3d
--- /dev/null
@@ -0,0 +1,162 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2001 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+// test the mumps sparse direct solver on a mass matrix
+// - check several ways to solve, also using block low rank preconditioning
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+solve_and_check(const SparseMatrix<double> &M,
+                const Vector<double>       &rhs,
+                const Vector<double>       &solution)
+{
+  {
+    int                              *icntl = nullptr;
+    SparseDirectMUMPS::AdditionalData data;
+    data.output_details = false;
+    SparseDirectMUMPS solver(data);
+    icntl = solver.get_icntl();
+    std::cout << "Elementi di icntl " << icntl[7] << " " << icntl[8] << " "
+              << icntl[9] << std::endl;
+    solver.initialize(M);
+    Vector<double> dst(rhs.size());
+    solver.vmult(dst, rhs);
+    dst -= solution;
+    Assert(dst.l2_norm() < 1e-9, ExcInternalError());
+  }
+}
+
+
+
+template <int dim>
+void
+test()
+{
+  deallog << dim << 'd' << std::endl;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1);
+  tria.refine_global(1);
+
+  // destroy the uniformity of the matrix by
+  // refining one cell
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(8 - 2 * dim);
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  deallog << "Number of dofs = " << dof_handler.n_dofs() << std::endl;
+
+  SparsityPattern sparsity_pattern;
+  sparsity_pattern.reinit(dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  SparseMatrix<double> B;
+  B.reinit(sparsity_pattern);
+
+  QGauss<dim> qr(2);
+  MatrixTools::create_mass_matrix(dof_handler, qr, B);
+
+  // compute a decomposition of the matrix
+  SparseDirectMUMPS::AdditionalData data;
+  data.blr_factorization = true;
+  SparseDirectMUMPS::AdditionalData::BlockLowRank blr;
+  blr.blr_ucfs          = false;
+  blr.lowrank_threshold = 1e-8;
+  data.output_details   = false;
+  data.error_statistics = false;
+  SparseDirectMUMPS Binv(data);
+
+  Binv.initialize(B);
+
+  // for a number of different solution
+  // vectors, make up a matching rhs vector
+  // and check what the solver finds
+  for (unsigned int i = 0; i < 3; ++i)
+    {
+      Vector<double> solution(dof_handler.n_dofs());
+      Vector<double> x(dof_handler.n_dofs());
+      Vector<double> b(dof_handler.n_dofs());
+
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        solution(j) = j + j * (i + 1) * (i + 1);
+
+      B.vmult(b, solution);
+
+      SolverControl            solver_control(1000, 1e-8);
+      SolverCG<Vector<double>> solver(solver_control);
+
+      solver.solve(B, x, b, Binv);
+
+      deallog << "Number of preconditioned CG iterations = "
+              << solver_control.last_step() << std::endl;
+
+      x -= solution;
+      deallog << "relative norm distance = " << x.l2_norm() / solution.l2_norm()
+              << std::endl;
+      deallog << "absolute norms = " << x.l2_norm() << ' ' << solution.l2_norm()
+              << std::endl;
+      Assert(x.l2_norm() / solution.l2_norm() < 1e-8, ExcInternalError());
+
+      solve_and_check(B, b, solution);
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, numbers::invalid_unsigned_int);
+
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/mumps/mumps_04.output b/tests/mumps/mumps_04.output
new file mode 100644 (file)
index 0000000..d0af0b4
--- /dev/null
@@ -0,0 +1,33 @@
+DEAL::1d
+DEAL::Number of dofs = 193
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 2.47865e-16
+DEAL::absolute norms = 7.64416e-13 3084.00
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 2.34248e-16
+DEAL::absolute norms = 1.80605e-12 7709.99
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 2.34248e-16
+DEAL::absolute norms = 3.61211e-12 15420.0
+DEAL::2d
+DEAL::Number of dofs = 1889
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 5.48214e-16
+DEAL::absolute norms = 5.19511e-11 94764.3
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 4.62197e-16
+DEAL::absolute norms = 1.09499e-10 236911.
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 4.68810e-16
+DEAL::absolute norms = 2.22132e-10 473822.
+DEAL::3d
+DEAL::Number of dofs = 1333
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 1.21396e-15
+DEAL::absolute norms = 6.81827e-11 56165.6
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 1.43949e-15
+DEAL::absolute norms = 2.02124e-10 140414.
+DEAL::Number of preconditioned CG iterations = 1
+DEAL::relative norm distance = 1.44178e-15
+DEAL::absolute norms = 4.04892e-10 280828.
\ No newline at end of file
diff --git a/tests/mumps/mumps_05.cc b/tests/mumps/mumps_05.cc
new file mode 100644 (file)
index 0000000..97802a0
--- /dev/null
@@ -0,0 +1,145 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2001 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+// test the mumps sparse direct solver on a mass matrix
+// tests of all the symmetric matrix options
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+void
+solve_and_check(const SparseMatrix<double> &M,
+                const Vector<double>       &rhs,
+                const Vector<double>       &solution)
+{
+  SparseDirectMUMPS::AdditionalData data;
+  data.output_details = false;
+  data.symmetric      = true;
+  data.posdef         = true;
+  SparseDirectMUMPS solver(data);
+  solver.initialize(M);
+  Vector<double> dst(rhs.size());
+  solver.vmult(dst, rhs);
+  dst -= solution;
+  Assert(dst.l2_norm() < 1e-8, ExcInternalError());
+}
+
+template <int dim>
+void
+test()
+{
+  deallog << dim << 'd' << std::endl;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, -1, 1);
+  tria.refine_global(1);
+
+  // destroy the uniformity of the matrix by
+  // refining one cell
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(8 - 2 * dim);
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  deallog << "Number of dofs = " << dof_handler.n_dofs() << std::endl;
+
+  SparsityPattern sparsity_pattern;
+  sparsity_pattern.reinit(dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  SparseMatrix<double> B;
+  B.reinit(sparsity_pattern);
+
+  QGauss<dim> qr(2);
+  MatrixTools::create_mass_matrix(dof_handler, qr, B);
+
+  // compute a decomposition of the matrix
+
+  SparseDirectMUMPS::AdditionalData data;
+  data.output_details   = true;
+  data.error_statistics = true;
+  data.symmetric        = true;
+  data.posdef           = false;
+  SparseDirectMUMPS Binv(data);
+  Binv.initialize(B);
+
+  // for a number of different solution
+  // vectors, make up a matching rhs vector
+  // and check what the solver finds
+  for (unsigned int i = 0; i < 3; ++i)
+    {
+      Vector<double> solution(dof_handler.n_dofs());
+      Vector<double> x(dof_handler.n_dofs());
+      Vector<double> b(dof_handler.n_dofs());
+
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        solution(j) = j + j * (i + 1) * (i + 1);
+
+      B.vmult(b, solution);
+
+      Binv.vmult(x, b);
+
+      x -= solution;
+      deallog << "relative norm distance = " << x.l2_norm() / solution.l2_norm()
+              << std::endl;
+      deallog << "absolute norms = " << x.l2_norm() << ' ' << solution.l2_norm()
+              << std::endl;
+      Assert(x.l2_norm() / solution.l2_norm() < 1e-8, ExcInternalError());
+
+      // check with also the posdef option
+      solve_and_check(B, b, solution);
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, numbers::invalid_unsigned_int);
+
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/mumps/mumps_05.output b/tests/mumps/mumps_05.output
new file mode 100644 (file)
index 0000000..3fd778e
--- /dev/null
@@ -0,0 +1,24 @@
+DEAL::1d
+DEAL::Number of dofs = 193
+DEAL::relative norm distance = 2.44327e-16
+DEAL::absolute norms = 7.53504e-13 3084.00
+DEAL::relative norm distance = 2.25536e-16
+DEAL::absolute norms = 1.73888e-12 7709.99
+DEAL::relative norm distance = 2.25536e-16
+DEAL::absolute norms = 3.47776e-12 15420.0
+DEAL::2d
+DEAL::Number of dofs = 1889
+DEAL::relative norm distance = 5.02782e-16
+DEAL::absolute norms = 4.76458e-11 94764.3
+DEAL::relative norm distance = 4.72365e-16
+DEAL::absolute norms = 1.11908e-10 236911.
+DEAL::relative norm distance = 4.74917e-16
+DEAL::absolute norms = 2.25026e-10 473822.
+DEAL::3d
+DEAL::Number of dofs = 1333
+DEAL::relative norm distance = 1.29200e-15
+DEAL::absolute norms = 7.25659e-11 56165.6
+DEAL::relative norm distance = 1.42228e-15
+DEAL::absolute norms = 1.99708e-10 140414.
+DEAL::relative norm distance = 1.42006e-15
+DEAL::absolute norms = 3.98792e-10 280828.
\ No newline at end of file

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.