]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Address some comments
authorMarco Feder <marco.feder@sissa.it>
Fri, 6 Jun 2025 12:21:44 +0000 (14:21 +0200)
committerMarco Feder <marco.feder@sissa.it>
Fri, 6 Jun 2025 12:21:44 +0000 (14:21 +0200)
include/deal.II/lac/sparse_direct.h
source/lac/sparse_direct.cc
tests/mumps/mumps_06.cc

index 6f0f46232ea04814bffef06ee23b5e3ab0370640..2a75a0493cfce17178577501f9d968fa0e7e5ef4 100644 (file)
@@ -526,12 +526,12 @@ public:
     {}
 
     /**
-     * Print details of the execution
+     * If true, the MUMPS solver will print out details of the execution.
      */
     bool output_details;
 
     /**
-     * Print error statistics
+     * If true, the MUMPS solver will print out error statistics.
      */
     bool error_statistics;
 
@@ -672,7 +672,7 @@ private:
   /**
    * IndexSet storing the locally owned rows of the matrix.
    */
-  IndexSet row_range;
+  IndexSet locally_owned_rows;
 
   /**
    * This function initializes a MUMPS instance and hands over the system's
index 2d02cd554b5c52d4c7d33b69b05a4ec6cdf6652c..3bbf2ceb9262717eaca6b4e6b2c90b1309cc73f9 100644 (file)
@@ -931,15 +931,14 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
   n    = matrix.n();
   id.n = n;
 
-  if constexpr (std::is_same_v<SparseMatrix<double>, Matrix>)
+  if constexpr (std::is_same_v<Matrix, SparseMatrix<double>>)
     {
       // Serial matrix: hand over matrix to MUMPS as centralized assembled
       // matrix
       if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
         {
           // number of nonzero elements in matrix
-          if constexpr (std::is_same_v<Matrix, SparseMatrix<double>>)
-            nnz = matrix.n_actually_nonzero_elements();
+          nnz = matrix.n_actually_nonzero_elements();
 
           // representation of the matrix
           a = std::make_unique<double[]>(nnz);
@@ -994,9 +993,17 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
           id.a   = a.get();
         }
     }
-  else if constexpr (std::is_same_v<TrilinosWrappers::SparseMatrix, Matrix> ||
-                     std::is_same_v<PETScWrappers::MPI::SparseMatrix, Matrix>)
+  else if constexpr (std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix> ||
+                     std::is_same_v<Matrix, PETScWrappers::MPI::SparseMatrix>)
     {
+      int result;
+      MPI_Comm_compare(mpi_communicator,
+                       matrix.get_mpi_communicator(),
+                       &result);
+      AssertThrow(result == MPI_IDENT,
+                  ExcMessage("The matrix communicator must match the MUMPS "
+                             "communicator."));
+
       // Distributed matrix case
       id.icntl[17]               = 3; // distributed matrix assembly
       id.nnz                     = matrix.n_nonzero_elements();
@@ -1004,34 +1011,36 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
       size_type n_non_zero_local = 0;
 
       // Get the range of rows owned by this process
-      row_range                 = matrix.locally_owned_range_indices();
+      locally_owned_rows        = matrix.locally_owned_range_indices();
       size_type local_non_zeros = 0;
 
-      if constexpr (std::is_same_v<TrilinosWrappers::SparseMatrix, Matrix>)
+      if constexpr (std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix>)
         {
           const auto &trilinos_matrix = matrix.trilinos_matrix();
           local_non_zeros             = trilinos_matrix.NumMyNonzeros();
         }
-      else if constexpr (std::is_same_v<PETScWrappers::MPI::SparseMatrix,
-                                        Matrix>)
+      else if constexpr (std::is_same_v<Matrix,
+                                        PETScWrappers::MPI::SparseMatrix>)
         {
           Mat &petsc_matrix =
             const_cast<PETScWrappers::MPI::SparseMatrix &>(matrix)
               .petsc_matrix();
           MatInfo info;
           MatGetInfo(petsc_matrix, MAT_LOCAL, &info);
-          local_non_zeros = (PetscInt)info.nz_used;
+          local_non_zeros = (size_type)info.nz_used;
         }
 
 
+      // We allocate enough entries for the general, nonsymmetric, case. If case
+      // of a symmetric matrix, we will end up with fewer entries.
       irn = std::make_unique<MUMPS_INT[]>(local_non_zeros);
       jcn = std::make_unique<MUMPS_INT[]>(local_non_zeros);
       a   = std::make_unique<double[]>(local_non_zeros);
-      irhs_loc.resize(row_range.n_elements());
+      irhs_loc.resize(locally_owned_rows.n_elements());
 
       if (additional_data.symmetric == true)
         {
-          for (const auto &row : row_range)
+          for (const auto &row : locally_owned_rows)
             {
               for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
                 if (std::abs(it->value()) > 0.0 && it->column() >= row)
@@ -1049,13 +1058,13 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
                   }
 
               const types::global_cell_index local_index =
-                row_range.index_within_set(row);
+                locally_owned_rows.index_within_set(row);
               irhs_loc[local_index] = row + 1;
             }
         }
       else
         {
-          for (const auto &row : row_range)
+          for (const auto &row : locally_owned_rows)
             {
               // Loop over columns
               for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
@@ -1075,7 +1084,7 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
                 }
 
               const types::global_cell_index local_index =
-                row_range.index_within_set(row);
+                locally_owned_rows.index_within_set(row);
               irhs_loc[local_index] = row + 1;
             }
         }
@@ -1092,7 +1101,7 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
       id.icntl[20] = 0;  // centralized solution, stored on rank 0 by MUMPS
       id.nrhs      = 1;
       id.lrhs_loc  = n;
-      id.nloc_rhs  = row_range.n_elements();
+      id.nloc_rhs  = locally_owned_rows.n_elements();
     }
   else
     {
@@ -1220,6 +1229,10 @@ SparseDirectMUMPS::vmult(VectorType &dst, const VectorType &src) const
 
       rhs.resize(0); // remove rhs again
     }
+  else
+    {
+      DEAL_II_NOT_IMPLEMENTED();
+    }
 }
 
 
index b07422dee47a43cfeedd6afeb9d6c19f40a6952b..d899174489b3e1665732b163400878d6ddadea01 100644 (file)
@@ -56,7 +56,7 @@ solve_and_check(const MatrixType &M,
   data.output_details = false;
   data.symmetric      = true;
   data.posdef         = true;
-  SparseDirectMUMPS solver(data);
+  SparseDirectMUMPS solver(data, M.get_mpi_communicator());
   solver.initialize(M);
   VectorType dst(rhs);
   solver.vmult(dst, rhs);
@@ -110,7 +110,7 @@ test()
   SparseDirectMUMPS::AdditionalData data;
   data.output_details   = true;
   data.error_statistics = true;
-  SparseDirectMUMPS Binv(data);
+  SparseDirectMUMPS Binv(data, B.get_mpi_communicator());
   Binv.initialize(B);
 
   // for a number of different solution

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.