LU decomposition: Difference between revisions
Content added Content deleted
(Updated D entry) |
(Initial PL/I version) |
||
Line 1,531: | Line 1,531: | ||
lu_backsub(lup, transpose([1, 1, -1, -1])); |
lu_backsub(lup, transpose([1, 1, -1, -1])); |
||
/* matrix([-204], [2100], [-4740], [2940]) */</lang> |
/* 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}}== |
=={{header|Python}}== |