From: Davide Polverino Date: Tue, 20 May 2025 13:30:26 +0000 (+0200) Subject: Add additional data and tests final commit X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F18466%2Fhead;p=dealii.git Add additional data and tests final commit 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 --- diff --git a/doc/news/changes/minor/20250520Polverino b/doc/news/changes/minor/20250520Polverino new file mode 100644 index 0000000000..f21ef769da --- /dev/null +++ b/doc/news/changes/minor/20250520Polverino @@ -0,0 +1,4 @@ +Improved: The SparseDirectMUMPS class now takes an AdditionalData to +control the MUMPS execution. +
+(Davide Polverino, 2025/05/20) diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index c3f0c7df76..af369422c8 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -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 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 AdditionalData 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 matrix. + * Exception. + */ + DeclException0(ExcInitializeAlreadyCalled); + + /** + * This function computes the LU factorization + * of the system's matrix matrix with the options given in the + * constructor's AdditionalData. */ template 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 src and the solution is written into the output vector + * dst. */ void - copy_solution(Vector &vector); + vmult(Vector &dst, const Vector &src) const; /** - * + * A function in which the inverse of the transposed matrix is applied to the + * input vector src and the solution is written into the output + * vector dst. */ void - copy_rhs_to_mumps(const Vector &rhs); + Tvmult(Vector &dst, const Vector &src) const; /** - * Flags storing whether the function initialize () 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 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 matrix and right-hand side rhs_vector to the - * solver. + * The number of rows of the matrix. The matrix is square. */ - template - void - initialize(const Matrix &matrix, const Vector &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 matrix. + * This function hands over to MUMPS the system's matrix. */ template 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 vector. The right-hand side need to be - * supplied in initialize(matrix, vector); + * Copy the computed solution into the solution vector. */ void - solve(Vector &vector); + copy_solution(Vector &vector) const; /** - * A function in which the inverse of the matrix is applied to the input - * vector src and the solution is written into the output vector - * dst. + * Copy the right-hand side vector into the MUMPS instance. */ void - vmult(Vector &dst, const Vector &src); + copy_rhs_to_mumps(const Vector &rhs) const; + + /** + * Struct that holds the additional data for the MUMPS solver. + */ + AdditionalData additional_data; }; DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index f59d6a9316..48f9d47633 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -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 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 -void -SparseDirectMUMPS::initialize(const Matrix &matrix, - const Vector &vector) -{ - // Hand over matrix and right-hand side - initialize_matrix(matrix); - copy_rhs_to_mumps(vector); -} void -SparseDirectMUMPS::copy_rhs_to_mumps(const Vector &new_rhs) +SparseDirectMUMPS::copy_rhs_to_mumps(const Vector &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 &new_rhs) } } + + void -SparseDirectMUMPS::copy_solution(Vector &vector) +SparseDirectMUMPS::copy_solution(Vector &vector) const { Assert(n == vector.size(), ExcMessage("Matrix size and solution vector length must be equal.")); @@ -998,48 +1035,52 @@ SparseDirectMUMPS::copy_solution(Vector &vector) } } + + template 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 &vector) +SparseDirectMUMPS::vmult(Vector &dst, const Vector &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 &dst, const Vector &src) +SparseDirectMUMPS::Tvmult(Vector &dst, const Vector &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 &dst, const Vector &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>); // explicit instantiations for SparseDirectMUMPS #ifdef DEAL_II_WITH_MUMPS -# define InstantiateMUMPS(MATRIX) \ - template void SparseDirectMUMPS::initialize(const MATRIX &, \ - const Vector &); \ +# define InstantiateMUMPS(MATRIX) \ template void SparseDirectMUMPS::initialize(const MATRIX &); InstantiateMUMPS(SparseMatrix) InstantiateMUMPS(SparseMatrix) diff --git a/tests/mumps/mumps_01.cc b/tests/mumps/mumps_01.cc index 380ef42a74..5420fbf24a 100644 --- a/tests/mumps/mumps_01.cc +++ b/tests/mumps/mumps_01.cc @@ -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 #include @@ -56,18 +54,8 @@ solve_and_check(const SparseMatrix &M, dst -= solution; Assert(dst.l2_norm() < 1e-9, ExcInternalError()); } - { - SparseDirectMUMPS solver; - solver.initialize(M, rhs); - Vector dst(rhs.size()); - solver.solve(dst); - dst -= solution; - Assert(dst.l2_norm() < 1e-9, ExcInternalError()); - } } - - template void test() diff --git a/tests/mumps/mumps_01.output b/tests/mumps/mumps_01.output index bb98ce4852..9a03e0c969 100644 --- a/tests/mumps/mumps_01.output +++ b/tests/mumps/mumps_01.output @@ -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 diff --git a/tests/mumps/mumps_02.cc b/tests/mumps/mumps_02.cc index cda7a2d9d4..c6e38f366f 100644 --- a/tests/mumps/mumps_02.cc +++ b/tests/mumps/mumps_02.cc @@ -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() diff --git a/tests/mumps/mumps_02.output b/tests/mumps/mumps_02.output index 8df3f884c1..bc8c8288be 100644 --- a/tests/mumps/mumps_02.output +++ b/tests/mumps/mumps_02.output @@ -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 index 0000000000..e85c602ae0 --- /dev/null +++ b/tests/mumps/mumps_03.cc @@ -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 +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + + +void +solve_and_check(const SparseMatrix &M, + const Vector &rhs, + const Vector &solution) +{ + { + SparseDirectMUMPS::AdditionalData data; + SparseDirectMUMPS solver(data); + data.output_details = false; + solver.initialize(M); + Vector dst(rhs.size()); + solver.Tvmult(dst, rhs); + dst -= solution; + Assert(dst.l2_norm() < 1e-9, ExcInternalError()); + } +} + +template +void +test() +{ + deallog << dim << 'd' << std::endl; + + Triangulation 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 fe(1); + DoFHandler 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 B; + B.reinit(sparsity_pattern); + + QGauss 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 solution(dof_handler.n_dofs()); + Vector x(dof_handler.n_dofs()); + Vector 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 index 0000000000..09220a8526 --- /dev/null +++ b/tests/mumps/mumps_03.output @@ -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 index 0000000000..06c2b3d343 --- /dev/null +++ b/tests/mumps/mumps_04.cc @@ -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 +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + + +void +solve_and_check(const SparseMatrix &M, + const Vector &rhs, + const Vector &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 dst(rhs.size()); + solver.vmult(dst, rhs); + dst -= solution; + Assert(dst.l2_norm() < 1e-9, ExcInternalError()); + } +} + + + +template +void +test() +{ + deallog << dim << 'd' << std::endl; + + Triangulation 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 fe(1); + DoFHandler 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 B; + B.reinit(sparsity_pattern); + + QGauss 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 solution(dof_handler.n_dofs()); + Vector x(dof_handler.n_dofs()); + Vector 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> 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 index 0000000000..d0af0b4893 --- /dev/null +++ b/tests/mumps/mumps_04.output @@ -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 index 0000000000..97802a03ec --- /dev/null +++ b/tests/mumps/mumps_05.cc @@ -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 +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + +void +solve_and_check(const SparseMatrix &M, + const Vector &rhs, + const Vector &solution) +{ + SparseDirectMUMPS::AdditionalData data; + data.output_details = false; + data.symmetric = true; + data.posdef = true; + SparseDirectMUMPS solver(data); + solver.initialize(M); + Vector dst(rhs.size()); + solver.vmult(dst, rhs); + dst -= solution; + Assert(dst.l2_norm() < 1e-8, ExcInternalError()); +} + +template +void +test() +{ + deallog << dim << 'd' << std::endl; + + Triangulation 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 fe(1); + DoFHandler 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 B; + B.reinit(sparsity_pattern); + + QGauss 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 solution(dof_handler.n_dofs()); + Vector x(dof_handler.n_dofs()); + Vector 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 index 0000000000..3fd778ead1 --- /dev/null +++ b/tests/mumps/mumps_05.output @@ -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