added symmetric matrix type for the equilibrium enforcement

Author: rrsettgast

src/common/DirectSystemSolve.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "macros.hpp"
+#include "symmetricMatrix.hpp"
 #include <cmath>
 
 namespace hpcReact
@@ -75,6 +76,47 @@ void solveNxN_Cholesky(REAL_TYPE const (&A)[N][N], REAL_TYPE const (&b)[N], REAL
 }
 
 
+
+template< typename REAL_TYPE, int N >
+void solveNxN_Cholesky( symmetricMatrix<REAL_TYPE,int,N> const & A, 
+                        REAL_TYPE const (&b)[N], 
+                        REAL_TYPE (&x)[N]) 
+{
+  symmetricMatrix<REAL_TYPE,int,N> L = {0};
+
+  // **Cholesky Decomposition**
+  for (int i = 0; i < N; i++) {
+      for (int j = 0; j <= i; j++) {
+          REAL_TYPE sum = 0;
+          for (int k = 0; k < j; k++) {
+              sum += L(i,k) * L(j,k);
+          }
+          if (i == j)
+              L(i,j) = sqrt(A(i,i) - sum);
+          else
+              L(i,j) = (A(i,j) - sum) / L[j][j];
+      }
+  }
+
+  // **Forward Substitution: Solve L y = b**
+  REAL_TYPE y[N];
+  for (int i = 0; i < N; i++) {
+      y[i] = b[i];
+      for (int j = 0; j < i; j++)
+          y[i] -= L(i,j) * y[j];
+      y[i] /= L(i,i);
+  }
+
+  // **Backward Substitution: Solve L^T x = y**
+  for (int i = N - 1; i >= 0; i--) {
+      x[i] = y[i];
+      for (int j = i + 1; j < N; j++)
+          x[i] -= L(j,i) * x[j];
+      x[i] /= L(i,i);
+  }
+}
+
+
 template< typename REAL_TYPE, int N >
 HPCREACT_HOST_DEVICE
 void solveNxN_pivoted( REAL_TYPE (&A)[N][N], REAL_TYPE (&b)[N], REAL_TYPE (&x)[N] )