@@ -8,118 +8,130 @@ namespace hpcReact
88{
99
1010template < typename REAL_TYPE, int N >
11- bool isPositiveDefinite ( REAL_TYPE const (&A)[N][N])
11+ bool isPositiveDefinite ( REAL_TYPE const (&A)[N][N] )
1212{
1313 REAL_TYPE temp[N][N];
1414
1515 // Copy A into temp to avoid modifying original
16- for ( int i = 0 ; i < N; i++)
16+ for ( int i = 0 ; i < N; i++ )
1717 {
18- for ( int j = 0 ; j < N; j++)
18+ for ( int j = 0 ; j < N; j++ )
1919 {
2020 temp[i][j] = A[i][j];
2121 }
2222 }
2323
24- for ( int k = 0 ; k < N; k++)
24+ for ( int k = 0 ; k < N; k++ )
2525 {
2626 // Compute determinant of k-th leading principal minor using Gaussian elimination
27- if (temp[k][k] <= 0 ) return false ; // Must be positive
27+ if ( temp[k][k] <= 0 )
28+ return false ; // Must be positive
2829
29- for ( int i = k + 1 ; i < N; i++)
30+ for ( int i = k + 1 ; i < N; i++ )
3031 {
3132 REAL_TYPE factor = temp[i][k] / temp[k][k];
32- for ( int j = k; j < N; j++)
33+ for ( int j = k; j < N; j++ )
3334 {
3435 temp[i][j] -= factor * temp[k][j];
3536 }
3637 }
3738 }
38- return true ;
39+ return true ;
3940}
4041
4142
4243template < typename REAL_TYPE, int N >
43- void solveNxN_Cholesky (REAL_TYPE const (&A)[N][N], REAL_TYPE const (&b)[N], REAL_TYPE (&x)[N]) {
44+ void solveNxN_Cholesky ( REAL_TYPE const (&A)[N][N], REAL_TYPE const (&b)[N], REAL_TYPE (& x)[N] )
45+ {
4446 REAL_TYPE L[N][N] = {{0 }};
4547
4648 // **Cholesky Decomposition**
47- for (int i = 0 ; i < N; i++) {
48- for (int j = 0 ; j <= i; j++) {
49- REAL_TYPE sum = 0 ;
50- for (int k = 0 ; k < j; k++) {
51- sum += L[i][k] * L[j][k];
52- }
53- if (i == j)
54- L[i][j] = sqrt (A[i][i] - sum);
55- else
56- L[i][j] = (A[i][j] - sum) / L[j][j];
49+ for ( int i = 0 ; i < N; i++ )
50+ {
51+ for ( int j = 0 ; j <= i; j++ )
52+ {
53+ REAL_TYPE sum = 0 ;
54+ for ( int k = 0 ; k < j; k++ )
55+ {
56+ sum += L[i][k] * L[j][k];
5757 }
58+ if ( i == j )
59+ L[i][j] = sqrt ( A[i][i] - sum );
60+ else
61+ L[i][j] = (A[i][j] - sum) / L[j][j];
62+ }
5863 }
5964
6065 // **Forward Substitution: Solve L y = b**
6166 REAL_TYPE y[N];
62- for (int i = 0 ; i < N; i++) {
63- y[i] = b[i];
64- for (int j = 0 ; j < i; j++)
65- y[i] -= L[i][j] * y[j];
66- y[i] /= L[i][i];
67+ for ( int i = 0 ; i < N; i++ )
68+ {
69+ y[i] = b[i];
70+ for ( int j = 0 ; j < i; j++ )
71+ y[i] -= L[i][j] * y[j];
72+ y[i] /= L[i][i];
6773 }
6874
6975 // **Backward Substitution: Solve L^T x = y**
70- for (int i = N - 1 ; i >= 0 ; i--) {
71- x[i] = y[i];
72- for (int j = i + 1 ; j < N; j++)
73- x[i] -= L[j][i] * x[j];
74- x[i] /= L[i][i];
76+ for ( int i = N - 1 ; i >= 0 ; i-- )
77+ {
78+ x[i] = y[i];
79+ for ( int j = i + 1 ; j < N; j++ )
80+ x[i] -= L[j][i] * x[j];
81+ x[i] /= L[i][i];
7582 }
7683}
7784
7885
7986
8087template < typename REAL_TYPE, int N >
81- void solveNxN_Cholesky ( symmetricMatrix<REAL_TYPE,int ,N > const & A,
82- REAL_TYPE const (&b)[N],
83- REAL_TYPE (&x)[N])
88+ void solveNxN_Cholesky ( symmetricMatrix< REAL_TYPE, int , N > const & A,
89+ REAL_TYPE const (&b)[N],
90+ REAL_TYPE (& x)[N] )
8491{
85- symmetricMatrix<REAL_TYPE,int ,N > L = {0 };
92+ symmetricMatrix< REAL_TYPE, int , N > L = {0 };
8693
8794 // **Cholesky Decomposition**
88- for (int i = 0 ; i < N; i++) {
89- for (int j = 0 ; j <= i; j++) {
90- REAL_TYPE sum = 0 ;
91- for (int k = 0 ; k < j; k++) {
92- sum += L (i,k) * L (j,k);
93- }
94- if (i == j)
95- L (i,j) = sqrt (A (i,i) - sum);
96- else
97- L (i,j) = (A (i,j) - sum) / L[j][j];
95+ for ( int i = 0 ; i < N; i++ )
96+ {
97+ for ( int j = 0 ; j <= i; j++ )
98+ {
99+ REAL_TYPE sum = 0 ;
100+ for ( int k = 0 ; k < j; k++ )
101+ {
102+ sum += L ( i, k ) * L ( j, k );
98103 }
104+ if ( i == j )
105+ L ( i, j ) = sqrt ( A ( i, i ) - sum );
106+ else
107+ L ( i, j ) = (A ( i, j ) - sum) / L[j][j];
108+ }
99109 }
100110
101111 // **Forward Substitution: Solve L y = b**
102112 REAL_TYPE y[N];
103- for (int i = 0 ; i < N; i++) {
104- y[i] = b[i];
105- for (int j = 0 ; j < i; j++)
106- y[i] -= L (i,j) * y[j];
107- y[i] /= L (i,i);
113+ for ( int i = 0 ; i < N; i++ )
114+ {
115+ y[i] = b[i];
116+ for ( int j = 0 ; j < i; j++ )
117+ y[i] -= L ( i, j ) * y[j];
118+ y[i] /= L ( i, i );
108119 }
109120
110121 // **Backward Substitution: Solve L^T x = y**
111- for (int i = N - 1 ; i >= 0 ; i--) {
112- x[i] = y[i];
113- for (int j = i + 1 ; j < N; j++)
114- x[i] -= L (j,i) * x[j];
115- x[i] /= L (i,i);
122+ for ( int i = N - 1 ; i >= 0 ; i-- )
123+ {
124+ x[i] = y[i];
125+ for ( int j = i + 1 ; j < N; j++ )
126+ x[i] -= L ( j, i ) * x[j];
127+ x[i] /= L ( i, i );
116128 }
117129}
118130
119131
120132template < typename REAL_TYPE, int N >
121133HPCREACT_HOST_DEVICE
122- void solveNxN_pivoted ( REAL_TYPE (&A)[N][N], REAL_TYPE (&b)[N], REAL_TYPE (&x)[N] )
134+ void solveNxN_pivoted ( REAL_TYPE (& A)[N][N], REAL_TYPE (& b)[N], REAL_TYPE (& x)[N] )
123135{
124136
125137
0 commit comments