LA: add a fast path for updating matrix values

The fast path will use the PETSc function MatSetValuesRow to update
all the values of a row at once (without the need to get the row
pattern).

Fast update can be performed if:
 - the matrix has already been assembled;
 - the assembler is providing all the values of the row;
 - values provided by the assembler are sorted by ascending column.

src/LA/system_solvers_large.cpp

@@ -58,6 +58,20 @@ SystemSparseMatrixAssembler::SystemSparseMatrixAssembler(const SparseMatrix *mat
 {
 }
 
+/*!
+ * Get the assembly options.
+ *
+ * \result The assembly options that will be used.
+ */
+SystemMatrixAssembler::AssemblyOptions SystemSparseMatrixAssembler::getOptions() const
+{
+    AssemblyOptions options;
+    options.full   = true;
+    options.sorted = false;
+
+    return options;
+}
+
 /*!
  * Get the number of rows of the matrix.
  *
@@ -1094,6 +1108,25 @@ void SystemSolver::matrixUpdate(long nRows, const long *rows, const SystemMatrix
     PetscInt rowGlobalOffset;
     MatGetOwnershipRangeColumn(m_A, &rowGlobalOffset, nullptr);
 
+    // Get the options for assembling the matrix
+    SystemMatrixAssembler::AssemblyOptions assemblyOptions = assembler.getOptions();
+
+    // Check if the assembler provides full and sorted rows
+    //
+    // The option MAT_SORTED_FULL means that each process provides exactly its
+    // local rows; all column indices for a given row are passed in a single call
+    // to MatSetValues(), preallocation is perfect, row oriented, INSERT_VALUES
+    // is used.
+    //
+    // This options needs at least PETSc 3.12.
+    PetscBool matrixSortedFull;
+#if PETSC_VERSION_GE(3, 12, 0)
+    matrixSortedFull = (assemblyOptions.full && assemblyOptions.sorted) ? PETSC_TRUE : PETSC_FALSE;
+    MatSetOption(m_A, MAT_SORTED_FULL, matrixSortedFull);
+#else
+    matrixSortedFull = PETSC_FALSE;
+#endif
+
     // Update element values
     //
     // If the sizes of PETSc data types match the sizes of data types expected by
@@ -1127,11 +1160,20 @@ void SystemSolver::matrixUpdate(long nRows, const long *rows, const SystemMatrix
 
         const PetscInt globalRow = rowGlobalOffset + row;
 
+        // Check if it possible to perform a fast update
+        //
+        // A fast update allows to set all the values of a row at once (without
+        // the need to get the row pattern), it can be performed if:
+        //  - the matrix has already been assembled;
+        //  - the assembler is providing all the values of the row;
+        //  - values provided by the assembler are sorted by ascending column.
+        bool fastUpdate = isAssembled() && (matrixSortedFull == PETSC_TRUE);
+
         // Get row data
-        assembler.getRowData(assemblerRow, &rowPattern, &rowValues);
-        const int nRowValues = rowValues.size();
-        if (nRowValues == 0) {
-            continue;
+        if (fastUpdate) {
+            assembler.getRowValues(assemblerRow, &rowValues);
+        } else {
+            assembler.getRowData(assemblerRow, &rowPattern, &rowValues);
         }
 
         // Get values in PETSc format
@@ -1141,22 +1183,33 @@ void SystemSolver::matrixUpdate(long nRows, const long *rows, const SystemMatrix
             std::copy(rowValues.cbegin(), rowValues.cend(), petscRowValuesStorage.begin());
         }
 
-        // Get pattern in PETSc format
-        for (int k = 0; k < nRowValues; ++k) {
-            long globalCol = rowPattern[k];
-            if (m_colPermutation) {
-                if (globalCol >= colGlobalBegin && globalCol < colGlobalEnd) {
-                    long col = globalCol - colGlobalBegin;
-                    col = colInvRanks[col];
-                    globalCol = colGlobalBegin + col;
+        if (fastUpdate) {
+            // Set values
+            MatSetValuesRow(m_A, globalRow, petscRowValues);
+        } else {
+            // Count the values that will be updated
+            const int nRowValues = rowValues.size();
+            if (nRowValues == 0) {
+                continue;
+            }
+
+            // Get pattern in PETSc format
+            for (int k = 0; k < nRowValues; ++k) {
+                long globalCol = rowPattern[k];
+                if (m_colPermutation) {
+                    if (globalCol >= colGlobalBegin && globalCol < colGlobalEnd) {
+                        long col = globalCol - colGlobalBegin;
+                        col = colInvRanks[col];
+                        globalCol = colGlobalBegin + col;
+                    }
                 }
+
+                petscRowPattern[k] = globalCol;
             }
 
-            petscRowPattern[k] = globalCol;
+            // Set data
+            MatSetValues(m_A, 1, &globalRow, nRowValues, petscRowPattern.data(), petscRowValues, INSERT_VALUES);
         }
-
-        // Set data
-        MatSetValues(m_A, 1, &globalRow, nRowValues, petscRowPattern.data(), petscRowValues, INSERT_VALUES);
     }
 
     // Let petsc assembly the matrix after the update

src/LA/system_solvers_large.hpp

@@ -70,8 +70,15 @@ struct KSPStatus {
 class SystemMatrixAssembler {
 
 public:
+    struct AssemblyOptions {
+        bool full; //!< Controls if the assembler is providing all the non-zero values of a row
+        bool sorted; //! Controls if the values provided by the assembler are sorted by ascending column
+    };
+
     virtual ~SystemMatrixAssembler() = default;
 
+    virtual AssemblyOptions getOptions() const = 0;
+
     virtual long getRowCount() const = 0;
     virtual long getColCount() const = 0;
 
@@ -100,6 +107,8 @@ class SystemSparseMatrixAssembler : public SystemMatrixAssembler {
 public:
     SystemSparseMatrixAssembler(const SparseMatrix *matrix);
 
+    AssemblyOptions getOptions() const override;
+
     long getRowCount() const override;
     long getColCount() const override;
 

src/discretization/stencil_solver.hpp

@@ -122,6 +122,8 @@ class DiscretizationStencilSolverAssembler : public StencilSolverAssembler {
     DiscretizationStencilSolverAssembler(MPI_Comm communicator, bool partitioned, const stencil_container_t *stencils);
 #endif
 
+    AssemblyOptions getOptions() const override;
+
     int getBlockSize() const;
 
     long getRowCount() const override;

src/discretization/stencil_solver.tpp

@@ -298,6 +298,21 @@ DiscretizationStencilSolverAssembler<stencil_t>::DiscretizationStencilSolverAsse
 {
 }
 
+/*!
+ * Get the assembly options.
+ *
+ * \result The assembly options that will be used.
+ */
+template<typename stencil_t>
+SystemMatrixAssembler::AssemblyOptions DiscretizationStencilSolverAssembler<stencil_t>::getOptions() const
+{
+    AssemblyOptions options;
+    options.full   = true;
+    options.sorted = false;
+
+    return options;
+}
+
 /*!
  * Set block size.
  *