Cholesky decomposition written in Fortran

cholesky_d.f

! cholesky_d.f -*-f90-*-
! Using Cholesky decomposition, cholesky_d.f solve a linear equation Ax=b,
! where A is a n by n positive definite real symmetric matrix, x and b are
! real*8 vectors length n.
!
! Time-stamp: <2015-06-25 18:05:47 takeshi>
! Author: Takeshi NISHIMATSU
! Licence: GPLv3
!
! [1] A = G tG, where G is a lower triangular matrix and tG is transpose of G.
! [2] Solve  Gy=b with  forward elimination
! [3] Solve tGx=y with backward elimination
!
! Reference: Taketomo MITSUI: Solvers for linear equations [in Japanese]
!            http://www2.math.human.nagoya-u.ac.jp/~mitsui/syllabi/sis/info_math4_chap2.pdf
!
! Comment:   This Cholesky decomposition is used in src/elastic.F and
!            src/optimize-inho-strain.F of feram http://loto.sourceforge.net/feram/ .
!!
subroutine cholesky_d(n, A, G, b)
  implicit none
  integer,    intent(in)    :: n
  real*8,     intent(in)    :: A(n,n)
  real*8,     intent(out)   :: G(n,n)
  real*8, intent(inout) :: b(n)
  real*8 :: tmp
  integer    :: i,j

  ! Light check of positive definite
  do i = 1, n
     if (A(i,i).le.0.0d0) then
        b(:) = 0.0d0
        return
     end if
  end do

  ! [1]
  G(:,:)=0.0d0
  do j = 1, n
     G(j,j) = sqrt( A(j,j) - dot_product(G(j,1:j-1),G(j,1:j-1)) )
     do i = j+1, n
        G(i,j)  = ( A(i,j) - dot_product(G(i,1:j-1),G(j,1:j-1)) ) / G(j,j)
     end do
  end do
  ! write(6,'(3f10.5)') (G(i,:), i=1,n)

  ! [2]
  do i = 1, n
     tmp = 0.0d0
     do j = 1, i-1
        tmp = tmp + G(i,j)*b(j)
     end do
     b(i) = (b(i)-tmp)/G(i,i)
  end do

  ! [3]
  do i = n, 1, -1
     tmp = 0.0d0
     do j = i+1, n
        tmp = tmp + b(j)*G(j,i)
     end do
     b(i) = (b(i)-tmp)/G(i,i)
  end do
end subroutine cholesky_d

cholesky_test.f

! cholesky_test.f -*-f90-*-
!!
program cholesky_test
  implicit none
  real*8 :: A(5,5)
  real*8 :: G(5,5)
  real*8 :: b(5)
  A = reshape((/ 5.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0, &
       &         4.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0, &
       &         3.0d0, 3.0d0, 3.0d0, 2.0d0, 1.0d0, &
       &         2.0d0, 2.0d0, 2.0d0, 2.0d0, 1.0d0, &
       &         1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0/), (/5,5/))
  b = (/ 0.5d0, 0.1d0, 0.8d0, 0.3d0, 0.6d0 /)
  call cholesky_d(5, A, G, b)
  write(6,'(5f6.2)') b
  write(6,'(5f6.2)') matmul(A,b)
end program cholesky_test
!Local variables:
!  compile-command: "gfortran -ffree-form -c cholesky_d.f && gfortran -ffree-form -c cholesky_test.f && gfortran cholesky_d.o cholesky_test.o -o cholesky_test && ./cholesky_test"
!End:

example.txt

$ gfortran -ffree-form -c cholesky_d.f
$ gfortran -ffree-form -c cholesky_test.f
$ gfortran cholesky_d.o cholesky_test.o -o cholesky_test
$ ./cholesky_test
  0.40 -1.10  1.20 -0.80  0.90
  0.50  0.10  0.80  0.30  0.60