Merge pull request arcaneframework#313 from arcaneframework/dev/gg-separate-matrix-and-rhs-transformation

Separate Matrix and RHS transformation for `AlephDoFLinearSystemImpl` et `HypreDoFLinearSystemImpl`

Author: grospelliergilles

femutils/AlephDoFLinearSystem.cc

@@ -199,8 +199,9 @@ class AlephDoFLinearSystemImpl
  private:
 
   AlephParams* _createAlephParam() const;
-  void _fillMatrix();
+  void _applyMatrixTransformationAndFillAlephMatrix();
   void _fillRHSVector();
+  void _applyRHSTransformationAndFillAlephRHS();
   void _setMatrixValue(DoF row, DoF column, Real value)
   {
     if (m_do_print_filling)
@@ -265,7 +266,6 @@ _fillRowColumnEliminationInfos()
 
   auto& dof_elimination_info = getEliminationInfo();
   auto& dof_elimination_value = getEliminationValue();
-  auto& rhs_variable = rhsVariable();
 
   for (const auto& rc_value : m_values_map) {
     RowColumn rc = rc_value.first;
@@ -283,7 +283,7 @@ _fillRowColumnEliminationInfos()
 /*---------------------------------------------------------------------------*/
 
 void AlephDoFLinearSystemImpl::
-_fillMatrix()
+_applyMatrixTransformationAndFillAlephMatrix()
 {
   _fillRowColumnEliminationInfos();
   OrderedRowColumnMap& rc_elimination_map = _rowColumnEliminationMap();
@@ -293,7 +293,6 @@ _fillMatrix()
 
   auto& dof_elimination_info = getEliminationInfo();
   auto& dof_elimination_value = getEliminationValue();
-  auto& rhs_variable = rhsVariable();
 
   // Fill the matrix from the values of \a m_values_map
   // Skip (row,column) values which are part of an elimination.
@@ -322,14 +321,44 @@ _fillMatrix()
     _setMatrixValue(dof_row, dof_column, value);
   }
 
+  const bool do_print_filling = m_do_print_filling;
+
+  // Apply Row or Row+Column elimination on Matrix
+  // Phase 2: set the diagonal value for elimination row to 1.0
+  ENUMERATE_ (DoF, idof, dof_family->allItems()) {
+    DoF dof = *idof;
+    if (!dof.isOwn())
+      continue;
+    Byte elimination_info = dof_elimination_info[dof];
+    if (elimination_info == ELIMINATE_ROW || elimination_info == ELIMINATE_ROW_COLUMN) {
+      Real elimination_value = dof_elimination_value[dof];
+      if (do_print_filling)
+        info() << "EliminateMatrix info=" << static_cast<int>(elimination_info) << " row="
+               << std::setw(4) << dof.localId() << " value=" << elimination_value;
+      _setMatrixValue(dof, dof, 1.0);
+    }
+  }
+}
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
+void AlephDoFLinearSystemImpl::
+_applyRHSTransformationAndFillAlephRHS()
+{
+  const bool do_print_filling = m_do_print_filling;
+
   // Apply Row+Column elimination
   // Phase 1:
   // - subtract values of the RHS vector if Row+Column elimination
-  _applyRowColumnEliminationToRHS(m_do_print_filling);
+  _applyRowColumnEliminationToRHS(do_print_filling);
+
+  IItemFamily* dof_family = dofFamily();
+
+  auto& dof_elimination_info = getEliminationInfo();
+  auto& dof_elimination_value = getEliminationValue();
+  auto& rhs_variable = rhsVariable();
 
-  // Apply Row or Row+Column elimination
-  // Phase 2: set the value of the RHS
-  // Phase 2: fill the diagonal with 1.0
+  // Apply Row or Row+Column elimination on RHS
   ENUMERATE_ (DoF, idof, dof_family->allItems()) {
     DoF dof = *idof;
     if (!dof.isOwn())
@@ -338,11 +367,13 @@ _fillMatrix()
     if (elimination_info == ELIMINATE_ROW || elimination_info == ELIMINATE_ROW_COLUMN) {
       Real elimination_value = dof_elimination_value[dof];
       rhs_variable[dof] = elimination_value;
-      info() << "Eliminate info=" << static_cast<int>(elimination_info) << " row="
-             << std::setw(4) << dof.localId() << " value=" << elimination_value;
-      _setMatrixValue(dof, dof, 1.0);
+      if (do_print_filling)
+        info() << "EliminateRHS info=" << static_cast<int>(elimination_info) << " row="
+               << std::setw(4) << dof.localId() << " value=" << elimination_value;
     }
   }
+
+  _fillRHSVector();
 }
 
 /*---------------------------------------------------------------------------*/
@@ -468,15 +499,16 @@ solve()
 {
   UniqueArray<Real> aleph_result;
 
-  // _fillMatrix() may change the values of RHS vector
-  // with row or row-column elimination so we have to fill the RHS vector
-  // before the matrix.
-  _fillMatrix();
-  _fillRHSVector();
-
   info() << "[AlephFem] Assemble matrix ptr=" << m_aleph_matrix;
+
+  // Matrix transformation
+  _applyMatrixTransformationAndFillAlephMatrix();
   m_aleph_matrix->assemble();
+
+  // RHS Transformation
+  _applyRHSTransformationAndFillAlephRHS();
   m_aleph_rhs_vector->assemble();
+
   auto* aleph_solution_vector = m_aleph_solution_vector;
   IItemFamily* dof_family = dofFamily();
   DoFGroup own_dofs = dof_family->allItems().own();

femutils/CsrDoFLinearSystemImpl.cc

@@ -76,18 +76,18 @@ _applyForcedValuesToLhs()
 /*---------------------------------------------------------------------------*/
 
 void CsrDoFLinearSystemImpl::
-_applyRowOrRowColumnElimination()
+_applyRowOrRowColumnEliminationOnMatrix()
 {
-  _applyRowElimination();
+  _applyRowEliminationOnMatrix();
   if (m_has_row_column_elimination)
-    _applyRowColumnElimination();
+    _applyRowColumnEliminationOnMatrix();
 }
 
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 
 void CsrDoFLinearSystemImpl::
-_applyRowColumnElimination()
+_applyRowColumnEliminationOnMatrix()
 {
   IItemFamily* dof_family = dofFamily();
 
@@ -118,18 +118,14 @@ _applyRowColumnElimination()
         }
       }
     }
-    if (is_row_elimination) {
-      auto elimination_value = in_elimination_value[dof_row];
-      in_out_rhs_variable[dof_row] = elimination_value;
-    }
   };
 }
 
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 
 void CsrDoFLinearSystemImpl::
-_applyRowElimination()
+_applyRowEliminationOnMatrix()
 {
   IItemFamily* dof_family = dofFamily();
 
@@ -152,7 +148,37 @@ _applyRowElimination()
       auto elimination_value = in_elimination_value[dof_id];
       for (CsrRowColumnIndex csr_index : csr_view.rowRange(dof_id))
         csr_view.value(csr_index) = (csr_view.column(csr_index) == dof_id) ? 1.0 : 0.0;
-      in_out_rhs_variable[dof_id] = elimination_value;
+    }
+  };
+}
+
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
+void CsrDoFLinearSystemImpl::
+_applyRowOrRowColumnEliminationOnRHS()
+{
+  IItemFamily* dof_family = dofFamily();
+
+  auto nb_dof = dof_family->nbItem();
+
+  RunQueue queue = makeQueue(runner());
+  auto command = makeCommand(queue);
+
+  auto in_elimination_info = Accelerator::viewIn(command, getEliminationInfo());
+  auto in_elimination_value = Accelerator::viewIn(command, getEliminationValue());
+
+  auto in_out_rhs_variable = Accelerator::viewInOut(command, rhsVariable());
+  auto csr_view = m_csr_view;
+  command << RUNCOMMAND_LOOP1(iter, nb_dof)
+  {
+    auto [row_index] = iter();
+    DoFLocalId dof_row(row_index);
+    auto row_elimination_info = in_elimination_info[dof_row];
+    bool is_row_elimination = (row_elimination_info == ELIMINATE_ROW) || (row_elimination_info == ELIMINATE_ROW_COLUMN);
+    if (is_row_elimination) {
+      auto elimination_value = in_elimination_value[dof_row];
+      in_out_rhs_variable[dof_row] = elimination_value;
     }
   };
 }

femutils/HypreDoFLinearSystem.cc

@@ -258,11 +258,17 @@ namespace
 void HypreDoFLinearSystemImpl::
 solve()
 {
+  const bool do_debug_print = false;
+
+  // Matrix transformation
   _fillRowColumnEliminationInfos();
-  _applyRowColumnEliminationToRHS(true);
-  _applyRowOrRowColumnElimination();
+  _applyRowOrRowColumnEliminationOnMatrix();
   _applyForcedValuesToLhs();
 
+  // RHS transformation
+  _applyRowColumnEliminationToRHS(do_debug_print);
+  _applyRowOrRowColumnEliminationOnRHS();
+
 #if HYPRE_RELEASE_NUMBER >= 22700
   HYPRE_MemoryLocation hypre_memory = HYPRE_MEMORY_HOST;
   HYPRE_ExecutionPolicy hypre_exec_policy = HYPRE_EXEC_HOST;
@@ -343,7 +349,6 @@ solve()
   HYPRE_IJMatrix ij_A = nullptr;
   HYPRE_ParCSRMatrix parcsr_A = nullptr;
 
-  const bool do_debug_print = false;
   const bool do_dump_matrix = false;
 
   Span<const Int32> rows_index_span = m_dof_matrix_numbering.asArray();

femutils/internal/CsrDoFLinearSystemImpl.h

@@ -100,12 +100,13 @@ class CsrDoFLinearSystemImpl
   // These methods should be private but has to be public because of NVidia compiler
   void _applyForcedValuesToLhs();
   void _fillRowColumnEliminationInfos();
-  void _applyRowElimination();
-  void _applyRowColumnElimination();
+  void _applyRowEliminationOnMatrix();
+  void _applyRowColumnEliminationOnMatrix();
+  void _applyRowOrRowColumnEliminationOnRHS();
 
  protected:
 
-  void _applyRowOrRowColumnElimination();
+  void _applyRowOrRowColumnEliminationOnMatrix();
 };
 
 /*---------------------------------------------------------------------------*/