LU decomposition: Difference between revisions

Content added Content deleted
(Updated D code)
(Improved D code)
Line 589: Line 589:
std.array, std.conv;
std.array, std.conv;


bool isRectangular(T)(in T[][] M) pure nothrow {
bool isRectangular(T)(in T[][] m) pure nothrow {
//return !canFind!((r){ return r.length != M[0].length; })(M);
//return !canFind!((r){ return r.length != m[0].length; })(m);
foreach (row; M)
foreach (row; m)
if (row.length != M[0].length)
if (row.length != m[0].length)
return false;
return false;
return true;
return true;
}
}


string prettyPrint(T)(in T[][] M) {
bool isSquare(T)(in T[][] m) pure nothrow {
return "[" ~ array(map!text(M)).join(",\n ") ~ "]";
return isRectangular(m) && m[0].length == m.length;
}

string prettyPrint(T)(in T[][] m) {
return "[" ~ array(map!text(m)).join(",\n ") ~ "]";
}
}


Line 622: Line 626:
in {
in {
assert(n >= 0);
assert(n >= 0);
} out(result) {
assert(isSquare(result));
} body {
} body {
auto m = new typeof(return)(n, n);
auto m = new typeof(return)(n, n);
Line 634: Line 640:
T[][] pivotize(T)(in T[][] m) pure nothrow
T[][] pivotize(T)(in T[][] m) pure nothrow
in {
in {
// is square
assert(isSquare(m));
assert(isRectangular(m) && m[0].length == m.length);
} body {
} body {
immutable int n = m.length;
immutable n = m.length;
auto P = identityMatrix!T(n);
auto P = identityMatrix!T(n);


foreach (i; 0 .. n) {
foreach (i; 0 .. n) {
T max = m[i][i];
T max = m[i][i];
int row = i;
auto row = i;
foreach (j; i .. n)
foreach (j; i .. n)
if (m[j][i] > max) {
if (m[j][i] > max) {
Line 659: Line 664:
lu(T)(in T[][] A) pure nothrow
lu(T)(in T[][] A) pure nothrow
in {
in {
// is square
assert(isSquare(A));
assert(isRectangular(A) && A[0].length == A.length);
} body {
} body {
immutable int n = A.length;
immutable n = A.length;
auto L = new T[][](n, n);
auto L = new T[][](n, n);
auto U = new T[][](n, n);
auto U = new T[][](n, n);
Line 669: Line 673:
U[i][0 .. i] = 0;
U[i][0 .. i] = 0;
}
}

⚫
const auto P = pivotize!T(A);
auto A2 = matrixMul!T(P, A);
/*immutable*/ const P = pivotize!T(A);
⚫
/*immutable*/ const A2 = matrixMul!T(P, A);


foreach (j; 0 .. n) {
foreach (j; 0 .. n) {
Retrieved from "https://rosettacode.org/wiki/LU_decomposition"