Anonymous user
LU decomposition: Difference between revisions
→{{header|Fortran}}
m (→{{header|zkl}}: remove verbage) |
|||
Line 826:
=={{header|Fortran}}==
<lang Fortran>
program lu1
implicit none
call check( reshape([real(8)::1,2,1,3,4,1,5,7,0 ],[3,3]) )
call check( reshape([real(8)::11,1,3,2,9,5,17,5,24,2,18,7,2,6,1,1],[4,4]) )
contains
subroutine check(a)
real(8), intent(in) :: a(:,:)
integer :: i,j,n
real(8), allocatable :: aa(:,:),l(:,:),u(:,:)
integer, allocatable :: p(:,:)
integer, allocatable :: ipiv(:)
n = size(a,1)
allocate(aa(n,n),l(n,n),u(n,n),p(n,n),ipiv(n))
forall (j=1:n,i=1:n)
aa(i,j) = a(i,j)
u (i,j) = 0d0
p (i,j) = merge(1 ,0 ,i.eq.j)
l (i,j) = merge(1d0,0d0,i.eq.j)
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
call mat_print('a',a)
call mat_print('p',p)
call mat_print('l',l)
call mat_print('u',u)
print *, "residual"
print *, "|| P.A - L.U || = ", maxval(abs(matmul(p,a)-matmul(l,u)))
end subroutine
subroutine lu(a,p)
! in situ decomposition, corresponds to LAPACK's dgebtrf
real(8), intent(inout) :: a(:,:)
integer, intent(out ) :: p(:)
integer :: n, i,j,k,kmax
n = size(a,1)
p = [ ( i, i=1,n ) ]
do k = 1,n-1
kmax = maxloc(abs(a(p(k:),k)),1) + k-1
if (kmax /= k ) p([k, kmax]) = p([kmax, k])
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 do
end subroutine
subroutine mat_print(amsg,a)
character(*), intent(in) :: amsg
class (*), intent(in) :: a(:,:)
integer :: i
print*,' '
print*,amsg
do i=1,size(a,1)
select type (a)
type is (real(8)) ; print'(100f8.2)',a(i,:)
type is (integer) ; print'(100i8 )',a(i,:)
end select
end do
print*,' '
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>
<!--
<lang Fortran>
program lu1
implicit none
Line 901 ⟶ 1,017:
end subroutine
end program></lang>
--!>
=={{header|Go}}==
| |||