Jump to content

LU decomposition: Difference between revisions

LU decomposition (view source)

Revision as of 17:33, 29 July 2012

464 bytes removed ,  11 years ago
Udated D entry
m (→‎{{header|REXX}}: updated the '''output'''s with new version (trailing zeroes suppressed). -- ~~~~)
(Udated D entry)
Line 587:
{{trans|Common Lisp}}
<lang d>import std.stdio, std.algorithm, std.typecons, std.numeric,
std.array, std.conv, std.string, std.range;
 
bool isRectangular(T)(in T[][] m) pure /*nothrow*/ {
Line 598:
}
 
immutable(T[][]) matrixMul(T)(immutablein T[][] A, in T[][] B) pure nothrow
in {▼
immutable T[][] B) pure nothrow
assert(A[0].length ==&& B.length);▼
▲ in {
assert(isRectangular(A.length) && isRectangular(B.length &&));
assert(A[0].length == B.length);
isRectangular(A) && isRectangular(B) &&
} body {▼
▲ A[0].length == B.length);
auto result = new T[][](A.length, B[0].length);
▲ } body {
auto resultaux = new T[][](AB.length, B[0].length);
foreach (j; 0 .. auto aux = new TB[B0].length];) {
foreach (jk, row; 0 .. B[0].length) {
foreach (aux[k,] = row[j]; B)
foreach (i, ai; aux[k] = row[j];A)
foreachresult[i][j] = dotProduct(iai, ai; Aaux);
result[i][j] = dotProduct(ai, aux);▼
}▼
return result;▼
}▼
 
/// Creates a nxn identity matrix.
T[][] identityMatrix(T)(in int n) pure nothrow
in {
assert(n >= 0);▼
} out(result) {
assert(isSquare(result));
} body {
auto m = new typeof(return)(n, n);
foreach (i, row; m) {▼
row[] = 0;▼
row[i] = 1;▼
}
return m;▼
}
▲ return result;
}
 
/// Creates the pivoting matrix for m.
immutable(T[][]) pivotize(T)(immutable T[][] m) /*pure nothrow*/
in {
assert(isSquare(m));
} body {
immutableconst n = m.length;
auto Pid = identityMatrix!Tiota(n);
.map!(j=> iota(n).map!(i => cast(T)(i == j))().array())()
▲ row[] = 0.array();
 
foreach (i; 0 .. n) {
// immutable row = Tiota(i, n).max!(j => m[ij][i])();
T auto rowmaxm = m[i][i];
auto row = foreach (ji; i .. n)
ifforeach (m[j][; i] >.. maxn) {
if max = (m[j][i]; > maxm) {
rowmaxm = m[j][i];
}row = j;
if (i != row)}
swap(P[i], P[row]);▼
}
 
returnif P;(i != row)
▲ rowswap(id[i], = 1id[row]);
}
 
▲ return mid;
}
 
/// Decomposes a square matrix A by PA=LU and returns L, U and P.
Tuple!(T[][],"L", T[][],"U", const T[][],"P")
lu(T)(immutable T[][] A) /*pure nothrow*/
in {
assert(isSquare(A));
} body {
immutable 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;
}
 
immutableconst P = pivotize!T(A);
immutableconst 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];
}
}
▲ foreach (i, row; mj .. n) {
▲ assert(n > T s2 = 0);
▲ foreach (ik; j0 .. nj) {
▲ results2 += U[ik][j] =* dotProduct(ai, aux)L[i][k];
▲ L[i][j] s2 += U(A2[ki][j] *- Ls2) / U[ij][kj];
▲ }
▲ }
 
return typeof(return)(L, U, P);
}
}
 
void main() {
immutable a = [[1.0, 3, 5],
[2.0, 4, 7],
[1.0, 1, 0]];
immutable b = [[11.0, 9, 24, 2],
foreach (part; lu(a))
writeln(" ["1.0, join(map!text(part), "5,\n ") 2, "6]\n");,
▲ T s2 = [3.0;, 17, 18, 1],
writeln();
▲ swap(P [i]2.0, P[row 5, 7, 1]]);
 
immutableauto bf = std.array.replicate("[%([11%(%.1f, 9%)],\n 24, 2%)]]\n\n", 3);
foreach (m; [1.a, 5, 2, 6b],)
writefln(f, [3lu(m)., 17, 18, 1],tupleof);
[2., 5, 7, 1]];
foreach (part; lu(b))
writeln("[", join(map!text(part), ",\n "), "]\n");
}</lang>
{{out}}
Output:
<pre>[[1.0, 0.0, 0.0],
[0.5, 1.0, 0.0],
[0.5, -1.0, 1.0]]
 
[[2.0, 4.0, 7.0],
[0.0, 1.0, 1.5],
[0.0, 0.0, -2.0]]
 
[[0.0, 1.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 0.0, 1.0]]
 
 
[[1.0, 0.0, 0.0, 0.0],
[0.2727273, 1.0, 0.0, 0.0],
[0.0909090, 0.28753, 1.0, 0.0],
[0.1818182, 0.231252, 0.003597120, 1.0]]
 
[[11.0, 9.0, 24.0, 2.0],
[0.0, 14.54555, 11.45455, 0.4545455],
[0.0, 0.0, -3.4755, 5.68757],
[0.0, 0.0, 0.0, 0.5107915]]
 
[[1.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0]]</pre>
</pre>
 
=={{header|Go}}==
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/175222"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information