From c997cb94da8336b80941ccf6019f79c536955dff Mon Sep 17 00:00:00 2001 From: Davide Polverino Date: Tue, 20 May 2025 15:30:26 +0200 Subject: [PATCH] Avoid raw pointers in SparseDirectMUMPS --- include/deal.II/lac/sparse_direct.h | 10 ++++------ source/lac/sparse_direct.cc | 29 +++++++++++------------------ tests/mumps/mumps_03.output | 2 +- tests/mumps/mumps_04.output | 2 +- 4 files changed, 17 insertions(+), 26 deletions(-) diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index af369422c8..236230f1ae 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -458,6 +458,7 @@ public: */ using size_type = types::global_dof_index; + /** * The AdditionalData struct contains data for controlling the MUMPS * execution. @@ -480,7 +481,6 @@ public: * algorithm, an algorithm with higher compression than the standard one. */ bool blr_ucfs; - /** * Threshold for the low-rank truncation of the blocks. */ @@ -510,13 +510,11 @@ public: * 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. @@ -598,7 +596,7 @@ private: * 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. */ - double *a; + std::unique_ptr a; /** * The right-hand side vector. This is the right hand side of the system. @@ -608,12 +606,12 @@ private: /** * irn contains the row indices of the non-zero entries of the matrix. */ - int *irn; + std::unique_ptr irn; /** * jcn contains the column indices of the non-zero entries of the matrix. */ - int *jcn; + std::unique_ptr jcn; /** * The number of rows of the matrix. The matrix is square. diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index 48f9d47633..f4d6e6377f 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -884,7 +884,7 @@ SparseDirectMUMPS::SparseDirectMUMPS(const AdditionalData &data) id.sym = 0; // Use MPI_COMM_WORLD as communicator - id.comm_fortran = -987654; + id.comm_fortran = (MUMPS_INT)MPI_Comm_c2f(MPI_COMM_WORLD); dmumps_c(&id); if (additional_data.output_details == false) @@ -914,16 +914,9 @@ SparseDirectMUMPS::SparseDirectMUMPS(const AdditionalData &data) SparseDirectMUMPS::~SparseDirectMUMPS() { + // MUMPS destructor id.job = -2; dmumps_c(&id); - - // Do some cleaning - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - { - delete[] a; - delete[] irn; - delete[] jcn; - } } @@ -933,6 +926,9 @@ void SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) { Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square.")); + // Here we should be checking if the matrix is respecting the symmetry given + //(I.E. sym = 0 for non-symmetric matrix, sym = 1 for posdef matrix, sym = 2 + // for general symmetric). // Hand over matrix to MUMPS as centralized assembled matrix if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) @@ -944,13 +940,13 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) nnz = matrix.n_actually_nonzero_elements(); // representation of the matrix - a = new double[nnz]; + a = std::make_unique(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[nnz]; - jcn = new int[nnz]; + irn = std::make_unique(nnz); + jcn = std::make_unique(nnz); size_type n_non_zero_elements = 0; @@ -991,9 +987,9 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix) } id.n = n; id.nnz = n_non_zero_elements; - id.irn = irn; - id.jcn = jcn; - id.a = a; + id.irn = irn.get(); + id.jcn = jcn.get(); + id.a = a.get(); } } @@ -1049,14 +1045,11 @@ SparseDirectMUMPS::initialize(const Matrix &matrix) dmumps_c(&id); } - - void SparseDirectMUMPS::vmult(Vector &dst, const Vector &src) const { // and that 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.")); diff --git a/tests/mumps/mumps_03.output b/tests/mumps/mumps_03.output index 09220a8526..08e46485ff 100644 --- a/tests/mumps/mumps_03.output +++ b/tests/mumps/mumps_03.output @@ -21,4 +21,4 @@ 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 +DEAL::absolute norms = 3.15569e-10 280828. diff --git a/tests/mumps/mumps_04.output b/tests/mumps/mumps_04.output index d0af0b4893..c6c94f2d4a 100644 --- a/tests/mumps/mumps_04.output +++ b/tests/mumps/mumps_04.output @@ -30,4 +30,4 @@ 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 +DEAL::absolute norms = 4.04892e-10 280828. -- 2.39.5