c
c  Program to test p.d. of a matrix
c
c   Updated: May 31, 2011  LRS
c
      parameter (nt=50,ntp=2500)
      common /mats/v(ntp),u(ntp),ss(nt),tt(nt),uu(nt)
      real*8 v,u,b,x,z,ss,tt
c
      open(10,file='vcvR.d',form='formatted',status='old')
      open(17,file='EiglogR.d',form='formatted',status='unknown')
c
c  read in a matrix, residual matrix for fertility traits
c
      v = 0.0
      mt = 16
      ss=0.0
      write(17,8001)
 8001 format(/'  Residual Matrix')
 10   read(10,1001,end=20)ir,ic,x
 1001 format(2i3,1x,f10.0)
      if(ir.gt.mt)go to 10
      if(ic.gt.mt)go to 10
      if(ir.lt.1)go to 20
      if(ir.eq.ic)then
         ss(ir)=dsqrt(x)
         m = ihmssf(ir,ic,mt)
         v(m)=1.0
       else
          m=ihmssf(ir,ic,mt)
          v(m)=x
       endif   
       write(17,3002)ir,ic,x
 3002  format(1x,2i4,1x,d20.10)
      go to 10
c
c  compute the eigen values of the matrix
c
 20   b=0.0
      m = 0
      do 12 ir=1,mt
         do 13 ic=ir,mt
            m=m+1
            v(m) = v(m)*ss(ir)*ss(ic)
 13      continue
 12   continue
c   
      call PDcheq(v,mt)
      close(10)
      open(14,file='vcvR2.d',form='formatted',status='unknown')
c
c  transform all matrices, so residual diagonals = 1
c
      uu = 0.0
      do 121 ir=1,mt
         m=ihmssf(ir,ir,mt)
         uu(ir) = 1.d0/(dsqrt(v(m)))
 121  continue   
c
       m=0
      do 22 ir=1,mt
         do 22 ic=ir,mt
            m=m+1
            z = v(m)*uu(ir)*uu(ic)
            write(14,1402)ir,ic,z,v(m)
 1402 format(1x,2i3,1x,2d20.10)       
 22   continue
      close(14)
c
c  read in PE matrix
c
      mt = 8
      v = 0.0
      open(10,file='vcvPE.d',form='formatted',status='old')
      write(17,8002)
 8002 format(/'  PERMANENT ENVIRONMENTAL VCV')
c
      tt = 0.0
 30   read(10,1001,end=40)ir,ic,x
        ir = ir - 8
        ic = ic - 8
      if(ir.gt.mt)go to 30
      if(ic.gt.mt)go to 30
      if(ir.lt.1)go to 40
      if(ir.eq.ic)then
         tt(ir)=dsqrt(x)
         m = ihmssf(ir,ic,mt)
         v(m) = 1.0
       else
          m=ihmssf(ir,ic,mt)
          v(m)=x
       endif   
      write(17,3002)ir,ic,x
      go to 30
c
 40   m=0
      do 31 ir=1,mt
         do 33 ic=ir,mt
            m=m+1
            v(m)=v(m)*tt(ir)*tt(ic)
 33      continue
 31   continue
c   
      call PDcheq(v,mt)
      close(10)
c
      open(14,file='vcvPE2.d',form='formatted',status='unknown')
      m = 0
      do 42 ir=1,mt
         do 42 ic=ir,mt
            m=m+1
            z = v(m)*uu(ir)*uu(ic)
            write(14,1402)ir,ic,z,v(m)
 42   continue
      close(14)
c
c  Genetic and RYSH matrices (correlations)
c
      open(11,file='vcvGH.d',form='formatted',status='old')
      write(17,8003)
 8003 format(/'  GENETIC and RYSH vcv ')
      mt = 16
      v = 0.0
      u = 0.0
      ss = 0.0
      tt = 0.0
 50   read(11,1101,end=60)ir,ic,x,z
 1101 format(2i3,1x,f8.0,1x,f12.0)
      if(ir.gt.mt)go to 50
      if(ic.gt.mt)go to 50
      if(ir.lt.1)go to 60
      if(ir.eq.ic)then
         ss(ir)=dsqrt(x)
         tt(ir)=dsqrt(z)
         m=ihmssf(ir,ic,mt)
         v(m)=1.0
         u(m)=1.0
      else
         m=ihmssf(ir,ic,mt)
         v(m) = x
         u(m) = z
      endif
      write(17,3007)ir,ic,x,z
 3007 format(1x,2i4,1x,2d20.10)
      go to 50
c
c  change to covariance matrices
c
 60   close(11)
      m = 0
      do 61 ir = 1,mt
         do 62 ic = ir,mt
            m = m + 1
            v(m) = v(m)*ss(ir)*ss(ic)
            u(m) = u(m)*tt(ir)*tt(ic)
 62      continue
 61   continue
      call PDcheq(v,mt)
      call PDcheq(u,mt)
      open(14,file='vcvGH2.d',form='formatted',status='unknown')
      m = 0
      do 64 ir=1,mt
         do 64 ic=ir,mt
            m=m+1
            b = v(m)*uu(ir)*uu(ic)
            z = u(m)*uu(ir)*uu(ic)
            write(14,1402)ir,ic,b,z
 64   continue
      close(14)

      stop
      end


      include 'PDsub.f'
      include 'dkmvhf.f'