Lines Matching +full:- +full:a
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.
4 // SPDX-License-Identifier: BSD-2-Clause
29 ///@brief Find the off-diagonal index in row i whose absolute value is largest
31 /// @param[in] *A matrix
33 /// @returns Index of absolute largest off-diagonal element in row i
34 CEED_QFUNCTION_HELPER CeedInt MaxEntryRow(const CeedScalar *A, CeedInt N, CeedInt i) { in MaxEntryRow() argument
37 if (fabs(A[i * N + j]) > fabs(A[i * N + j_max])) j_max = j; in MaxEntryRow()
45 /// @param[in] *A matrix
48 CEED_QFUNCTION_HELPER void MaxEntry(const CeedScalar *A, CeedInt N, CeedInt *max_idx_row, CeedInt *… in MaxEntry() argument
51 CeedScalar max_entry = fabs(A[*i_max * N + *j_max]); in MaxEntry()
52 for (CeedInt i = 1; i < N - 1; i++) { in MaxEntry()
54 if (fabs(A[i * N + j]) > max_entry) { in MaxEntry()
55 max_entry = fabs(A[i * N + j]); in MaxEntry()
62 /// @brief Calculate the components of a rotation matrix which performs a
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
68 /// @param[in] *A matrix
71 CEED_QFUNCTION_HELPER void CalcRot(const CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedScalar… in CalcRot() argument
73 CeedScalar A_jj_ii = (A[j * N + j] - A[i * N + i]); in CalcRot()
75 // kappa = (A[j][j] - A[i][i]) / (2*A[i][j]) in CalcRot()
78 CeedScalar A_ij = A[i * N + j]; in CalcRot()
81 // t satisfies: t^2 + 2*t*kappa - 1 = 0 in CalcRot()
84 if (kappa < 0.0) rotmat_cst[2] = -rotmat_cst[2]; in CalcRot()
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
98 /// below the diagonal (ie. A[u][v] where u>v) are not computed.
100 /// A' = R^T * A * R
112 /// | -s ... c |
119 /// Let A' denote the matrix A after multiplication by R^T and R.
120 /// The components of A' are:
123 /// A'_uv = Σ_w Σ_z R_wu * A_wz * R_zv
126 /// Note that a the rotation at location i,j will modify all of the matrix
128 /// such as: A[w][i], A[i][w], A[w][j], A[j][w].
133 /// matrix elements in the upper-right triangle strictly above the diagonal.
147 /// A = | . X |
156 /// @param[in] *A matrix
159 CEED_QFUNCTION_HELPER void ApplyRot(CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedInt *max_id… in ApplyRot() argument
160 // Compute the diagonal elements of A which have changed: in ApplyRot()
161 A[i * N + i] -= rotmat_cst[2] * A[i * N + j]; in ApplyRot()
162 A[j * N + j] += rotmat_cst[2] * A[i * N + j]; in ApplyRot()
164 // A[i][i] = c*c*A[i][i] + s*s*A[j][j] - 2*s*c*A[i][j] in ApplyRot()
165 // A[j][j] = s*s*A[i][i] + c*c*A[j][j] + 2*s*c*A[i][j] in ApplyRot()
167 // Update the off-diagonal elements of A which will change (above the diagonal) in ApplyRot()
169 A[i * N + j] = 0.0; in ApplyRot()
171 // compute A[w][i] and A[i][w] for all w!=i,considering above-diagonal elements in ApplyRot()
173 …A[i * N + w] = A[w * N + i]; // backup the previou… in ApplyRot()
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] f… in ApplyRot()
175 if (i == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(A, N, w); in ApplyRot()
176 else if (fabs(A[w * N + i]) > fabs(A[w * N + max_idx_row[w]])) max_idx_row[w] = i; in ApplyRot()
179 …A[w * N + i] = A[i * N + w]; // backup the previou… in ApplyRot()
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] f… in ApplyRot()
183 …A[w * N + i] = A[i * N + w]; // backup the previou… in ApplyRot()
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] f… in ApplyRot()
188 max_idx_row[i] = MaxEntryRow(A, N, i); in ApplyRot()
190 // compute A[w][j] and A[j][w] for all w!=j,considering above-diagonal elements in ApplyRot()
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] f… in ApplyRot()
193 if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(A, N, w); in ApplyRot()
194 else if (fabs(A[w * N + j]) > fabs(A[w * N + max_idx_row[w]])) max_idx_row[w] = j; in ApplyRot()
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] f… in ApplyRot()
198 if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(A, N, w); in ApplyRot()
199 else if (fabs(A[w * N + j]) > fabs(A[w * N + max_idx_row[w]])) max_idx_row[w] = j; in ApplyRot()
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] f… in ApplyRot()
205 max_idx_row[j] = MaxEntryRow(A, N, j); in ApplyRot()
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
212 /// A'_uv = Σ_w R_wu * A_wv
215 /// @param[in] *A matrix
218 CEED_QFUNCTION_HELPER void ApplyRotLeft(CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedScalar … in ApplyRotLeft() argument
221 CeedScalar Aiv = A[i * N + v]; in ApplyRotLeft()
222 A[i * N + v] = rotmat_cst[0] * A[i * N + v] - rotmat_cst[1] * A[j * N + v]; in ApplyRotLeft()
223 A[j * N + v] = rotmat_cst[1] * Aiv + rotmat_cst[0] * A[j * N + v]; in ApplyRotLeft()
236 for (CeedInt i = 0; i < N - 1; i++) { in SortRows()
262 /// @brief Calculate all the eigenvalues and eigevectors of a symmetric matrix
270 /// on it, so divergence from normalized is due to finite-precision
273 // @param[in] A the matrix you wish to diagonalize (size NxN)
280 … max_num_sweeps maximum number of iterations = max_num_sweeps * number of off-diagonals (N*(N-1)/2)
281 CEED_QFUNCTION_HELPER CeedInt Diagonalize(CeedScalar *A, CeedInt N, CeedScalar *eval, CeedScalar *e… in Diagonalize() argument
289 for (CeedInt i = 0; i < N - 1; i++) max_idx_row[i] = MaxEntryRow(A, N, i); in Diagonalize()
291 // -- Iteration -- in Diagonalize()
293 CeedInt max_num_iters = max_num_sweeps * N * (N - 1) / 2; in Diagonalize()
296 MaxEntry(A, N, max_idx_row, &i, &j); in Diagonalize()
298 // If A[i][j] is small compared to A[i][i] and A[j][j], set it to 0. in Diagonalize()
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])… in Diagonalize()
300 A[i * N + j] = 0.0; in Diagonalize()
301 max_idx_row[i] = MaxEntryRow(A, N, i); in Diagonalize()
304 if (A[i * N + j] == 0.0) break; in Diagonalize()
306 …CalcRot(A, N, i, j, rotmat_cst); // Calculate the parameters of the rotation matrix. in Diagonalize()
307 ApplyRot(A, N, i, j, max_idx_row, rotmat_cst); // Apply this rotation to the A matrix. in Diagonalize()
311 for (CeedInt i = 0; i < N; i++) eval[i] = A[i * N + i]; in Diagonalize()
323 CEED_QFUNCTION_HELPER CeedInt Diagonalize3(CeedScalar A[3][3], CeedScalar eval[3], CeedScalar evec[… in Diagonalize3()
325 …return Diagonalize((CeedScalar *)A, 3, (CeedScalar *)eval, (CeedScalar *)evec, (CeedInt *)max_idx_… in Diagonalize3()