]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use directly petsc and trilinos interfaces 18497/head
authorMarco Feder <marco.feder@sissa.it>
Sat, 21 Jun 2025 17:04:37 +0000 (19:04 +0200)
committerMarco Feder <marco.feder@sissa.it>
Sat, 21 Jun 2025 17:04:37 +0000 (19:04 +0200)
source/lac/sparse_direct.cc

index 3bbf2ceb9262717eaca6b4e6b2c90b1309cc73f9..75e46da54b105fba109915acff5ddd7fbe641a68 100644 (file)
@@ -1040,52 +1040,171 @@ SparseDirectMUMPS::initialize_matrix(const Matrix &matrix)
 
       if (additional_data.symmetric == true)
         {
-          for (const auto &row : locally_owned_rows)
+          if constexpr (std::is_same_v<Matrix,
+                                       PETScWrappers::MPI::SparseMatrix>)
             {
-              for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
-                if (std::abs(it->value()) > 0.0 && it->column() >= row)
-                  {
-                    const int global_row = row + 1;
-                    const int global_col = it->column() + 1;
-
-                    // Store this non-zero entry
-                    irn[n_non_zero_local] = global_row;
-                    jcn[n_non_zero_local] = global_col;
-                    a[n_non_zero_local]   = it->value();
-
-                    // Count local non-zeros
-                    n_non_zero_local++;
-                  }
-
-              const types::global_cell_index local_index =
-                locally_owned_rows.index_within_set(row);
-              irhs_loc[local_index] = row + 1;
+              Mat &petsc_matrix =
+                const_cast<PETScWrappers::MPI::SparseMatrix &>(matrix)
+                  .petsc_matrix();
+
+              PetscInt rstart, rend;
+              MatGetOwnershipRange(petsc_matrix, &rstart, &rend);
+              for (PetscInt i = rstart; i < rend; i++)
+                {
+                  PetscInt           n_cols;
+                  const PetscInt    *cols;
+                  const PetscScalar *values;
+                  MatGetRow(petsc_matrix, i, &n_cols, &cols, &values);
+
+                  for (PetscInt j = 0; j < n_cols; j++)
+                    {
+                      if (cols[j] >= i)
+                        {
+                          irn[n_non_zero_local] = i + 1;
+                          jcn[n_non_zero_local] = cols[j] + 1;
+                          a[n_non_zero_local]   = values[j];
+
+                          // Count local non-zeros
+                          n_non_zero_local++;
+                        }
+                    }
+                  MatRestoreRow(petsc_matrix, i, &n_cols, &cols, &values);
+
+                  // Store the row index for the rhs vector
+                  const types::global_cell_index local_index =
+                    locally_owned_rows.index_within_set(i);
+                  irhs_loc[local_index] = i + 1;
+                }
+
+              id.a_loc = a.get();
+            }
+          else if constexpr (std::is_same_v<Matrix,
+                                            TrilinosWrappers::SparseMatrix>)
+            {
+              const auto &trilinos_matrix = matrix.trilinos_matrix();
+              for (int local_row = 0; local_row < trilinos_matrix.NumMyRows();
+                   ++local_row)
+                {
+                  int     num_entries;
+                  double *values;
+                  int    *local_cols;
+                  int     ierr = trilinos_matrix.ExtractMyRowView(local_row,
+                                                              num_entries,
+                                                              values,
+                                                              local_cols);
+                  (void)ierr;
+                  Assert(
+                    ierr == 0,
+                    ExcMessage(
+                      "Error extracting global row view from Trilinos matrix. Error code " +
+                      std::to_string(ierr) + "."));
+
+                  int global_row = trilinos_matrix.GRID(local_row);
+
+                  for (int j = 0; j < num_entries; ++j)
+                    {
+                      if (trilinos_matrix.GCID(local_cols[j]) >= global_row)
+                        {
+                          irn[n_non_zero_local] = global_row + 1;
+                          jcn[n_non_zero_local] =
+                            trilinos_matrix.GCID(local_cols[j]) + 1;
+                          a[n_non_zero_local] = values[j];
+
+                          // Count local non-zeros
+                          n_non_zero_local++;
+                        }
+                    }
+
+                  // Store the row index for the rhs vector
+                  irhs_loc[local_row] = global_row + 1;
+                }
+              id.a_loc = a.get();
+            }
+          else
+            {
+              DEAL_II_NOT_IMPLEMENTED();
             }
         }
       else
         {
-          for (const auto &row : locally_owned_rows)
+          // Unsymmetric case
+          if constexpr (std::is_same_v<Matrix,
+                                       PETScWrappers::MPI::SparseMatrix>)
             {
-              // Loop over columns
-              for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
+              Mat &petsc_matrix =
+                const_cast<PETScWrappers::MPI::SparseMatrix &>(matrix)
+                  .petsc_matrix();
+
+              PetscInt rstart, rend;
+              MatGetOwnershipRange(petsc_matrix, &rstart, &rend);
+              for (PetscInt i = rstart; i < rend; i++)
                 {
-                  if (std::abs(it->value()) > 0.0)
-                    {
-                      const int global_row = row + 1;
-                      const int global_col = it->column() + 1;
+                  PetscInt           n_cols;
+                  const PetscInt    *cols;
+                  const PetscScalar *values;
+                  MatGetRow(petsc_matrix, i, &n_cols, &cols, &values);
 
-                      irn[n_non_zero_local] = global_row;
-                      jcn[n_non_zero_local] = global_col;
-                      a[n_non_zero_local]   = it->value();
+                  for (PetscInt j = 0; j < n_cols; j++)
+                    {
+                      irn[n_non_zero_local] = i + 1;
+                      jcn[n_non_zero_local] = cols[j] + 1;
+                      a[n_non_zero_local]   = values[j];
 
                       // Count local non-zeros
                       n_non_zero_local++;
                     }
+                  MatRestoreRow(petsc_matrix, i, &n_cols, &cols, &values);
+
+                  // Store the row index for the rhs vector
+                  const types::global_cell_index local_index =
+                    locally_owned_rows.index_within_set(i);
+                  irhs_loc[local_index] = i + 1;
                 }
 
-              const types::global_cell_index local_index =
-                locally_owned_rows.index_within_set(row);
-              irhs_loc[local_index] = row + 1;
+              id.a_loc = a.get();
+            }
+          else if constexpr (std::is_same_v<Matrix,
+                                            TrilinosWrappers::SparseMatrix>)
+            {
+              const auto &trilinos_matrix = matrix.trilinos_matrix();
+              for (int local_row = 0; local_row < trilinos_matrix.NumMyRows();
+                   ++local_row)
+                {
+                  int     num_entries;
+                  double *values;
+                  int    *local_cols;
+                  int     ierr = trilinos_matrix.ExtractMyRowView(local_row,
+                                                              num_entries,
+                                                              values,
+                                                              local_cols);
+                  (void)ierr;
+                  Assert(
+                    ierr == 0,
+                    ExcMessage(
+                      "Error extracting global row view from Trilinos matrix. Error code " +
+                      std::to_string(ierr) + "."));
+
+                  int global_row = trilinos_matrix.GRID(local_row);
+
+                  for (int j = 0; j < num_entries; ++j)
+                    {
+                      irn[n_non_zero_local] = global_row + 1;
+                      jcn[n_non_zero_local] =
+                        trilinos_matrix.GCID(local_cols[j]) + 1;
+                      a[n_non_zero_local] = values[j];
+
+                      // Count local non-zeros
+                      n_non_zero_local++;
+                    }
+
+                  // Store the row index for the rhs vector
+                  irhs_loc[local_row] = global_row + 1;
+                }
+              id.a_loc = a.get();
+            }
+          else
+            {
+              DEAL_II_NOT_IMPLEMENTED();
             }
         }
 

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.