1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Eigen system solver for symmetric NxN matrices. Modified from the CC0 code provided at https://github.com/jewettaij/jacobi_pd
10 #pragma once
11
12 #include <ceed/types.h>
13 #ifndef CEED_RUNNING_JIT_PASS
14 #include <math.h>
15 #include <stdbool.h>
16 #endif
17
18 #include "utils.h"
19
20 // @typedef choose the criteria for sorting eigenvalues and eigenvectors
21 typedef enum eSortCriteria {
22 SORT_NONE,
23 SORT_DECREASING_EVALS,
24 SORT_INCREASING_EVALS,
25 SORT_DECREASING_ABS_EVALS,
26 SORT_INCREASING_ABS_EVALS
27 } SortCriteria;
28
29 ///@brief Find the off-diagonal index in row i whose absolute value is largest
30 ///
31 /// @param[in] *A matrix
32 /// @param[in] i row index
33 /// @returns Index of absolute largest off-diagonal element in row i
MaxEntryRow(const CeedScalar * A,CeedInt N,CeedInt i)34 CEED_QFUNCTION_HELPER CeedInt MaxEntryRow(const CeedScalar *A, CeedInt N, CeedInt i) {
35 CeedInt j_max = i + 1;
36 for (CeedInt j = i + 2; j < N; j++)
37 if (fabs(A[i * N + j]) > fabs(A[i * N + j_max])) j_max = j;
38 return j_max;
39 }
40
41 /// @brief Find the indices (i_max, j_max) marking the location of the
42 /// entry in the matrix with the largest absolute value. This
43 /// uses the max_idx_row[] array to find the answer in O(n) time.
44 ///
45 /// @param[in] *A matrix
46 /// @param[inout] i_max row index
47 /// @param[inout] j_max column index
MaxEntry(const CeedScalar * A,CeedInt N,CeedInt * max_idx_row,CeedInt * i_max,CeedInt * j_max)48 CEED_QFUNCTION_HELPER void MaxEntry(const CeedScalar *A, CeedInt N, CeedInt *max_idx_row, CeedInt *i_max, CeedInt *j_max) {
49 *i_max = 0;
50 *j_max = max_idx_row[*i_max];
51 CeedScalar max_entry = fabs(A[*i_max * N + *j_max]);
52 for (CeedInt i = 1; i < N - 1; i++) {
53 CeedInt j = max_idx_row[i];
54 if (fabs(A[i * N + j]) > max_entry) {
55 max_entry = fabs(A[i * N + j]);
56 *i_max = i;
57 *j_max = j;
58 }
59 }
60 }
61
62 /// @brief Calculate the components of a rotation matrix which performs a
63 /// rotation in the i,j plane by an angle (θ) that (when multiplied on
64 /// both sides) will zero the ij'th element of A, so that afterwards
65 /// A[i][j] = 0. The results will be stored in c, s, and t
66 /// (which store cos(θ), sin(θ), and tan(θ), respectively).
67 ///
68 /// @param[in] *A matrix
69 /// @param[in] i row index
70 /// @param[in] j column index
CalcRot(const CeedScalar * A,CeedInt N,CeedInt i,CeedInt j,CeedScalar * rotmat_cst)71 CEED_QFUNCTION_HELPER void CalcRot(const CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedScalar *rotmat_cst) {
72 rotmat_cst[2] = 1.0; // = tan(θ)
73 CeedScalar A_jj_ii = (A[j * N + j] - A[i * N + i]);
74 if (A_jj_ii != 0.0) {
75 // kappa = (A[j][j] - A[i][i]) / (2*A[i][j])
76 CeedScalar kappa = A_jj_ii;
77 rotmat_cst[2] = 0.0;
78 CeedScalar A_ij = A[i * N + j];
79 if (A_ij != 0.0) {
80 kappa /= (2.0 * A_ij);
81 // t satisfies: t^2 + 2*t*kappa - 1 = 0
82 // (choose the root which has the smaller absolute value)
83 rotmat_cst[2] = 1.0 / (sqrt(1 + kappa * kappa) + fabs(kappa));
84 if (kappa < 0.0) rotmat_cst[2] = -rotmat_cst[2];
85 }
86 }
87 rotmat_cst[0] = 1.0 / sqrt(1 + rotmat_cst[2] * rotmat_cst[2]);
88 rotmat_cst[1] = rotmat_cst[0] * rotmat_cst[2];
89 }
90
91 /// @brief Perform a similarity transformation by multiplying matrix A on both
92 /// sides by a rotation matrix (and its transpose) to eliminate A[i][j].
93 /// @details This rotation matrix performs a rotation in the i,j plane by
94 /// angle θ. This function assumes that c=cos(θ). s=sin(θ), t=tan(θ)
95 /// have been calculated in advance (using the CalcRot() function).
96 /// It also assumes that i<j. The max_idx_row[] array is also updated.
97 /// To save time, since the matrix is symmetric, the elements
98 /// below the diagonal (ie. A[u][v] where u>v) are not computed.
99 /// @verbatim
100 /// A' = R^T * A * R
101 /// where R the rotation in the i,j plane and ^T denotes the transpose.
102 /// i j
103 /// _ _
104 /// | 1 |
105 /// | . |
106 /// | . |
107 /// | 1 |
108 /// | c ... s |
109 /// | . . . |
110 /// R = | . 1 . |
111 /// | . . . |
112 /// | -s ... c |
113 /// | 1 |
114 /// | . |
115 /// | . |
116 /// |_ 1 _|
117 /// @endverbatim
118 ///
119 /// Let A' denote the matrix A after multiplication by R^T and R.
120 /// The components of A' are:
121 ///
122 /// @verbatim
123 /// A'_uv = Σ_w Σ_z R_wu * A_wz * R_zv
124 /// @endverbatim
125 ///
126 /// Note that a the rotation at location i,j will modify all of the matrix
127 /// elements containing at least one index which is either i or j
128 /// such as: A[w][i], A[i][w], A[w][j], A[j][w].
129 /// Check and see whether these modified matrix elements exceed the
130 /// corresponding values in max_idx_row[] array for that row.
131 /// If so, then update max_idx_row for that row.
132 /// This is somewhat complicated by the fact that we must only consider
133 /// matrix elements in the upper-right triangle strictly above the diagonal.
134 /// (ie. matrix elements whose second index is > the first index).
135 /// The modified elements we must consider are marked with an "X" below:
136 ///
137 /// @verbatim
138 /// i j
139 /// _ _
140 /// | . X X |
141 /// | . X X |
142 /// | . X X |
143 /// | . X X |
144 /// | X X X X X 0 X X X X | i
145 /// | . X |
146 /// | . X |
147 /// A = | . X |
148 /// | . X |
149 /// | X X X X X | j
150 /// | . |
151 /// | . |
152 /// | . |
153 /// |_ . _|
154 /// @endverbatim
155 ///
156 /// @param[in] *A matrix
157 /// @param[in] i row index
158 /// @param[in] j column index
ApplyRot(CeedScalar * A,CeedInt N,CeedInt i,CeedInt j,CeedInt * max_idx_row,CeedScalar * rotmat_cst)159 CEED_QFUNCTION_HELPER void ApplyRot(CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedInt *max_idx_row, CeedScalar *rotmat_cst) {
160 // Compute the diagonal elements of A which have changed:
161 A[i * N + i] -= rotmat_cst[2] * A[i * N + j];
162 A[j * N + j] += rotmat_cst[2] * A[i * N + j];
163 // Note: This is algebraically equivalent to:
164 // A[i][i] = c*c*A[i][i] + s*s*A[j][j] - 2*s*c*A[i][j]
165 // A[j][j] = s*s*A[i][i] + c*c*A[j][j] + 2*s*c*A[i][j]
166
167 // Update the off-diagonal elements of A which will change (above the diagonal)
168
169 A[i * N + j] = 0.0;
170
171 // compute A[w][i] and A[i][w] for all w!=i,considering above-diagonal elements
172 for (CeedInt w = 0; w < i; w++) { // 0 <= w < i < j < N
173 A[i * N + w] = A[w * N + i]; // backup the previous value. store below diagonal (i>w)
174 A[w * N + i] = rotmat_cst[0] * A[w * N + i] - rotmat_cst[1] * A[w * N + j]; // A[w][i], A[w][j] from previous iteration
175 if (i == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(A, N, w);
176 else if (fabs(A[w * N + i]) > fabs(A[w * N + max_idx_row[w]])) max_idx_row[w] = i;
177 }
178 for (CeedInt w = i + 1; w < j; w++) { // 0 <= i < w < j < N
179 A[w * N + i] = A[i * N + w]; // backup the previous value. store below diagonal (w>i)
180 A[i * N + w] = rotmat_cst[0] * A[i * N + w] - rotmat_cst[1] * A[w * N + j]; // A[i][w], A[w][j] from previous iteration
181 }
182 for (CeedInt w = j + 1; w < N; w++) { // 0 <= i < j+1 <= w < N
183 A[w * N + i] = A[i * N + w]; // backup the previous value. store below diagonal (w>i)
184 A[i * N + w] = rotmat_cst[0] * A[i * N + w] - rotmat_cst[1] * A[j * N + w]; // A[i][w], A[j][w] from previous iteration
185 }
186
187 // now that we're done modifying row i, we can update max_idx_row[i]
188 max_idx_row[i] = MaxEntryRow(A, N, i);
189
190 // compute A[w][j] and A[j][w] for all w!=j,considering above-diagonal elements
191 for (CeedInt w = 0; w < i; w++) { // 0 <= w < i < j < N
192 A[w * N + j] = rotmat_cst[1] * A[i * N + w] + rotmat_cst[0] * A[w * N + j]; // A[i][w], A[w][j] from previous iteration
193 if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(A, N, w);
194 else if (fabs(A[w * N + j]) > fabs(A[w * N + max_idx_row[w]])) max_idx_row[w] = j;
195 }
196 for (CeedInt w = i + 1; w < j; w++) { // 0 <= i+1 <= w < j < N
197 A[w * N + j] = rotmat_cst[1] * A[w * N + i] + rotmat_cst[0] * A[w * N + j]; // A[w][i], A[w][j] from previous iteration
198 if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(A, N, w);
199 else if (fabs(A[w * N + j]) > fabs(A[w * N + max_idx_row[w]])) max_idx_row[w] = j;
200 }
201 for (CeedInt w = j + 1; w < N; w++) { // 0 <= i < j < w < N
202 A[j * N + w] = rotmat_cst[1] * A[w * N + i] + rotmat_cst[0] * A[j * N + w]; // A[w][i], A[j][w] from previous iteration
203 }
204 // now that we're done modifying row j, we can update max_idx_row[j]
205 max_idx_row[j] = MaxEntryRow(A, N, j);
206 }
207
208 ///@brief Multiply matrix A on the LEFT side by a transposed rotation matrix R^T
209 /// This matrix performs a rotation in the i,j plane by angle θ (where
210 /// the arguments "s" and "c" refer to cos(θ) and sin(θ), respectively).
211 /// @verbatim
212 /// A'_uv = Σ_w R_wu * A_wv
213 /// @endverbatim
214 ///
215 /// @param[in] *A matrix
216 /// @param[in] i row index
217 /// @param[in] j column index
ApplyRotLeft(CeedScalar * A,CeedInt N,CeedInt i,CeedInt j,CeedScalar * rotmat_cst)218 CEED_QFUNCTION_HELPER void ApplyRotLeft(CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedScalar *rotmat_cst) {
219 // Recall that c = cos(θ) and s = sin(θ)
220 for (CeedInt v = 0; v < N; v++) {
221 CeedScalar Aiv = A[i * N + v];
222 A[i * N + v] = rotmat_cst[0] * A[i * N + v] - rotmat_cst[1] * A[j * N + v];
223 A[j * N + v] = rotmat_cst[1] * Aiv + rotmat_cst[0] * A[j * N + v];
224 }
225 }
226
227 /// @brief Sort the rows in evec according to the numbers in v (also sorted)
228 ///
229 /// @param[inout] *eval vector containing the keys used for sorting
230 /// @param[inout] *evec matrix whose rows will be sorted according to v
231 /// @param[in] n size of the vector and matrix
232 /// @param[in] s sort decreasing order?
SortRows(CeedScalar * eval,CeedScalar * evec,CeedInt N,SortCriteria sort_criteria)233 CEED_QFUNCTION_HELPER void SortRows(CeedScalar *eval, CeedScalar *evec, CeedInt N, SortCriteria sort_criteria) {
234 if (sort_criteria == SORT_NONE) return;
235
236 for (CeedInt i = 0; i < N - 1; i++) {
237 CeedInt i_max = i;
238 for (CeedInt j = i + 1; j < N; j++) {
239 // find the "maximum" element in the array starting at position i+1
240 switch (sort_criteria) {
241 case SORT_DECREASING_EVALS:
242 if (eval[j] > eval[i_max]) i_max = j;
243 break;
244 case SORT_INCREASING_EVALS:
245 if (eval[j] < eval[i_max]) i_max = j;
246 break;
247 case SORT_DECREASING_ABS_EVALS:
248 if (fabs(eval[j]) > fabs(eval[i_max])) i_max = j;
249 break;
250 case SORT_INCREASING_ABS_EVALS:
251 if (fabs(eval[j]) < fabs(eval[i_max])) i_max = j;
252 break;
253 default:
254 break;
255 }
256 }
257 SwapScalar(&eval[i], &eval[i_max]);
258 for (CeedInt k = 0; k < N; k++) SwapScalar(&evec[i * N + k], &evec[i_max * N + k]);
259 }
260 }
261
262 /// @brief Calculate all the eigenvalues and eigevectors of a symmetric matrix
263 /// using the Jacobi eigenvalue algorithm:
264 /// https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
265 /// @returns The number of Jacobi iterations attempted, which should be > 0.
266 /// If the return value is not strictly > 0 then convergence failed.
267 /// @note To reduce the computation time further, set calc_evecs=false.
268 /// Additionally, note that the output evecs should be normalized. It
269 /// simply takes the Identity matrix and performs (isometric) rotations
270 /// on it, so divergence from normalized is due to finite-precision
271 /// arithmetic of the rotations.
272 //
273 // @param[in] A the matrix you wish to diagonalize (size NxN)
274 // @param[in] N size of the matrix
275 // @param[out] eval store the eigenvalues here (size N)
276 // @param[out] evec store the eigenvectors here (in rows, size NxN)
277 // @param[out] max_idx_row work vector of size N
278 // @param[in] sort_criteria sort results?
279 // @param[in] calc_evecs calculate the eigenvectors?
280 // @param[in] max_num_sweeps maximum number of iterations = max_num_sweeps * number of off-diagonals (N*(N-1)/2)
Diagonalize(CeedScalar * A,CeedInt N,CeedScalar * eval,CeedScalar * evec,CeedInt * max_idx_row,SortCriteria sort_criteria,bool calc_evec,const CeedInt max_num_sweeps)281 CEED_QFUNCTION_HELPER CeedInt Diagonalize(CeedScalar *A, CeedInt N, CeedScalar *eval, CeedScalar *evec, CeedInt *max_idx_row,
282 SortCriteria sort_criteria, bool calc_evec, const CeedInt max_num_sweeps) {
283 CeedScalar rotmat_cst[3] = {0.}; // cos(θ), sin(θ), and tan(θ),
284
285 if (calc_evec)
286 for (CeedInt i = 0; i < N; i++)
287 for (CeedInt j = 0; j < N; j++) evec[i * N + j] = (i == j) ? 1.0 : 0.0; // Set evec equal to the identity matrix
288
289 for (CeedInt i = 0; i < N - 1; i++) max_idx_row[i] = MaxEntryRow(A, N, i);
290
291 // -- Iteration --
292 CeedInt n_iters;
293 CeedInt max_num_iters = max_num_sweeps * N * (N - 1) / 2;
294 for (n_iters = 1; n_iters <= max_num_iters; n_iters++) {
295 CeedInt i, j;
296 MaxEntry(A, N, max_idx_row, &i, &j);
297
298 // If A[i][j] is small compared to A[i][i] and A[j][j], set it to 0.
299 if ((A[i * N + i] + A[i * N + j] == A[i * N + i]) && (A[j * N + j] + A[i * N + j] == A[j * N + j])) {
300 A[i * N + j] = 0.0;
301 max_idx_row[i] = MaxEntryRow(A, N, i);
302 }
303
304 if (A[i * N + j] == 0.0) break;
305
306 CalcRot(A, N, i, j, rotmat_cst); // Calculate the parameters of the rotation matrix.
307 ApplyRot(A, N, i, j, max_idx_row, rotmat_cst); // Apply this rotation to the A matrix.
308 if (calc_evec) ApplyRotLeft(evec, N, i, j, rotmat_cst);
309 }
310
311 for (CeedInt i = 0; i < N; i++) eval[i] = A[i * N + i];
312
313 // Optional: Sort results by eigenvalue.
314 SortRows(eval, evec, N, sort_criteria);
315
316 if ((n_iters > max_num_iters) && (N > 1)) // If we exceeded max_num_iters,
317 return 0; // indicate an error occured.
318
319 return n_iters;
320 }
321
322 // @brief Interface to Diagonalize for 3x3 systems
Diagonalize3(CeedScalar A[3][3],CeedScalar eval[3],CeedScalar evec[3][3],CeedInt max_idx_row[3],SortCriteria sort_criteria,bool calc_evec,const CeedInt max_num_sweeps)323 CEED_QFUNCTION_HELPER CeedInt Diagonalize3(CeedScalar A[3][3], CeedScalar eval[3], CeedScalar evec[3][3], CeedInt max_idx_row[3],
324 SortCriteria sort_criteria, bool calc_evec, const CeedInt max_num_sweeps) {
325 return Diagonalize((CeedScalar *)A, 3, (CeedScalar *)eval, (CeedScalar *)evec, (CeedInt *)max_idx_row, sort_criteria, calc_evec, max_num_sweeps);
326 }
327