Anonymous user
LU decomposition: Difference between revisions
D entry updated
m (Updated D entry) |
(D entry updated) |
||
Line 687:
bool isRectangular(T)(in T[][] m) /*pure nothrow*/ {
}
Line 695 ⟶ 694:
}
T[][] matrixMul(T)(in T[][] A, in T[][] B) /*pure
in {
assert(A.
▲ assert(A[0].length == B.length);
} body {
auto result = new T[][](A.length, B[0].length);
auto aux = new T[B.length];
foreach (immutable k, const row; B)
aux[k] = row[j];
foreach (immutable i, const ai; A)
result[i][j] = dotProduct(ai, aux);
}
return result;
}
Line 717:
assert(isSquare(m));
} body {
auto id = iota(n)
.map!(j=> n.iota
.array
foreach (immutable i; 0 .. n) {
// immutable row = iota(i, n).max!(j => m[j][i])();
T maxm = m[i][i];
foreach (immutable j; i .. n)
if (m[j][i] > maxm) {
maxm = m[j][i];
Line 748:
auto L = new T[][](n, n);
auto U = new T[][](n, n);
foreach (immutable i; 0 .. n) {
L[i][i .. $] = 0;
U[i][0 .. i] = 0;
Line 756:
const A2 = matrixMul!T(P, A);
foreach (immutable j; 0 .. n) {
L[j][j] = 1;
foreach (immutable i; 0 .. j+1) {
T s1 = 0;
foreach (immutable k; 0 .. i)
s1 += U[k][j] * L[i][k];
U[i][j] = A2[i][j] - s1;
}
foreach (immutable i; j .. n) {
T s2 = 0;
foreach (immutable k; 0 .. j)
s2 += U[k][j] * L[i][k];
L[i][j] = (A2[i][j] - s2) / U[j][j];
| |||