Lines Matching +full:- +full:j
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
33 /// @returns Index of absolute largest off-diagonal element in row i
36 for (CeedInt j = i + 2; j < N; j++) in MaxEntryRow() local
37 if (fabs(A[i * N + j]) > fabs(A[i * N + j_max])) j_max = j; in MaxEntryRow()
52 for (CeedInt i = 1; i < N - 1; i++) { in MaxEntry()
53 CeedInt j = max_idx_row[i]; in MaxEntry() local
54 if (fabs(A[i * N + j]) > max_entry) { in MaxEntry()
55 max_entry = fabs(A[i * N + j]); in MaxEntry()
57 *j_max = j; in MaxEntry()
63 /// rotation in the i,j plane by an angle (θ) that (when multiplied on
65 /// A[i][j] = 0. The results will be stored in c, s, and t
70 /// @param[in] j column index
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()
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
96 /// It also assumes that i<j. The max_idx_row[] array is also updated.
101 /// where R the rotation in the i,j plane and ^T denotes the transpose.
102 /// i j
112 /// | -s ... c |
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].
133 /// matrix elements in the upper-right triangle strictly above the diagonal.
138 /// i j
149 /// | X X X X X | j
158 /// @param[in] j column index
159 CEED_QFUNCTION_HELPER void ApplyRot(CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedInt *max_id… in ApplyRot() argument
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()
172 …(CeedInt w = 0; w < i; w++) { // 0 <= w < i < j < N 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()
178 …for (CeedInt w = i + 1; w < j; w++) { // 0 <= i < w < … 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()
182 …for (CeedInt w = j + 1; w < N; w++) { // 0 <= i < j+1 <= … 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()
190 // compute A[w][j] and A[j][w] for all w!=j,considering above-diagonal elements in ApplyRot()
191 …CeedInt w = 0; w < i; w++) { // 0 <= w < i < j < N 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()
196 …for (CeedInt w = i + 1; w < j; w++) { // 0 <= i+1 <= w < … 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()
201 …for (CeedInt w = j + 1; w < N; w++) { // 0 <= i < 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()
204 // now that we're done modifying row j, we can update max_idx_row[j] in ApplyRot()
205 max_idx_row[j] = MaxEntryRow(A, N, j); in ApplyRot()
209 /// This matrix performs a rotation in the i,j plane by angle θ (where
217 /// @param[in] j column index
218 CEED_QFUNCTION_HELPER void ApplyRotLeft(CeedScalar *A, CeedInt N, CeedInt i, CeedInt j, CeedScalar … in ApplyRotLeft() argument
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()
238 for (CeedInt j = i + 1; j < N; j++) { in SortRows() local
242 if (eval[j] > eval[i_max]) i_max = j; in SortRows()
245 if (eval[j] < eval[i_max]) i_max = j; in SortRows()
248 if (fabs(eval[j]) > fabs(eval[i_max])) i_max = j; in SortRows()
251 if (fabs(eval[j]) < fabs(eval[i_max])) i_max = j; in SortRows()
270 /// on it, so divergence from normalized is due to finite-precision
280 … max_num_sweeps maximum number of iterations = max_num_sweeps * number of off-diagonals (N*(N-1)/2)
287 …for (CeedInt j = 0; j < N; j++) evec[i * N + j] = (i == j) ? 1.0 : 0.0; // Set evec equal to the … in Diagonalize() local
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()
295 CeedInt i, j; in Diagonalize() local
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()
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()
308 if (calc_evec) ApplyRotLeft(evec, N, i, j, rotmat_cst); in Diagonalize()