41 label m = tmpMatrix.
m();
44 for (label i = 0; i < m; ++i)
47 scalar largestCoeff =
mag(tmpMatrix(iMax, i));
50 for (label j = i + 1; j < m; ++j)
52 if (
mag(tmpMatrix(j, i)) > largestCoeff)
55 largestCoeff =
mag(tmpMatrix(iMax, i));
61 for (label
k = i;
k < m; ++
k)
63 Swap(tmpMatrix(i,
k), tmpMatrix(iMax,
k));
65 Swap(sourceSol[i], sourceSol[iMax]);
69 if (
mag(tmpMatrix(i, i)) < 1
e-20)
77 for (label j = i + 1; j < m; ++j)
79 sourceSol[j] -= sourceSol[i]*(tmpMatrix(j, i)/tmpMatrix(i, i));
81 for (label
k = m - 1;
k >= i; --
k)
84 tmpMatrix(i,
k)*tmpMatrix(j, i)/tmpMatrix(i, i);
90 for (label j = m - 1; j >= 0; --j)
94 for (label
k = j + 1;
k < m; ++
k)
96 ntempvec += tmpMatrix(j,
k)*sourceSol[
k];
99 sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix(j, j);
126 label m = luMatrix.
m();
130 for (label i = 0; i < m; ++i)
132 label ip = pivotIndices[i];
133 Type
sum = sourceSol[ip];
134 sourceSol[ip] = sourceSol[i];
135 const scalar* __restrict__ luMatrixi = luMatrix[i];
139 for (label j = ii - 1; j < i; ++j)
141 sum -= luMatrixi[j]*sourceSol[j];
152 for (label i = m - 1; i >= 0; --i)
154 Type
sum = sourceSol[i];
155 const scalar* __restrict__ luMatrixi = luMatrix[i];
157 for (label j = i + 1; j < m; ++j)
159 sum -= luMatrixi[j]*sourceSol[j];
162 sourceSol[i] =
sum/luMatrixi[i];
174 label m = luMatrix.
m();
178 for (label i = 0; i < m; ++i)
180 Type
sum = sourceSol[i];
181 const scalar* __restrict__ luMatrixi = luMatrix[i];
185 for (label j = ii - 1; j < i; ++j)
187 sum -= luMatrixi[j]*sourceSol[j];
195 sourceSol[i] =
sum/luMatrixi[i];
198 for (label i = m - 1; i >= 0; --i)
200 Type
sum = sourceSol[i];
201 const scalar* __restrict__ luMatrixi = luMatrix[i];
203 for (label j = i + 1; j < m; ++j)
205 sum -= luMatrixi[j]*sourceSol[j];
208 sourceSol[i] =
sum/luMatrixi[i];
238 template<
class Form,
class Type>
249 <<
"A and B must have identical inner dimensions but A.n = "
250 <<
A.n() <<
" and B.m = " <<
B.m()
256 for (label i = 0; i <
A.m(); ++i)
258 for (label j = 0; j <
B.n(); ++j)
260 for (label l = 0; l <
B.m(); ++l)
262 ans(i, j) +=
A(i, l)*
B(l, j);