Anonymous user
LU decomposition: Difference between revisions
Initial PL/I version
(Updated D entry) |
(Initial PL/I version) |
||
Line 1,531:
lu_backsub(lup, transpose([1, 1, -1, -1]));
/* matrix([-204], [2100], [-4740], [2940]) */</lang>
=={{header|PL/I}}==
<lang PL/I>(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;
</lang>
Derived from Fortran version above.
Results:
<pre>
A
1.00000 3.00000 5.00000
2.00000 4.00000 7.00000
1.00000 1.00000 0.00000
P
0 1 0
1 0 0
0 0 1
L
1.00000 0.00000 0.00000
0.50000 1.00000 0.00000
0.50000 -1.00000 1.00000
U
2.00000 4.00000 7.00000
0.00000 1.00000 1.50000
0.00000 0.00000 -2.00000
A
11.00000 9.00000 24.00000 2.00000
1.00000 5.00000 2.00000 6.00000
3.00000 17.00000 18.00000 1.00000
2.00000 5.00000 7.00000 1.00000
P
1 0 0 0
0 0 1 0
0 1 0 0
0 0 0 1
L
1.00000 0.00000 0.00000 0.00000
0.27273 1.00000 0.00000 0.00000
0.09091 0.28750 1.00000 0.00000
0.18182 0.23125 0.00360 1.00000
U
11.00000 9.00000 24.00000 2.00000
0.00000 14.54545 11.45455 0.45455
0.00000 0.00000 -3.47500 5.68750
0.00000 0.00000 0.00000 0.51079
</pre>
=={{header|Python}}==
| |||