********************************************************************
c Calculates variance components for multiple trait problems.
c It uses the canonical transformation. Consequently, only models-with 
c single random factor(including animal with groups) are supported.
c Computational costs are reduced by calculating the traces for a 
c range of values once and interpolating the other values later.
c Computationally intensive subroutines are optimized for Cray
c supercomputers. The program is built on another single-trait program.
c
c copyright (C) 1989 Ignacy Misztal 9-1-1989 -- 12/27/94.
c changed to double precision 4-13-92
c changed to true sparse inversion 5-30-92
c changed to FSPARSE sparse package 12-14-92
c added low-memory inversion option 12-13-1993
*********************************************************************

* REML variance component estimation using sparse matrix techniques.
* Supports a class of mixed models with fixed and random factors,
* including animal models. 
*
* This version supports regressions (identified by number of levels <0).
* The coefficient matrix must be of full rank. It can be made such by
* declaring the number of levels for appropriate fixed factors as smaller by
* one. Then, the last levels of these factors are automatically assumed 0.
* 
* Description of parameters
*
* nf     - maximum number of factors
* ntrait - maximum number of traits
* mx     - maximum number of equations
* nhash  - 1.2*maximum number of non-zero upper-diagonal entries in the
*          coefficient matrix.
* nhashr - 1.2*maximum number of non-zero entries of the animal
*          relationship matrix, including Westell's group effects
* ntab   - maximum number of tabulated trace points
*
* optim  - 0 for speed optimization and 1 for memory optimization
*
* After the hash storage is converted to SMPAK form, the first one is 
* reused for the sparse matrix solver
**********************************************************************
      integer nhash,mx,maxss,m1,nhashr,m1r,nf,ntrait,ntab,optim
      parameter (nhash=50000,mx=2600,m1=.9*nhash,ntrait=18,
     +           ntab=80,optim=1)
      parameter (nhashr=.7*nhash,m1r=.9*nhashr)
      integer storfull,storhalf
      parameter (storfull=1,storhalf=2)
      parameter (nf=5)
 
*--------------storage for the coefficient matrix----------------
      real*8 x(nhash,3),smpa(m1)
      integer smpia(mx+1),smpja(m1),ix,nsmp
*----------------------------------------------------------------
 
*--------------storage for the relationship matrix---------------
      real*8 rx(nhashr,3),ra(m1r)
      integer ria(mx+1),rja(m1r),ixr,inorec
*----------------------------------------------------------------
      
*--------------storage for fspak---------------------------------
      real*8 smprsp(1),tmp(mx)
      integer smpisp(1),smpesp,fratio
 
c_______________ storage for variance components________________
      integer icanon,irank
      real*8 w(ntrait,ntrait),winv(ntrait,ntrait),red,
     +     uiAuj(ntrait,ntrait),redij(ntrait,ntrait),
     +     G(ntrait,ntrait),R(ntrait,ntrait),V(ntrait,ntrait),
     +     scr1(ntrait,ntrait),scr2(ntrait,ntrait),scr3(ntrait,ntrait),
     +     uau,yy(ntrait,ntrait),y11y(ntrait),tracematrix
      real*8 trtab(ntab,3),tr(ntrait),var,varold,varw(ntrait),
     +     convr,denom(ntrait),var0(ntrait),step,from,to,
     + convg,sumr,sumg,sigmae(ntrait,ntrait),sigma(ntrait,ntrait) 
      integer ipnts
c_______________________________________________________________
 
c---------------storage for extrapolated traces-----------------
      real*8 adj,adjrel
c---------------------------------------------------------------
       
c--------------- general storage--------------------------------
      real*8 sol(mx,ntrait),rs(mx,ntrait),relsor
      real second
      integer ifact(nf),ilev(nf),ilev1(nf),ipos(nf),itrno
      integer nooffa,itr,iposmx,nfix,nrand,itmax,
     $ i,j,k,j1,j2,l,iflag,n,nr,igroup,mg,isize
      logical ireg(nf)
      character*15 npar
c---------------------------------------------------------------
c storage for hash tables is reused by the sparse matrix solver
      COMMON /A/X,RX
      common /b/smpa,smpia,smpja
      equivalence (smprsp(1),smpisp(1),x(1,1))

c amount of storage for fspak, reused from hash tables
      maxss=3*fratio()*(nhash+nhashr)
*
c input parameters
      call parinp(nooffa,itr,itrno,iposmx,nf,ipos,ilev,ilev1,ireg,
     +            isize,npar,mx,ntrait,optim)
c first round relaxation factor=0
      relsor=0
c 
c initialize 
26    do 20 i=1,mx
20        ria(i)=0
      do 28 i=1,3
        do 28 j=1,nhashr
28        rx(j,i)=0          
      do 29 i=1,ntrait
         y11y(i)=0
         do 29 j=1,ntrait
            yy(i,j)=0
29          g(i,j)=0
      n=0
      itmax=1
      ix=0
      ixr=0
      tmp(1)=second()

c create the coefficient matrix in hash table + V
      call CoefV(n,nf,iposmx,nooffa,ilev,ilev1,ireg,ipos,ntrait,yy,
     + y11y,V,itrno,itr,x,nhash,ix)
 
* add relationship matrix
      varold=1.0
      nr=0
75       read(3,*,end=100)(ifact(i),i=1,4)
         nr=nr+1
         call relhash(ifact,rx,nhashr,ixr,1.0d0,storfull)
         do 80 i=1,3
80         ifact(i)=ifact(i)+ilev1(nooffa)            
         call relhash(ifact,x,nhash,ix,varold,storhalf)
         goto 75
*         
100   continue
      igroup=ilev(nooffa)-nr
      print 102,nr,igroup
102   format(' read',i7,' pedigree records',i8,' unknown parent groups')
      print*,'filled',ix,' out of',nhash,'  elements in hash table 1'
      print*,'filled',ixr,' out of',nhashr,'  elements in hash table 2'
      print*,second(),'s'

c number of last equations to zero out (one group if exists)
      mg=min(1,igroup)
      irank=isize-mg
      nrand=ilev(nooffa)-igroup
      nfix=irank-nrand

c     copy everything to smpak format
      call hashia(x,nhash,irank,smpia,smpja,smpa,m1,tmp)
      print*,'smpa used ',2*smpia(irank+1)-2-isize+mg,
     $        ' out of ',2*m1,' words'
      call hashia(rx,nhashr,ilev(nooffa)-mg,ria,rja,ra,m1r,tmp)
      print*,'ra used ',ria(ilev(nooffa)-mg+1),
     $            ' out of ',m1r,' words'
      print*,'before odrv',second()

c verify that matrices have positive diagonals
      call checkdiag(irank,'coefficient',smpia,smpja,smpa)
      call checkdiag(ilev(nooffa)-mg,'relationship',ria,rja,ra)
c      call prsparse(ra,ria,rja,ilev(nooffa)-mg)
c      call prsparse(smpa,smpia,smpja,irank)
c open temporary file for storing the coefficient matrix
      open(10,form='unformatted')      
      nsmp=smpia(isize+1-mg)-1
      write(10)(smpa(i),i=1,nsmp),(smpia(i),i=1,isize+1),(smpja(i),
     $          i=1,nsmp) 
c========================================================
c prepare for inversion by smpak
c========================================================
c find ordering if it has not been found before
      open(30,file='ord'//npar)
      call fspak(70,irank,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i)
c if different rank, compute ordering
      if (iflag.eq.0) then
            print*,'read ordering information from file ','ord'//npar
         else
            call fspak(10,irank,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i)
             call fspakerr(iflag)
             print*,'ordering time =',second()
      endif

c symbolic factorization
      call fspak(20,irank,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i)

      print*,'symbolic factorization time=',second()
      print '('' used'',f7.2,''% memory'')',smpesp*100./maxss
      call fspakerr(iflag)
c
c restore traces from disk, if exist
      open(29,file='tra'//npar)
      read(29,*,end=105)step,from,to,ipnts
      print 107,from,to,step                                              
107   format(' traces exist from ',f7.1,' to ',f7.1,' step ',f7.2,
     + '  restore? (1/yes, 0/no)')
      read*,j
      if (j.eq.1) then
         read(29,'(2g20.12)')((trtab(i,j),j=1,2),i=1,ipnts)
         goto 141
      endif
c tabulate trace from FROM to TO step STEP
105   print*,'step,from,to?'
      read*,step,from,to
      var=from
c  read the coefficient matrix, the previous one has been overwritten
c by SMPAK. Add the new relationship matrix
      ipnts=0
220      rewind 10
         read(10)(smpa(i),i=1,nsmp),(smpia(i),i=1,isize+1),(smpja(i),
     $          i=1,nsmp)
         j1=ilev1(nooffa)+1
         j2=ilev(nooffa)+ilev1(nooffa)
         call adjadd(smpa,smpia,smpja,j1,j2-mg,ra,ria,rja,1,
     $                  var-varold,storhalf,storfull)

c numerically factorize
      call fspak(40,irank,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i)
      if (ipnts.eq.0) then
        print*,'numerical factorization time=',second()
        print '('' used'',f7.2,''% memory'')',smpesp*100./maxss
      endif
      call fspakerr(iflag)

c Invert
      rewind 10
      read(10)(smpa(i),i=1,nsmp),(smpia(i),i=1,isize+1),(smpja(i),
     $          i=1,nsmp)
c check if sufficient memory for inverse in smpia-smpja-smpa
      j=smpia(irank+1)-1
      k=(j-irank)*2+irank
      if (k.gt.m1 .and. optim.eq.0) then
          print*,'insufficent space for the inverse in smpia-smpja-smpa'
          print*,' increase m1 from ',m1,' to ',k
          stop
      endif

      call fspak(60+optim,irank,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i)
      if (ipnts.eq.0) then
        print*,'sparse inversion time=',second()
        print '('' used'',f7.2,''% memory'')',smpesp*100./maxss
      endif
      call fspakerr(iflag)

      tr(1)=tracematrix(irank,smpia,smpja,smpa,ilev(nooffa)-mg,
     $          ria,rja,ra,ilev1(nooffa),ilev1(nooffa),tmp)
      if (ipnts.eq.0) print*,'trace ',second()

         print*,'var,tr',var,tr(1)
         trtab(ipnts+1,1)=var
         trtab(ipnts+1,2)=tr(1)
         ipnts=ipnts+1
         var=var*step
         if (ipnts.lt.ntab .and. var.lt.to) goto 220
       rewind 29
       write(29,*)step,from,to,ipnts
       write(29,'(2g20.12)')((trtab(i,j),j=1,2),i=1,ipnts)
c
c read g and r if exist, create otherwise
141    open(31,file='gr'//npar)
       read(31,*,end=145)((g(i,j),j=1,itrno),i=1,itrno)
       read(31,*,end=145)((r(i,j),j=1,itrno),i=1,itrno)
       print*,' initial G and R read from the parameters'' file'
       goto 152
c set initial G to .20% of diagonal V
145    do  i=1,itrno
         do  j=1,itrno
            g(i,j)=0
            r(i,j)=0
         enddo
         g(i,i)=.5
         r(i,i)=1
       enddo
c adjustment for EMC by Paul VanRaden
c calculate coefficients of the canonical transformation
152   icanon=0
      adjrel=.0
250   icanon=icanon+1
      call tcanon(ntrait,itrno,g,r,w,winv,scr1,scr2,scr3,varw)  
      print*,'variance ratios after transformation'
      print '(6g12.4)',(varw(i),i=1,itrno)
      if (icanon.eq.1) print*,'tcanon',second()
      write(7,'('' w'')')
      call prden(w,ntrait,itrno,7,'(6f13.5)')
c maximum varw is 10e5
      do 260 i=1,itrno
260       if (varw(i).gt.100000) varw(i)=1e5
c construct transformed right-hand sides for all traits + sums of squares
      call trrhs(n,isize,iposmx,nooffa,ilev,ilev1,ipos,ireg,
     + ntrait,itrno,itr,mx,rs,w,yy)
      if (icanon.eq.1) print*,'trhs',second()
c calculate variance components for all single traits
      do 350 i=1,itrno
         if (icanon.lt.15) then
            do 310 j=1,isize
310              sol(j,i)=0
         endif
c---set up the system for the correct variance ratio      
300      var=varw(i)
         rewind 10
         read(10)(smpa(j),j=1,nsmp),(smpia(j),j=1,isize+1),(smpja(j),
     $         j=1,nsmp)
         j1=ilev1(nooffa)+1
         j2=ilev(nooffa)+ilev1(nooffa)
         call adjadd(smpa,smpia,smpja,j1,j2-mg,ra,ria,rja,1,
     $                  var-varold,storhalf,storfull)

         call fspak(40,irank,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +               smpesp,smpisp,j,j,j,j,j,j)
         call fspakerr(iflag)
c obtain solutions
         do  j=1,irank
            tmp(j)=rs(j,i)
         enddo
         call fspak(50,irank,smpia,smpja,smpa,tmp,iflag,6,30,maxss,
     +               smpesp,smpisp,j,j,j,j,j,j)
         call fspakerr(iflag)
         do j=1,irank
            sol(j,i)=tmp(j)
         enddo
c---get uAu
         call quadruv (uau,ra,ria,rja,sol(1,i),sol(1,i),1,
     +                           ilev(nooffa)-mg,ilev1(nooffa))

c---get reduction
         red=0
         do 320 j=1,irank
320         red=red+sol(j,i)*rs(j,i)
c---get trace
        call trinp
     +      (trtab,ntab,ipnts,nrand,inorec,var,tr(i))
        var0(i)=var
        denom(i)=nrand-var*tr(i)
        adj=(nrand-var*tr(i))*var*tr(i)/nrand
c       print*,'denom,uau,adj,tr',denom(i),uau*var,adj,tr(i)
        adj=min(adjrel*denom(i),adjrel*uau*var,adj)
        sigma(i,i)=(uau-adj/var)/(denom(i)-adj)
        sigmae(i,i)=(yy(i,i)-red)/(n-nfix)
        varw(i)=1/sigma(i,i)
        if (varw(i).gt.500) then
          varw(i)=1e5
          sigma(i,i)=sigmae(i,i)/1e5
          var=varw(i)
        endif
           write(*,330)icanon,i,varw(i),sigma(i,i),sigmae(i,i)
     +                 ,second()
330     format(' ican=',i3,' trait=',i2,
     +  '  var=',g10.4,' sa=',g10.4,' se=',g10.4,'t=',f8.2)
350   continue
c transform back the R and G
      call quadrmt (uiauj,ra,ria,rja,sol,ntrait,itrno,mx,1,
     +                             ilev(nooffa)-mg,ilev1(nooffa))
      do 400 i=1,itrno
         do 400 j=1,itrno
            scr1(i,j)=0
            scr2(i,j)=0
            redij(i,j)=0
            do 405 k=1,isize-igroup
405            redij(i,j)=redij(i,j)+sol(k,i)*rs(k,j)
400     continue
c
c      print*,'after uiauj',second()
      do 408 i=1,itrno
         do 407 j=i+1,itrno
            sigma(i,j)=uiAuj(i,j)/sqrt(denom(i)*denom(j))
            sigma(j,i)=sigma(i,j)
            sigmae(i,j)=(yy(i,j)-redij(i,j)-sqrt(var0(i)*var0(j))
     +                     *uiauj(i,j))/n
            sigmae(j,i)=sigmae(i,j)
407      continue
408      continue
      do 410 i=1,itrno
         do 410 j=1,itrno
            do 410 k=1,itrno
               do 410 l=1,itrno
               scr1(l,i)=scr1(l,i)+winv(l,k)*sigma(k,j)*winv(i,j)            
410            scr2(l,i)=scr2(l,i)+winv(l,k)*sigmae(k,j)*winv(i,j)            
      convr=0
      convg=0
      sumr=0
      sumg=0
      do 420 i=1,itrno
         do 420 j=1,itrno
            convg=convg+(scr1(i,j)-g(i,j))**2
            convr=convr+(scr2(i,j)-r(i,j))**2
            sumg=sumg+g(i,j)**2
            sumr=sumr+r(i,j)**2
            g(i,j)=scr1(i,j)
420         r(i,j)=scr2(i,j)
      convg=convg/sumg
      convr=convr/sumr
      write(7,'(''sigma & sigmae'')')
      call prden(sigma,ntrait,itrno,7,'(8f13.5)')
      call prden(sigmae,ntrait,itrno,7,'(8f13.5)')
      print '('' convg='',g8.2,''      convr='',g8.2)',convg,convr
      write(7,*)'g,r, and heritabilities + correlations'
      call prden(g,ntrait,itrno,7,'(6f13.5)')
      call prden(r,ntrait,itrno,7,'(6f13.5)')
      rewind 31
      call prden(g,ntrait,itrno,31,'(6f13.5)')
      call prden(r,ntrait,itrno,31,'(6f13.5)')
      call corr(ntrait,itrno,g,r,scr1,7)
      
      write(7,*)'sigma and sigmae'
      call prden(sigma,ntrait,itrno,7,'(6f13.5)')
      call prden(sigmae,ntrait,itrno,7,'(6f13.5)')
c if not converged and below 50 rounds, repeat
c     if (icanon.lt.itmax .and. (convg .gt. 1e-9 .or. convr. gt. 1e-8)) 
      if (icanon.lt.itmax) 
     +                                                     goto 250
      print*,'finish? (rounds-more + adjrel)'
      read*,j,adjrel
      itmax=itmax+j
      if (icanon.lt.itmax .and. j.ne.0) goto 250
*
c calculate standard errors of the estimates
      call stderr(r,g,ntrait,itrno,trtab,ntab,ipnts,inorec,
     +                  n,nrand,nfix,7)
      print*,'heritabilities on diagonal, correlations g above,r below'
      call prden(scr1,ntrait,itrno,6,'(20f6.2)')
      call prden(scr1,ntrait,itrno,7,'(20f6.2)')
c
      close(10,status='delete')
      write(7,*)     
      write(7,*) n,' records',isize,' equations'                        
      end