]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid raw pointers in SparseDirectMUMPS
authorDavide Polverino <davide.polverino@phd.unipi.it>
Tue, 20 May 2025 13:30:26 +0000 (15:30 +0200)
committerMarco Feder <marco.feder@sissa.it>
Fri, 23 May 2025 12:25:04 +0000 (14:25 +0200)
include/deal.II/lac/sparse_direct.h
source/lac/sparse_direct.cc
tests/mumps/mumps_03.output
tests/mumps/mumps_04.output

index af369422c8406e484b8fa32a783d90e80d45fa84..236230f1ae19428485eb4ba423bf2f051843f86c 100644 (file)
@@ -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<double[]> 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<int[]> irn;
 
   /**
    * jcn contains the column indices of the non-zero entries of the matrix.
    */
-  int *jcn;
+  std::unique_ptr<int[]> jcn;
 
   /**
    * The number of rows of the matrix. The matrix is square.
index 48f9d4763356c72f2ca9851a04867ec60152a274..f4d6e6377f35c6b2f6edba82d0a4fba19fa160b2 100644 (file)
@@ -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<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[nnz];
-      jcn = new int[nnz];
+      irn = std::make_unique<int[]>(nnz);
+      jcn = std::make_unique<int[]>(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<double> &dst, const Vector<double> &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."));
 
index 09220a85268d17c4153972ccd91913aed6788522..08e46485ff0e29c580bf18848de84859929a0b92 100644 (file)
@@ -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.
index d0af0b4893c9f0e3b89862a66b6f90c049afe744..c6c94f2d4acdb61aaea42a240ca434af039592ae 100644 (file)
@@ -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.

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.