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