Anonymous user
LU decomposition: Difference between revisions
D implementation
m (add missing paren) |
(D implementation) |
||
Line 264:
#2A((11 9 24 2) (0 160/11 126/11 5/11) (0 0 -139/40 91/16) (0 0 0 71/139))
#2A((1 0 0 0) (0 0 1 0) (0 1 0 0) (0 0 0 1))</lang>
=={{header|D}}==
<lang d>import std.stdio, std.algorithm, std.typecons;
// plus functions from Matrix multiplication task
/// Creates a nxn identity matrix.
pure T[][] identityMatrix(T)(const int n)
in {
assert(n >= 0);
} body {
auto m = new typeof(return)(n, n);
foreach (i, row; m) {
row[] = 0;
row[i] = 1;
}
return m;
}
/// Swap two rows l and k of a 2D array matrix m.
pure void swapRows(T)(ref T[][] m, const int l, const int k) {
swap(m[l], m[k]);
}
/// Creates the pivoting matrix for m.
pure T[][] pivotize(T)(const T[][] m)
in {
// is square
assert(isRectangular(m) && m[0].length == m.length);
} body {
int n = m.length;
auto P = identityMatrix!T(n);
foreach (j; 0 .. n) {
T max = m[j][j];
int row = j;
foreach (i; j .. n) {
if (m[i][j] > max) {
max = m[i][j];
row = i;
}
}
if (j != row)
swapRows(P, j, row);
}
return P;
}
/// Decomposes a square matrix A by PA=LU and returns L, U and P.
/*pure*/ Tuple!(T[][], T[][], T[][]) lu(T)(const T[][] A)
in {
// is square
assert(isRectangular(A) && A[0].length == A.length);
} body {
int n = A.length;
auto L = new T[][](n, n);
auto U = new T[][](n, n);
foreach (i; 0 .. n) {
L[i][i .. $] = 0;
U[i][0 .. i] = 0;
}
auto P = pivotize!T(A);
auto A2 = matrixMul!T(P, A);
foreach (j; 0 .. n) {
L[j][j] = 1;
foreach (i; 0 .. j+1) {
T s1 = 0;
foreach (k; 0 .. i)
s1 += U[k][j] * L[i][k];
U[i][j] = A2[i][j] - s1;
}
foreach (i; j .. n) {
T s2 = 0;
foreach (k; 0 .. j)
s2 += U[k][j] * L[i][k];
L[i][j] = (A2[i][j] - s2) / U[j][j];
}
}
return tuple(L, U, P);
}
void main() {
enum double[][] a = [[1, 3, 5], [2, 4, 7], [1, 1, 0]];
auto lua = lu(a);
writeln(prettyPrint(lua[0]), "\n\n",
prettyPrint(lua[1]), "\n\n",
prettyPrint(lua[2]), "\n\n");
enum double[][] b = [[11, 9, 24, 2], [1, 5, 2, 6],
[3, 17, 18, 1], [2, 5, 7, 1]];
auto lub = lu(b);
writeln(prettyPrint(lub[0]), "\n\n",
prettyPrint(lub[1]), "\n\n",
prettyPrint(lub[2]));
}</lang>
Output:
<pre>[[1, 0, 0],
[0.5, 1, 0],
[0.5, -1, 1]]
[[2, 4, 7],
[0, 1, 1.5],
[0, 0, -2]]
[[0, 1, 0],
[1, 0, 0],
[0, 0, 1]]
[[1, 0, 0, 0],
[0.272727, 1, 0, 0],
[0.090909, 0.2875, 1, 0],
[0.181818, 0.23125, 0.00359712, 1]]
[[11, 9, 24, 2],
[0, 14.5455, 11.4545, 0.454545],
[0, 0, -3.475, 5.6875],
[0, 0, 0, 0.510791]]
[[1, 0, 0, 0],
[0, 0, 1, 0],
[0, 1, 0, 0],
[0, 0, 0, 1]]</pre>
| |||