--- /dev/null
+Improved: The SparseDirectMUMPS class now takes an AdditionalData to
+control the MUMPS execution.
+<br>
+(Davide Polverino, 2025/05/20)
*/
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
#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()
{
}
}
+
+
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."));
}
}
+
+
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."));
}
}
+
+
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);
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 &); \
// 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>)
-// ---------------------------------------------------------------------
+// ------------------------------------------------------------------------
//
-// 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>
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()
-
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
-// ---------------------------------------------------------------------
+// ------------------------------------------------------------------------
//
-// 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
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()
-
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.
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+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