40 for (
register label i=0; i<n; i++)
42 scalar largestCoeff = 0.0;
44 const scalar* __restrict__ matrixi = matrix[i];
46 for (
register label j=0; j<n; j++)
48 if ((temp =
mag(matrixi[j])) > largestCoeff)
54 if (largestCoeff == 0.0)
59 "(scalarSquareMatrix& matrix, labelList& rowIndices)"
63 vv[i] = 1.0/largestCoeff;
66 for (
register label j=0; j<n; j++)
68 scalar* __restrict__ matrixj = matrix[j];
70 for (
register label i=0; i<j; i++)
72 scalar* __restrict__ matrixi = matrix[i];
74 scalar
sum = matrixi[j];
75 for (
register label
k=0;
k<i;
k++)
77 sum -= matrixi[
k]*matrix[
k][j];
84 scalar largestCoeff = 0.0;
85 for (
register label i=j; i<n; i++)
87 scalar* __restrict__ matrixi = matrix[i];
88 scalar
sum = matrixi[j];
90 for (
register label
k=0;
k<j;
k++)
92 sum -= matrixi[
k]*matrix[
k][j];
98 if ((temp = vv[i]*
mag(sum)) >= largestCoeff)
105 pivotIndices[j] = iMax;
109 scalar* __restrict__ matrixiMax = matrix[iMax];
111 for (
register label
k=0;
k<n;
k++)
113 Swap(matrixj[
k], matrixiMax[k]);
119 if (matrixj[j] == 0.0)
126 scalar rDiag = 1.0/matrixj[j];
128 for (
register label i=j+1; i<n; i++)
130 matrix[i][j] *= rDiag;
151 "scalarRectangularMatrix& answer "
152 "const scalarRectangularMatrix& A, "
153 "const scalarRectangularMatrix& B)"
154 ) <<
"A and B must have identical inner dimensions but A.m = "
155 << A.m() <<
" and B.n = " << B.n()
161 for(
register label i = 0; i < A.n(); i++)
163 for(
register label j = 0; j < B.m(); j++)
165 for(
register label l = 0; l < B.n(); l++)
167 ans[i][j] += A[i][l]*B[l][j];
187 "const scalarRectangularMatrix& A, "
188 "const scalarRectangularMatrix& B, "
189 "const scalarRectangularMatrix& C, "
190 "scalarRectangularMatrix& answer)"
191 ) <<
"A and B must have identical inner dimensions but A.m = "
192 << A.m() <<
" and B.n = " << B.n()
201 "const scalarRectangularMatrix& A, "
202 "const scalarRectangularMatrix& B, "
203 "const scalarRectangularMatrix& C, "
204 "scalarRectangularMatrix& answer)"
205 ) <<
"B and C must have identical inner dimensions but B.m = "
206 << B.m() <<
" and C.n = " << C.n()
212 for(
register label i = 0; i < A.n(); i++)
214 for(
register label g = 0; g < C.m(); g++)
216 for(
register label l = 0; l < C.n(); l++)
219 for(
register label j = 0; j < A.m(); j++)
221 ab += A[i][j]*B[j][l];
223 ans[i][g] += C[l][g] * ab;
234 const DiagonalMatrix<scalar>& B,
238 if (A.m() != B.size())
243 "const scalarRectangularMatrix& A, "
244 "const DiagonalMatrix<scalar>& B, "
245 "const scalarRectangularMatrix& C, "
246 "scalarRectangularMatrix& answer)"
247 ) <<
"A and B must have identical inner dimensions but A.m = "
248 << A.m() <<
" and B.n = " << B.size()
252 if (B.size() != C.n())
257 "const scalarRectangularMatrix& A, "
258 "const DiagonalMatrix<scalar>& B, "
259 "const scalarRectangularMatrix& C, "
260 "scalarRectangularMatrix& answer)"
261 ) <<
"B and C must have identical inner dimensions but B.m = "
262 << B.size() <<
" and C.n = " << C.n()
268 for(
register label i = 0; i < A.n(); i++)
270 for(
register label g = 0; g < C.m(); g++)
272 for(
register label l = 0; l < C.n(); l++)
274 ans[i][g] += C[l][g] * A[i][l]*B[l];
287 SVD svd(A, minCondition);
288 return svd.VSinvUt();