program matrixSAD

       IMPLICIT NONE
       

!modif (david,2013) pour faire passer modele antedependant
	Integer :: io,ef,nb,n, eff,i,j
	integer :: x(40)
   character :: xc(40)*40, parfile*40
   real :: xx(40)
	DOUBLE PRECISION, allocatable :: ginv(:,:),gante(:,:)
integer,parameter::   io_p=40          ! unit number for parameter file
	




write(*,'(a)',advance='no')' name of .gdg file?'
   read '(a)',parfile
   parfile=adjustl(parfile)               ! ignore leading spaces
   open(io_p,file=parfile)
call readline(io_p,x,xc,n)
   if (n == -1) then
       print*,'INPUT FILE EMPTY'
       stop
   endif
   nb=x(1)  !pour avoir la taille de la matrice
call readline(io_p,x,xc,n)
allocate(ginv(nb,nb),gante(nb,nb))
ginv=0.0d0
gante=0.0d0
print*,'lecture valeur'
 do ef=1,nb*nb
      read(io_p,'(I5,I5,G18.10)')i,j,xx(1)
      ginv(i,j)=xx(1)
      ginv(j,i)=xx(1)
      if (i.eq.nb.and.j.eq.nb) then
      exit
      endif
 enddo
 
 
print*,'inverse of the covariance matrix'
do i=1,nb
write(*,'(40f16.4)'),ginv(i,:)
enddo
 
 call inverse(ginv,gante,nb)
  
open(22,file='resumat')
print*,'covariance matrix'
write(22, fmt=*),'covariance matrix'
do i=1,nb
write(*,'(40f16.4)'),gante(i,:)
write(22,'(40f18.10)'),gante(i,:)
enddo


do ef=2,nb
do eff=1,ef-1
gante(ef,eff)=gante(ef,eff)/sqrt(gante(ef,ef)*gante(eff,eff))
gante(eff,ef)=gante(ef,eff)
enddo
enddo
print*,'correlation matrix (variances on the diagonal)'
write(22, fmt=*),'correlation matrix (variances on the diagonal)'
do i=1,nb
write(*,'(40f16.4)'),gante(i,:)
write(22,'(40f18.10)'),gante(i,:)
enddo


contains
subroutine inverse(a,c,n)
!============================================================
! Inverse matrix
! Method: Based on Doolittle LU factorization for Ax=b
! Alex G. December 2009
!-----------------------------------------------------------
! input ...
! a(n,n) - array of coefficients for matrix A
! n      - dimension
! output ...
! c(n,n) - inverse matrix of A
! comments ...
! the original matrix a(n,n) will be destroyed 
! during the calculation
!===========================================================
implicit none 
integer n
double precision a(n,n), c(n,n)
double precision L(n,n), U(n,n), b(n), d(n), x(n)
double precision coeff
integer i, j, k

! step 0: initialization for matrices L and U and b
! Fortran 90/95 aloows such operations on matrices
L=0.0
U=0.0
b=0.0

! step 1: forward elimination
do k=1, n-1
   do i=k+1,n
      coeff=a(i,k)/a(k,k)
      L(i,k) = coeff
      do j=k+1,n
         a(i,j) = a(i,j)-coeff*a(k,j)
      end do
   end do
end do

! Step 2: prepare L and U matrices 
! L matrix is a matrix of the elimination coefficient
! + the diagonal elements are 1.0
do i=1,n
  L(i,i) = 1.0
end do
! U matrix is the upper triangular part of A
do j=1,n
  do i=1,j
    U(i,j) = a(i,j)
  end do
end do

! Step 3: compute columns of the inverse matrix C
do k=1,n
  b(k)=1.0
  d(1) = b(1)
! Step 3a: Solve Ld=b using the forward substitution
  do i=2,n
    d(i)=b(i)
    do j=1,i-1
      d(i) = d(i) - L(i,j)*d(j)
    end do
  end do
! Step 3b: Solve Ux=d using the back substitution
  x(n)=d(n)/U(n,n)
  do i = n-1,1,-1
    x(i) = d(i)
    do j=n,i+1,-1
      x(i)=x(i)-U(i,j)*x(j)
    end do
    x(i) = x(i)/u(i,i)
  end do
! Step 3c: fill the solutions x(n) into column k of C
  do i=1,n
    c(i,k) = x(i)
  end do
  b(k)=0.0
end do
end subroutine inverse



subroutine readline(unit,x,xc,n,nohash)
! Reads a line from unit and decomposes it into numeric fields in x and
! alphanumeric fields in xc; the number of fields is stored in n.
! If optional variable nohash is absent, 
!  all characters after character # are ignored.

  integer unit,x(40),n,stat,i
  integer,optional::nohash
  character xc(40)*(*), a*200
 
  do
     read(unit,'(a)',iostat=stat)a
     if (stat /= 0) then
         n=-1
         return
     endif
     a=adjustl(a)               ! ignore leading spaces
     
     if (present(nohash)) then 
        exit
       else	!ignore characters after #
        i=scan(a,'#')
        if (i == 0) exit
        if (i>1) then
            a=a(1:i-1)
            exit
        endif
     endif  	
  enddo
  
  call nums2(a,n,x,xc)
  return
end subroutine


subroutine nums2(a,n,x,xc)
! separates array a into items delimited by blanks. character elements are
! put into optional character vector xc, decoded numeric values 
! into optional real vector x, and n contains the number of items. The 
! dimension of x and xc can be lower than n.
! A modification of nums() from f77 to f90
! Now accepts real numbers
! 2/23/2000

character (*)::a
character (*),optional::xc(:)
integer,optional::x(:)
integer::n,curr,first,last,lena,stat,i

curr=1;lena=len(a);n=0

do 
  ! search for first nonspace
  first=0
  do i=curr,lena
     if (a(i:i) /= ' ') then
        first=i
	exit
     endif
  enddo
  if (first == 0) exit
  
  
  ! search for first space
  curr=first+1
  last=0
  do i=curr,lena
      if (a(i:i) == ' ') then
        last=i
        exit
      endif
  enddo
  
  if (last == 0) last=lena
  
  n=n+1
  if (present(xc)) then
     if (size(xc) >= n) xc(n)=a(first:last)
  endif   
  if (present(x)) then
     if (size(x) >= n) then
        read(a(first:last),'(i12)',iostat=stat)x(n)
        if (stat /=0) x(n)=0
     endif
  endif   
	
  curr=last+1
enddo
end subroutine


end program matrixSAD





