**********************************************************************
* REML variance component estimation using sparse matrix techniques.
* Supports a class of mixed models with fixed and random factors,
* including animal models.
*     Version for sparse package FSPAK
*
* This version supports regressions (identified by number of levels <0).
* The coefficient matrix need not be of full rank. 
*
* @ Ignacy Misztal, 6/29/1989 -- 7/20/1992 - 12/15/92 - 12/13/93 - 03/07/94
*
*
* Description of parameters
*
* nf     - maximum number of factors
* 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
*
* 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,optim
      parameter (nhash=50000,mx=5000,m1=.7*nhash,optim=1)
      parameter (nhashr=.8*nhash,m1r=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
*----------------------------------------------------------------

*--------------storage for fspak---------------------------------
      real*8  smprsp(1),tmp(mx)
      integer smpisp(1),smpesp,fratio

c_______________ storage for variance components________________
      integer ireml,ibeg,iend,rank
      real*8 red(nf),sigma(nf),denom(nf),tr(nf)
      real*8 sigmae,tracediag,tracematrix
      real*8 bb,uau,redtot,dentot,yy
      real*8 vex(nf),vprev(nf),valfa

c--------------- general storage--------------------------------
      real*8 var(nf),varold(nf),val(nf),buf(50),sol(mx),rs(mx),rel
      real second
      integer ifact(nf),ilev(nf),ilev1(nf),ipos(nf)
      integer nunin,nunout,nooffa,itr,imax,iposmx,
     $ i,j,k,j1,j2,k1,jj,iflag,n,nr,igroup,mg,isize
      character*15 npar,ndat,nrel,nout
      logical ireg(nf),relat
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 the parameters can be read either from the console or from a file
      write(*,*)'name of a file with parameters, con - if console:'
      read(*,3)npar
3     format(a15)
      if (npar.ne.'con') then
          nunin=1
          nunout=20
          open(nunin,file=npar)
          open(nunout,file='null')
        else
          nunin=5
          nunout=6
      endif
      write(nunout,*)'data name:'
      read(nunin,3)ndat
      write(nunout,*)'name of file with animal relationships:'
      read(nunin,3)nrel
      write(nunout,*)'name of output file:'
      read(nunin,3)nout
      open(2,file=ndat)
      open(3,file=nrel)
      open(7,file=nout)
5     write(nunout,*) 'no of factors, trait position, relax.factor ?'
      read(nunin,*)nooffa,itr,rel
      write(nunout,*)' for each factor give position, no of levels,'
      write(nunout,*)' and variance ratio, 0 if fixed, negative',
     * ' if relat. matrix'
      write(nunout,*)
      iposmx=itr
      isize=0
      relat=.false.
      do 10 i=1,nooffa
        read(nunin,*)ipos(i),ilev(i),var(i)
c temporarily, number of levels <0  denotes regression
        if (ilev(i).lt.0) then
              ilev(i)=1
              ireg(i)=.true.
           else
              ireg(i)=.false.
       endif
        varold(i)=var(i)
        if (var(i).lt.0) then
            if (relat) then
                print*,'only one relat. matrix allowed, type all again'
                stop
              else
              relat=.true.
            endif
        endif
        iposmx=max(iposmx,ipos(i))
        ilev1(i)=isize
10      isize=isize+ilev(i)
*
c some parameter validation
      if (isize.gt.mx) then
         print*,'too large model, increase parameter mx'
         stop
      endif
      if (nooffa.le.0) then
          print*,' number of effects =< 0'
          stop
      endif
      if (nooffa.gt.nf) then
          print*,'number of effects: ',nooffa,' too large,  increase NF'
          stop
      endif
      if (optim.eq.0) then
          print*,'Optimizing for speed'
        else
          print*,'Optimizing for memory'
      endif
*
c initialize
26    do 20 i=1,mx
        ria(i)=0
        rs(i)=0
20      sol(i)=0
      do 25 i=1,3
        do 25 j=1,nhash
25         x(j,i)=0
      do 28 i=1,3
        do 28 j=1,nhashr
28        rx(j,i)=0
      imax=1
      n=0
      ix=0
      ixr=0
      yy=0
      sigmae=0
      mg=0
      ireml=1
      tmp(1)=second()
*
c   ******** while not eof(2) do *************
c
30    read(2,*,end=65)(buf(i),i=1,iposmx)
        n=n+1
        do 50 i=1,nooffa
          ifact(i)=buf(ipos(i))
          if (ireg(i)) then
                 val(i)=ifact(i)
                 ifact(i)=ilev1(i)+1
             else
               if (ifact(i).le.ilev(i)) then
                     val(i)=1
                  else
                     val(i)=0
               endif
               ifact(i)=ifact(i)+ilev1(i)
          endif
50      continue
        yy=yy+buf(itr)**2
        do 60 j=1,nooffa
          j1=ifact(j)
          rs(j1)=rs(j1)+buf(itr)*val(j)
          do 60 k=1,nooffa
            k1=ifact(k)
60          if (j1.le.k1 .and. val(j)*val(k).ne.0)
     +               call hashmx(val(j)*val(k),j1,k1,x,nhash,ix)
        goto 30

c ******** end while ******

* add variance ratio to diagonals of random factors
65    do 100 jj=1,nooffa
        j1=ilev1(jj)+1
        j2=ilev(jj)+ilev1(jj)
        if (var(jj).gt.0) then
           do 70 k=j1,j2
70           call hashmx(var(jj),k,k,x,nhash,ix)
        endif

       if (var(jj).lt.0) then
           nr=0
           igroup=0
* add relationship matrix if such exists
75         read(3,*,end=95)(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(jj)
           call relhash(ifact,x,nhash,ix,-var(jj),storhalf)
           goto 75
95         igroup=ilev(jj)-nr
           print 97,nr,igroup
        endif
97    format(' read',i7,' pedigree records',i8,' unknown parent groups')

100   continue

      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
      mg=min(1,igroup)
c     matrix is assumed to be of full rank
      isize=isize-mg

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

c      if (relat) call prsparse(ra,ria,rja,ilev(nooffa)-mg)
c      call prsparse(smpa,smpia,smpja,isize)


c verify that matrices have positive diagonals
c      call checkdiag(isize,'coefficient',smpia,smpja,smpa)
c      call checkdiag(ilev(nooffa)-mg,'relationship',ria,rja,ra)
c      call prsparse(ra,ria,rja,ilev(nooffa)-mg)
c      call prsparse(smpa,smpia,smpja,isize)
c open temporary file for storing the coefficient matrix
      open(10,form='unformatted')
      nsmp=smpia(isize+1)-1
      write(10)isize,nsmp,(smpia(i),i=1,isize+1),(smpa(i),i=1,nsmp),
     $          (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,isize,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,isize,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i,rank)
            call fspakerr(iflag)
            print*,'ordering time =',second()
      endif

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

      print*,'symbolic factorization time=',second()
      print '('' used'',f7.2,''% memory'')',smpesp*100./maxss
      call fspakerr(iflag)

c numerical factorization

190   call fspak(40,isize,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +               smpesp,smpisp,j,j,j,j,j,j,rank)
      if (ireml.eq.1) then
         print*,'numerical factorization time=',second()
         print '('' used'',f7.2,''% memory'')',smpesp*100./maxss
         print*,'rank=',rank
      endif
      call fspakerr(iflag)

c obtain solutions
      do  j=1,isize
         sol(j)=rs(j)
      enddo
      call fspak(50,isize,smpia,smpja,smpa,sol,iflag,6,30,maxss,
     +               smpesp,smpisp,j,j,j,j,j,j,rank)
      call fspakerr(iflag)

c Invert
      rewind 10
      read(10)i,j,(smpia(i),i=1,isize+1),(smpa(i),i=1,nsmp),(smpja(i),
     $          i=1,nsmp)
c check if sufficient memory for inverse in smpia-smpja-smpa
      i=isize
      j=smpia(i+1)
      k=(j-i)*2+i
      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,isize,smpia,smpja,smpa,rs,iflag,6,30,maxss,
     +            smpesp,smpisp,i,i,i,i,i,i,rank)
      if (ireml.eq.1) then
        print*,'sparse inversion time=',second()
        print '('' used'',f7.2,''% memory'')',smpesp*100./maxss
      endif
      call fspakerr(iflag)
* ========================REML ==============================
      redtot=0
      dentot=0
      do 200 jj=1,nooffa
        red(jj)=0
        tr(jj)=0
        if (var(jj).eq.0) then
            do 210 i=ilev1(jj)+1,ilev1(jj)+ilev(jj)
210            red(jj)=red(jj)+sol(i)*rs(i)
            redtot=redtot+red(jj)
            goto 200
        endif
        if (var(jj).lt.0) then
c tr(A-1W)
             ibeg=ilev1(jj)+1
             iend=ilev1(jj)+ilev(jj)-igroup
             do 220 i=ibeg,iend
220             red(jj)=red(jj)+sol(i)*rs(i)
             tr(jj)=tracematrix(isize,smpia,smpja,smpa,ilev(jj)-mg,
     $              ria,rja,ra,ilev1(jj),ilev1(jj),tmp)
             print*,'trace ',second()
c  u'Au
             call quadrf(uau,ra,ria,rja,sol,1,ilev(jj)-mg,ilev1(jj))
             red(jj)=red(jj)+uAu*(-var(jj))
c             print*,'quadratic form',second()
             denom(jj)=ilev(jj)-igroup+tr(jj)*var(jj)
             sigma(jj)=(uau)/denom(jj)
             dentot=dentot+tr(jj)*(-var(jj))
          else
            bb=0
            ibeg=ilev1(jj)+1
            iend=ilev1(jj)+ilev(jj)
            do 280 i=ibeg,iend
              red(jj)=red(jj)+sol(i)*rs(i)
280           bb=bb+sol(i)**2
290           tr(jj)=tracediag(isize,smpia,smpja,smpa,ibeg,iend)
            denom(jj)=ilev(jj)
            sigma(jj)=bb/(denom(jj)-tr(jj)*var(jj))
            red(jj)=red(jj)+bb*var(jj)
            dentot=dentot+tr(jj)*var(jj)
        endif

        redtot=redtot+red(jj)

200   continue

      sigmae=(yy-redtot)/(n-rank+dentot)

c calculate variance ratios, compute the exponential adjustment
c if possible
      do 300 i=1,nooffa
      if (var(i).ne.0) then
          valfa=sign(sigmae/sigma(i),var(i))
          if (vprev(i).eq.0 .or. ireml .lt.2) then
               vex(i)=0
             else
               vex(i)=(var(i)**2-valfa*vprev(i))/
     $                             (2*var(i)-valfa-vprev(i))
          endif
          vprev(i)=var(i)
          var(i)=valfa

          write(*,310)ireml,sigmae,i,sigma(i),abs(var(i)),vex(i)
          write(7,310)ireml,sigmae,i,sigma(i),abs(var(i)),vex(i)
310       format(' vround',i3,'  se=',f12.4,'  s(',i2,')=',f12.4,
     $           2x,'alfa=',f9.4,'  ex=',f9.4)
      endif
300   continue
      ireml=ireml+1
      if (ireml.gt.imax) then
         print*,'how many iterations more?, extrapolate (1)?'
         read*,j,j1
         imax=imax+j
      endif
      if (j1.eq.1) then
c extrapolate
          do 320 i=1,nooffa
             if (var(i).ne.0) then
                var(i)=vex(i)
                vprev(i)=0
             endif
320       continue
      endif
      if (ireml.le.imax) then
c  read the coefficient matrix, the previous one has been overwritten
c by SMPAK.
c      print*,'before reloading',second()
      rewind 10
      read(10)i,j,(smpia(i),i=1,isize+1),(smpa(i),i=1,nsmp),(smpja(i),
     $          i=1,nsmp)
c        print*,'after reloading',second()
        do 400 jj=1,nooffa
          j1=ilev1(jj)+1
          j2=ilev(jj)+ilev1(jj)
          if (var(jj).gt.0) then
             do 410 k=j1,j2
410           call addspar(smpa,smpia,smpja,k,k,var(jj)-varold(jj))
          endif
         if (var(jj).lt.0) then
* add relationship matrix if such exists
            call adjadd(smpa,smpia,smpja,j1,j2-mg,ra,ria,rja,1,
     $                  varold(jj)-var(jj),storhalf,storfull)
         endif
400      continue
         goto 190
      endif
*
      close(10,status='delete')
* write down solutions
      write(7,*)
      write(7,*) n,' records',isize,' equations'
      do 550 i=1,nooffa
        write(7,*)
        write(7,*) 'solutions of factor',ipos(i),'  var. ratio',var(i)
550     if (ilev(i).lt.100) write(7,570)(j,sol(j+ilev1(i)),j=1,ilev(i))
570   format(4(3x,i4,':',f11.4))
      open(21,file='sol'//npar)
      do 580 i=1,isize
580     write(21,'(f11.4)')sol(i)
      print*,'solutions were stored in file sol'//npar
      end