yyhxwhd的个人博客分享 http://blog.sciencenet.cn/u/yyhxwhd

博文

Cholesky decomposition for positive definite matrix

已有 2101 次阅读 2014-9-21 22:09 |个人分类:Program;Script|系统分类:科研笔记

!  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+

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

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

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部