Quantization of the Jacobian matrix to speed up linear algebra operations

.github/workflows/regression.yml

@@ -233,7 +233,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
@@ -282,7 +282,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests"
+          args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} -a "--tapetests"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
@@ -330,7 +330,7 @@ jobs:
           PMIX_MCA_gds: hash
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan"
+          args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} -a "--tsan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054
         with:
@@ -375,7 +375,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan"
+          args: -b ${{github.ref}} -t develop -c pedro/quantization -s ${{matrix.testscript}} -a "--asan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054
         with:

Common/include/basic_types/datatype_structure.hpp

@@ -135,6 +135,10 @@ template <>
 struct Passive<su2double> {
   FORCEINLINE static passivedouble Value(const su2double& val) { return GetValue(val); }
 };
+template <class T>
+FORCEINLINE auto PassiveValue(const T& val) {
+  return Passive<T>::Value(val);
+}
 
 /*!
  * \brief Casts the primitive value to int (uses GetValue, already implemented for each type).

Common/include/code_config.hpp

@@ -26,6 +26,7 @@
  */
 #pragma once
 
+#include <cstdint>
 #include <type_traits>
 #include <cmath>
 
@@ -100,6 +101,10 @@ FORCEINLINE Out su2staticcast_p(In ptr) {
 #undef USE_SINGLE_PRECISION
 #endif
 
+/*--- Default integer types. Currently used for rank-local sparse patterns. ---*/
+using su2uint = uint32_t;
+using su2int = int32_t;
+
 /*--- This type can be used for (rare) compatibility cases or for
  * computations that are intended to be (always) passive. ---*/
 #ifdef USE_SINGLE_PRECISION

Common/include/geometry/CGeometry.hpp

@@ -205,10 +205,10 @@ class CGeometry {
   LDUSparsePattern finiteElementPatternFill0; /*!< \brief FEM sparsity with 0-fill (structural pattern). */
   LDUSparsePattern finiteElementPatternFillN; /*!< \brief FEM sparsity with N-fill (e.g. for ILU-N). */
 
-  su2vector<unsigned long> finiteVolumeLToUTranspMap;  /*!< \brief FVM L-entry -> U-entry of its transpose. */
-  su2vector<unsigned long> finiteVolumeUToLTranspMap;  /*!< \brief FVM U-entry -> L-entry of its transpose. */
-  su2vector<unsigned long> finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */
-  su2vector<unsigned long> finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */
+  su2vector<su2uint> finiteVolumeLToUTranspMap;  /*!< \brief FVM L-entry -> U-entry of its transpose. */
+  su2vector<su2uint> finiteVolumeUToLTranspMap;  /*!< \brief FVM U-entry -> L-entry of its transpose. */
+  su2vector<su2uint> finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */
+  su2vector<su2uint> finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */
 
   /*--- Edge and element colorings. ---*/
 
@@ -1892,15 +1892,15 @@ class CGeometry {
    * \param[in] type - Finite volume or finite element.
    * \return Reference to the l_to_u map.
    */
-  const su2vector<unsigned long>& GetLToUTransposeSparsePatternMap(ConnectivityType type);
+  const su2vector<su2uint>& GetLToUTransposeSparsePatternMap(ConnectivityType type);
 
   /*!
    * \brief Get the bijective map from U-entry indices to L-entry indices of their transposes.
    * \note Requires symmetric pattern. Builds both LU transpose maps if not already built.
    * \param[in] type - Finite volume or finite element.
    * \return Reference to the u_to_l map.
    */
-  const su2vector<unsigned long>& GetUToLTransposeSparsePatternMap(ConnectivityType type);
+  const su2vector<su2uint>& GetUToLTransposeSparsePatternMap(ConnectivityType type);
 
   /*!
    * \brief Get the edge coloring.

Common/include/linear_algebra/CPastixWrapper.hpp

@@ -28,6 +28,8 @@
 
 #pragma once
 
+#include "../code_config.hpp"
+
 #ifdef HAVE_PASTIX
 
 #ifdef CODI_FORWARD_TYPE
@@ -74,11 +76,11 @@ class CPastixWrapper {
     unsigned long nPointDomain = 0;
     unsigned long blkSz = 0; /*!< \brief Block size (nVar * nVar) for value assembly. */
 
-    const unsigned long* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */
-    const unsigned long* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */
-    const ScalarType* d = nullptr;            /*!< \brief Diagonal blocks (matrix-owned). */
-    const ScalarType* l = nullptr;            /*!< \brief Lower blocks (matrix-owned). */
-    const ScalarType* u = nullptr;            /*!< \brief Upper blocks (matrix-owned). */
+    const su2uint* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */
+    const su2uint* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */
+    const ScalarType* d = nullptr;      /*!< \brief Diagonal blocks (matrix-owned). */
+    const ScalarType* l = nullptr;      /*!< \brief Lower blocks (matrix-owned). */
+    const ScalarType* u = nullptr;      /*!< \brief Upper blocks (matrix-owned). */
 
     unsigned long size_rhs() const { return nPointDomain * nVar; }
   } matrix; /*!< \brief Dimensions and LDU pointers captured from the owning CSysMatrix. */
@@ -150,9 +152,9 @@ class CPastixWrapper {
    * \param[in] col_ind_l/u - LDU lower/upper column indices (geometry-owned).
    * \param[in] d/l/u - LDU value blocks (matrix-owned, must outlive wrapper).
    */
-  void SetLDU(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const unsigned long* row_ptr_l,
-              const unsigned long* col_ind_l, const unsigned long* row_ptr_u, const unsigned long* col_ind_u,
-              const ScalarType* d, const ScalarType* l, const ScalarType* u) {
+  void SetLDU(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const su2uint* row_ptr_l,
+              const su2uint* col_ind_l, const su2uint* row_ptr_u, const su2uint* col_ind_u, const ScalarType* d,
+              const ScalarType* l, const ScalarType* u) {
     if (issetup) return;
     matrix.nVar = nVar;
     matrix.nPoint = nPoint;

Common/include/linear_algebra/CPreconditioner.hpp

@@ -210,6 +210,37 @@ class CLU_SGSPreconditioner final : public CPreconditioner<ScalarType> {
   }
 };
 
+/*!
+ * \class CQuantizedLUSGSPreconditioner
+ * \brief Specialization of preconditioner that uses CSysMatrix class.
+ */
+template <class ScalarType>
+class CQuantizedLUSGSPreconditioner final : public CPreconditioner<ScalarType> {
+ private:
+  CSysMatrix<ScalarType>& sparse_matrix;
+  CGeometry* geometry;
+  const CConfig* config;
+
+ public:
+  inline CQuantizedLUSGSPreconditioner(CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref,
+                                       const CConfig* config_ref)
+      : sparse_matrix(matrix_ref) {
+    if ((geometry_ref == nullptr) || (config_ref == nullptr))
+      SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
+    geometry = geometry_ref;
+    config = config_ref;
+  }
+
+  CQuantizedLUSGSPreconditioner() = delete;
+
+  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
+    sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config);
+  }
+
+  /*! \brief Quantize the diagonal blocks (off diagonals are quantized on the fly). */
+  inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); }
+};
+
 /*!
  * \class CLineletPreconditioner
  * \brief Specialization of preconditioner that uses CSysMatrix class.
@@ -341,6 +372,9 @@ CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOL
     case LU_SGS:
       prec = new CLU_SGSPreconditioner<ScalarType>(jacobian, geometry, config);
       break;
+    case Q_LU_SGS:
+      prec = new CQuantizedLUSGSPreconditioner<ScalarType>(jacobian, geometry, config);
+      break;
     case ILU:
       prec = new CILUPreconditioner<ScalarType>(jacobian, geometry, config);
       break;