88 lines
2.3 KiB
Plaintext
88 lines
2.3 KiB
Plaintext
(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;
|