program matso1
!      metodo di Cholesky
real, dimension(3,3)	:: a, c, d
real, dimension(3)		:: b, x, diag, h
!
call matdef(a,b,3)
write(*,*) 'A e B originali'
call matwrt(a,b,3)
!
call chodec(a,diag,3)
write(*,*) 'A dopo CHODEC e diagonale calcolata'
call matwrt(a,diag,3)
!
call chomat(a,diag,c,3)
write(*,*) 'matrice L ricostruita e vettore B originale'
call matwrt(c,b,3)
!
call chtes1(c,d,3)
write(*,*) 'A=L.LT e vettore B originale'
call matwrt(d,b,3)
!
call chosol(a,b,x,diag,3)
write(*,*) 'A dopo CHOSOL e vettore soluzione X'
call matwrt(a,x,3)
!
call chtes2(d,x,h,3)
write(*,*) 'A=L.LT e vettore B=L.LT.X'
call matwrt(d,h,3)
!
stop
end program
 
subroutine chodec(a,diag,n)
implicit none
integer					:: n, i, j, k
real					:: sum
real, dimension(n,n)	:: a
real, dimension(n)		:: diag
!
do i=1,n
!
  do j=i,n
!
    sum = a(i,j)
    do k = i-1,1,-1
      sum = sum - a(i,k) * a(j,k)
    enddo
    if (i .eq. j) then
      if (sum .le. 0.) pause 'matrice non simetrica e positiva'
      diag(i) = sqrt(sum)
    else
      a(j,i) = sum / diag(i)
    endif
!
  enddo
!
enddo
!
return
end subroutine

subroutine chtes1(c,d,n)
implicit none
integer					:: n, i, j, k
real					:: sum
real, dimension(n,n)	:: c, d
do i=1,n
  do j = 1,n
    sum = 0.
    do k=1,n
      sum = sum + c(i,k) * c(j,k)
    enddo
    d(i,j) = sum
  enddo
enddo
return
end subroutine

subroutine chosol(a,b,x,diag,n)
implicit none
integer					:: n, i, k
real					:: sum
real, dimension(n,n)	:: a
real, dimension(n)		:: b, diag, x
!
do i=1,n
  sum = b(i)
  do k=i-1,1,-1
    sum = sum - a(i,k) * x(k)
  enddo 
  x(i) = sum / diag(i)
enddo
!
do i=n,1,-1
  sum = x(i)
  do k = i+1,n
    sum = sum - a(k,i) * x(k)
  enddo
  x(i) = sum / diag(i)
enddo
!
return
end subroutine

subroutine chomat(a,diag,c,n)
implicit none
integer					:: n, i, j
real, dimension(n,n)	:: a, c
real, dimension(n)		:: diag
!
do i=1,n
  do j=1,n
    if (i .gt. j) then
       c(i,j) = a(i,j)
    else if (i.eq. j) then
       c(i,j) = diag(i)
    else
       c(i,j) = 0.
    endif
  enddo
enddo
!
return
end subroutine

subroutine chtes2(a,x,ax,n)
implicit none
integer					:: n, i, j, k
real, dimension(n,n)	:: a
real, dimension(n)		:: x, ax
real					:: sum
!       
do i=1,n
  sum = 0.
  do j=1,n
    sum = sum + a(i,j) * x(j)
  enddo
  ax(i) = sum
enddo
!
return
end subroutine

subroutine matdef(a,b,n)
implicit none
integer					:: n, i, j
real, dimension(n,n)	:: a
real, dimension(n)		:: b
!       
do i=1,n
  do j=1,n
     write(*,*) 'elemento di matrice ',i,',',j
     read(*,*) a(i,j)
  enddo
enddo
!
do i=1,n
  write(*,*) 'elemento del vettore ',i
  read(*,*) b(i)
enddo
!
return
end subroutine

subroutine matwrt(a,b,n)
implicit none
integer					:: n, i, j
real, dimension(n,n)	:: a
real, dimension(n)		:: b
!
do i=1,n
  write(*,'(1X,5f12.4)') (a(i,j), j=1,n)
enddo 
do i=1,n
  write(*,'(1X,f12.4)') b(i)
enddo
return
end subroutine
