博文
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
上一篇：Bandstructure calculation with empirical method
下一篇：MINIMAL BASIS STO-3G CALCULATION ON HEH+
