Anonymous user
LU decomposition: Difference between revisions
Fortran added
(Racket) |
(Fortran added) |
||
Line 816:
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0]]</pre>
=={{header|Fortran}}==
<lang fortran>program lu1
implicit none
real*8 :: a1(3,3), a2(4,4)
a1 = reshape((/real*8::1,2,1,3,4,1,5,7,0/),(/3,3/))
call check(a1)
a2 = reshape((/real*8::11,1,3,2,9,5,17,5,24,2,18,7,2,6,1,1/),(/4,4/))
call check(a2)
contains
subroutine lu(a,p)
! in situ decomposition, correspondes to LAPACK's dgebtrf
implicit none
real*8, intent(inout) :: a(:,:)
integer, intent(out) :: p(:)
integer :: n, i,j,k,ii
n = ubound(a,1)
p = (/ ( i, i=1,n ) /)
do k = 1,n-1
ii = k-1+maxloc(abs(a(p(k:),k)),1)
if (ii /= k ) then
p((/k, ii/)) = p((/ii, k/))
end if
a(p(k+1:),k) = a(p(k+1:),k)/a(p(k),k)
forall (j = k+1:n)
a(p(k+1:),j) = a(p(k+1:),j)-a(p(k+1:),k)*a(p(k),j)
end forall
end do
end subroutine
subroutine check(a)
implicit none
real*8, intent(in) :: a(:,:)
real*8 :: aa(ubound(a,1), ubound(a,2))
real*8 :: l(ubound(a,1), ubound(a,2))
real*8 :: u(ubound(a,1), ubound(a,2))
integer :: p(ubound(a,1), ubound(a,2)), ipiv(ubound(a,1))
integer :: i, n
character(len=100) :: fmt
n = ubound(a,1)
aa = a ! work with a copy
p = 0; l=0; u = 0
forall (i=1:n)
p(i,i) = 1d0; l(i,i) = 1d0 ! convert permutation vector a matrix
end forall
call lu(aa, ipiv)
do i = 1,n
l(i,:i-1) = aa(ipiv(i),:i-1)
u(i,i:) = aa(ipiv(i),i:)
end do
p(ipiv,:) = p
write (fmt,"(a,i1,a)") "(",n,"(f8.2,1x))"
print *, "a"
print fmt, transpose(a)
print *, "p"
print fmt, transpose(dble(p))
print *, "l"
print fmt, transpose(l)
print *, "u"
print fmt, transpose(u)
print *, "residual"
print *, "|| P.A - L.U || = ", maxval(abs(matmul(p,a)-matmul(l,u)))
end subroutine
end program></lang>
{{out}}
<pre>
a
1.00 3.00 5.00
2.00 4.00 7.00
1.00 1.00 0.00
p
0.00 1.00 0.00
1.00 0.00 0.00
0.00 0.00 1.00
l
1.00 0.00 0.00
0.50 1.00 0.00
0.50 -1.00 1.00
u
2.00 4.00 7.00
0.00 1.00 1.50
0.00 0.00 -2.00
residual
|| P.A - L.U || = 0.0000000000000000
a
11.00 9.00 24.00 2.00
1.00 5.00 2.00 6.00
3.00 17.00 18.00 1.00
2.00 5.00 7.00 1.00
p
1.00 0.00 0.00 0.00
0.00 0.00 1.00 0.00
0.00 1.00 0.00 0.00
0.00 0.00 0.00 1.00
l
1.00 0.00 0.00 0.00
0.27 1.00 0.00 0.00
0.09 0.29 1.00 0.00
0.18 0.23 0.00 1.00
u
11.00 9.00 24.00 2.00
0.00 14.55 11.45 0.45
0.00 0.00 -3.47 5.69
0.00 0.00 0.00 0.51
residual
|| P.A - L.U || = 0.0000000000000000
</pre>
=={{header|Go}}==
| |||