Cholesky decomposition for positive definite matrix

!  Cholesky_decomposition.f90

!

!****************************************************************************

!

!  PROGRAM: Cholesky_decomposition

!

!  PURPOSE:  Use the Cholesky method to decopose a positive definite matrix ,S=L*L'

!  here L' is transpose matrix of L. We then obtain the inverse matrix of S and L

!  matrix as an output.

!

!  Written by Wang Haidi  at  2013-9-15 12:11:50

!****************************************************************************

program Cholesky_decomposition

implicit none

! Variables

integer::i,j,k

real(kind=8)::summ

real(kind=8)::S(3,3)=(/ 1.0,1.0,1.0,1.0,5.0,3.0,1.0,3.0,11.0 /)

real(kind=8)::L(3,3),L_inv(3,3),L_inv_T(3,3),S_inv(3,3)

real(kind=8),intrinsic::sqrt,transpose,matmul

! Body of Cholesky_decomposition for the positive definite matrix S

write(*,*),"S matrix"

write(*,"(3(f12.6,2x))"),S

!   L(1,1)=sqrt(S(1,1))

do i=1,3

do j=1,3

if(i>j)then

summ=0.0

do k=1,j-1

summ=summ+L(i,k)*L(j,k)

end do

L(i,j)=(S(i,j)-summ)/L(j,j)

elseif(i==j)then

summ=0.0

do k=1,i-1

summ=summ+L(i,k)**2

end do

L(i,j)=sqrt(S(i,j)-summ)

else

L(i,j)=0.0

end if

end do

end do

write(*,*)"L matrix"

write(*,"(3(f12.6,2x))"),((L(j,i),i=1,3),j=1,3)

!!calculate the inverse matrix of L matrix ,at the same time

!!  S matrix by be sovled with sqare root method

do i=1,3

do j=1,3

if(i==j)then

L_inv(i,j)=1/L(i,j)

elseif(i>j)then

summ=0.0

do k=1,i-1

summ=summ+L(i,k)*L_inv(k,j)

end do

L_inv(i,j)=-1*summ/L(i,i)

else

L_inv(i,j)=0.0

end if

end do

end do

write(*,*)"L_inv matrix"

write(*,"(3(f12.6,2x))"),((L_inv(j,i),i=1,3),j=1,3)

L_inv_T=transpose(L_inv)

S_inv=matmul(L_inv_T,L_inv)

write(*,*)"S_inv matrix"

write(*,"(3(f12.6,2x))"),((S_inv(j,i),i=1,3),j=1,3)

end program Cholesky_decomposition

https://wap.sciencenet.cn/blog-276702-829703.html

全部精选博文导读

GMT+8, 2022-1-19 17:54