Time for an 2014 update…
This commit is contained in:
parent
372c577f83
commit
09687c4926
2520 changed files with 34227 additions and 7318 deletions
|
|
@ -28,13 +28,13 @@ in {
|
|||
}
|
||||
|
||||
/// Creates the pivoting matrix for m.
|
||||
T[][] pivotize(T)(immutable T[][] m) /*pure nothrow*/
|
||||
T[][] pivotize(T)(immutable T[][] m) pure nothrow
|
||||
in {
|
||||
assert(m.isSquare);
|
||||
} body {
|
||||
immutable n = m.length;
|
||||
auto id = iota(n)
|
||||
.map!(j=> n.iota.map!(i => cast(T)(i == j)).array)
|
||||
.map!((in j) => n.iota.map!(i => cast(T)(i == j)).array)
|
||||
.array;
|
||||
|
||||
foreach (immutable i; 0 .. n) {
|
||||
|
|
@ -56,7 +56,7 @@ in {
|
|||
|
||||
/// 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*/
|
||||
lu(T)(immutable T[][] A) pure nothrow
|
||||
in {
|
||||
assert(A.isSquare);
|
||||
} body {
|
||||
|
|
@ -68,8 +68,8 @@ in {
|
|||
U[i][0 .. i] = 0;
|
||||
}
|
||||
|
||||
const P = A.pivotize!T;
|
||||
const A2 = matrixMul!T(P, A);
|
||||
immutable P = A.pivotize!T;
|
||||
immutable A2 = matrixMul!T(P, A);
|
||||
|
||||
foreach (immutable j; 0 .. n) {
|
||||
L[j][j] = 1;
|
||||
|
|
@ -99,8 +99,7 @@ void main() {
|
|||
[3.0, 17, 18, 1],
|
||||
[2.0, 5, 7, 1]];
|
||||
|
||||
//auto f = "[%([%(%.1f, %)],\n %)]]\n\n".replicate(3);
|
||||
auto f = std.array.replicate("[%([%(%.1f, %)],\n %)]]\n\n", 3);
|
||||
foreach (m; [a, b])
|
||||
auto f = "[%([%(%.1f, %)],\n %)]]\n\n".replicate(3);
|
||||
foreach (immutable m; [a, b])
|
||||
writefln(f, lu(m).tupleof);
|
||||
}
|
||||
|
|
|
|||
87
Task/LU-decomposition/PL-I/lu-decomposition.pli
Normal file
87
Task/LU-decomposition/PL-I/lu-decomposition.pli
Normal file
|
|
@ -0,0 +1,87 @@
|
|||
(subscriptrange, fofl, size): /* 2 Nov. 2013 */
|
||||
LU_Decomposition: procedure options (main);
|
||||
declare a1(3,3) float (18) initial ( 1, 3, 5,
|
||||
2, 4, 7,
|
||||
1, 1, 0);
|
||||
declare a2(4,4) float (18) initial (11, 9, 24, 2,
|
||||
1, 5, 2, 6,
|
||||
3, 17, 18, 1,
|
||||
2, 5, 7, 1);
|
||||
call check(a1);
|
||||
call check(a2);
|
||||
|
||||
|
||||
/* In-situ decomposition */
|
||||
LU: procedure(a, p);
|
||||
declare a(*,*) float (18);
|
||||
declare p(*) fixed binary;
|
||||
declare (maximum, rtemp) float (18);
|
||||
declare (n, i, j, k, ii, temp) fixed binary;
|
||||
|
||||
n = hbound(a,1);
|
||||
do i = 1 to n; p(i) = i; end;
|
||||
|
||||
do k = 1 to n-1;
|
||||
|
||||
maximum = 0; ii = k;
|
||||
do i = k to n;
|
||||
if maximum < abs(a(p(i),k)) then
|
||||
do; maximum = abs(a(p(i),k)); ii = i; end;
|
||||
end;
|
||||
if ii ^= k then do; temp = p(k); p(k) = p(ii); p(ii) = temp; end;
|
||||
|
||||
do i = k+1 to n; a(p(i),k) = a(p(i),k) / a(p(k),k); end;
|
||||
|
||||
do j = k+1 to n;
|
||||
do i = k+1 to n;
|
||||
a(p(i),j) = a(p(i),j) - a(p(i),k) * a(p(k),j);
|
||||
end;
|
||||
end;
|
||||
|
||||
end;
|
||||
end LU;
|
||||
|
||||
CHECK: procedure(a);
|
||||
declare a(*,*) float (18) nonassignable;
|
||||
|
||||
declare aa(hbound(a,1), hbound(a,2)) float (18);
|
||||
declare L(hbound(a,1), hbound(a,2)) float (18);
|
||||
declare U(hbound(a,1), hbound(a,2)) float (18);
|
||||
declare (p(hbound(a,1), hbound(a,2)), ipiv(hbound(a,1)) ) fixed binary;
|
||||
declare pp(hbound(a,1), hbound(a,2)) fixed binary;
|
||||
declare (i, j, n, temp(hbound(a,1))) fixed binary;
|
||||
|
||||
n = hbound(a,1);
|
||||
aa = A; /* work with a copy */
|
||||
P = 0; L = 0; U = 0;
|
||||
do i = 1 to n;
|
||||
p(i,i) = 1; L(i,i) = 1; /* convert permutation vector to a matrix */
|
||||
end;
|
||||
|
||||
call LU(aa, ipiv);
|
||||
|
||||
do i = 1 to n;
|
||||
do j = 1 to i-1; L(i,j) = aa(ipiv(i),j); end;
|
||||
do j = i to n; U(i,j) = aa(ipiv(i),j); end;
|
||||
end;
|
||||
|
||||
pp = p;
|
||||
do i = 1 to n;
|
||||
p(ipiv(i), *) = pp(i,*);
|
||||
end;
|
||||
|
||||
put skip list ('A');
|
||||
put edit (A) (skip, (n) f(10,5));
|
||||
|
||||
put skip list ('P');
|
||||
put edit (P) (skip, (n) f(11));
|
||||
|
||||
put skip list ('L');
|
||||
put edit (L) (skip, (n) f(10,5));
|
||||
|
||||
put skip list ('U');
|
||||
put edit (U) (skip, (n) f(10,5));
|
||||
|
||||
end CHECK;
|
||||
|
||||
end LU_Decomposition;
|
||||
|
|
@ -9,7 +9,7 @@ def pivotize(m):
|
|||
n = len(m)
|
||||
ID = [[float(i == j) for i in xrange(n)] for j in xrange(n)]
|
||||
for j in xrange(n):
|
||||
row = max(xrange(j, n), key=lambda i: m[i][j])
|
||||
row = max(xrange(j, n), key=lambda i: abs(m[i][j]))
|
||||
if j != row:
|
||||
ID[j], ID[row] = ID[row], ID[j]
|
||||
return ID
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue