c-----------------------------------------------------------------------      
c Program JAA version 2.0
c
c Computes solutions of mixed model equations directly from the data.    
c Supports animal models. Contains fast I/O subroutines.
c Version 2.0 can calculate repeatabilities of animals
c under the animal model for a model containing HYS (or equivalent),
c permanent environment (optional) and animal effects.
c
c Design of the program is based on the following papers:
c Iteration on data:
c   Misztal, I., and D. Gianola.  1987.  Indirect solution of
c       mixed model equations.  J. Dairy Sci. 70:716-723.
c   Misztal, I., D. Gianola, and J. L. Foulley.  1989.
c       Computing aspects of a nonlinear method of sire evaluation
c       for categorical data.  J. Dairy Sci. 72:1557-1568.
c Repeatability:
c   Misztal, I., and G.R. Wiggans. 1988. Approximation of
c       prediction error variance in large-scale animal models. J.
c       Dairy Sci. 71 (Suppl. 2):27-32.
c   Misztal, I., G. R. Wiggans, T.J. Lawlor, and T. Short. 1991.
c       Continuous genetic evaluation of dairy cattle. J. Dairy Sci.
c       74:2001-2009.
c
c  Copyright Ignacy Misztal & University of Illinois
c  Dec 31, 1988- Oct 2, 1992
c
c  General use of this program is allowed if acknowledged in publications.
c
c  Parameters
c     mx - total number of equations
c
c  Vectors
c     acc(mx)  - adjusted right-hand sides
c     diag(mx) - diagonals of the system
c     sol(mx)  - solutions
c     dsol(mx) - old solutions
c     elm(mx)  - effective number of records
c     pevn(mx) - current information on animals
c     pevo(mx) - previous round information on animal
c      
c-----------------------------------------------------------------------      
c
      parameter (mx=20000)                                               
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      real  acc(mx),diag(mx),sol(mx),dsol(mx) 
      real var(20),buf(50),val(20)
      integer ifact(10),ifact1(20),ilev(20),ilev1(20),ipos(20),ifil(9)
      logical ireg(20)
      character*15 npar,ndat,nrel,nout

*--------------storage for pev-----------------------------------
      real pevo(mx),pevn(mx),itot(3),ifra(3),elm(mx)
      integer icut,npev,ipev
*----------------------------------------------------------------
      equivalence (acc(1),pevo(1)),(dsol(1),pevn(1))

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 pos., rel. factor,'
     * ,'iterations?:'  
      read(nunin,*)nooffa,itr,rel,imax                                      
      write(nunout,*)' for each factor give position, no of levels,'   
      write(nunout,*)' and variance ratio, 0 if fixed, negative',
     * ' if relat. matrix'
      write(nunout,*)'  no of levels < 0 denotes covariable'
      write(nunout,*)
      iposmx=itr                                                        
      isize=0                                                           
      irelat=0                                                          
      do 10 i=1,nooffa                                                  
        read(nunin,*)ipos(i),ilev(i),var(i)                                   
c number of levels <0  denotes covariable
        if (ilev(i).lt.0) then
              ilev(i)=1
              ireg(i)=.true.
           else
              ireg(i)=.false.
        endif
        if (var(i).lt.0) then                                             
          if (irelat.eq.1) then                                        
            print*,'only one relat. matrix allowed, type all again' 
            stop                                                            
         else                                                              
          irelat=1                                                          
        endif                                                             
        endif
        iposmx=max(iposmx,ipos(i))                                        
        ilev1(i)=isize                                                    
10      isize=isize+ilev(i)                                               
*                                                                       
      if (isize.gt.mx) then
         print*,'too large model, increase parameter mx'
         stop
      endif
      write(*,*)'iteration no.  convergence criterion   last solution'
      write(7,*)'iteration no.  convergence criterion   last solution'
* rewrite the data to an unformatted file to gain on speed              
      call iof(iopen,ifil,n,4)
      call iof1(iopen,ifil,n,8)
      n=0                                                               
22    read(2,*,end=24)(buf(i),i=1,iposmx)                               
      n=n+1                                                             
      do 23 i=1,nooffa                                                  
23      ifil(i)= int(buf(ipos(i)))                                      
      ifil(nooffa+1)=buf(itr)                                         
      n1=nooffa+1
      call iof(iwrite,ifil,n1,4)     
      goto 22                                                           
*                                                                       
24    nr=1
      if (irelat.eq.1) then
27      read(3,*,end=26)(ifact(i),i=1,4)
        nr=nr+1
        call iof1(iwrite,ifact,4,8)
        goto 27
      endif
26    do 20 i=1,mx                                                   
        dsol(i)=0
20      sol(i)=0                                                          
      i=1                                                               
*                                                                       
25    adj=0
      adj0=0
      do 30 j=1,isize
        acc(j)=0
30      diag(j)=0

* do one round of gauss-seidel (jacobi if relationship matrix)
      do 110 jj=1,nooffa
        call iof(irew,ifact,0,4)
        do 60 ii=1,n
          n1=nooffa+1
          call iof(iread,ifact1,n1,4)
          y=float(ifact1(n1))
          do 50 j=1,nooffa
           if (ireg(j)) then
                val(j)=ifact1(j)
                ifact1(j)=1
            else
              if (ifact1(j).le.ilev(j) .and. ifact1(j).ge.1) then
                  val(j)=1
                 else
                  val(j)=0
              endif
           endif
           ifact(j)=ifact1(j)+ilev1(j)
50         continue
          j1=ifact(jj)
          diag(ifact(jj))=diag(ifact(jj))+val(jj)*val(jj)
          acc(ifact(jj))=acc(ifact(jj))-y*val(jj)
        do 60 k=1,nooffa
            k1=ifact(k)
60        if (k1 .ne. j1) acc(ifact(jj))=acc(ifact(jj))+
     +                         sol(k1)*val(jj)*val(k)

* solve each factor separately, only fixed factors during the first
* round
40    j1=ilev1(jj)+1
      j2=ilev(jj)+ilev1(jj)

* add variance ratio to diagonals of random factors
       if (var(jj).gt.0) then
           do 70 k=j1,j2
70           diag(k)=diag(k)+var(jj)
       endif

       if (var(jj).ge.0) goto 80                                 
* add animal relationship matrix if such exists
         call iof1(irew,ifact,0,8)
         do 75 k=1,nr-1
           call iof1(iread,ifact,4,8)
           ian=ifact(1)
           isr=ifact(2)
           idam=ifact(3)
           itype=ifact(4)
           if (itype.eq.1) type=-var(jj)
           if (itype.eq.2) type=-var(jj)*2./3
           if (itype.eq.3) type=-var(jj)*.5
           ian=ian+ilev1(jj)
           isr=isr+ilev1(jj)
           idam=idam+ilev1(jj)
           diag(ian)=diag(ian)+2*type                                
           diag(idam)=diag(idam)+type*0.5                                   
           diag(isr)=diag(isr)+type*0.5                                
           acc(ian)=acc(ian)-type*(sol(isr)+sol(idam))          
           acc(isr)=acc(isr)-type*(sol(ian)-.5*sol(idam))          
75         acc(idam)=acc(idam)-type*(sol(ian)-.5*sol(isr))         

* solve
80    do 110 k=j1,j2    
      rel1=0
        if (var(jj).lt.0) rel1=rel 
        if (diag(k).eq.0) goto 110                                        
        add=-acc(k)/diag(k)-sol(k)                                        
        adj=adj+add*add        
        xx=sol(k)                                               
        sol(k)=sol(k)+rel1*(sol(k)-dsol(k))+add                           
        dsol(k)=xx
105     adj0=adj0+sol(k)*sol(k)                                                                                 
110     continue
*                                                                       
*                                                                       
130   adj=adj/adj0                                                      
      write(7,*)i,adj,sol(isize)
      write(*,*)i,adj,sol(isize)                                        
      i=i+1                                                             
      if (i.lt.imax ) goto 25                                           
      write(*,135)                                                      
135   format(' how many more iterations, 0-none?')                       
      read(*,*)j                                                        
      if (j.gt.0) then                                                  
      imax=imax+j                                                       
      goto 25                                                           
      endif                                                             
*
* +++++++++++++   p e v begins here ++++++++++++++++++
      if (var(nooffa).ge.0) goto 500
      print 200
200   format(' In approximate Prediction Error Variance for animals',
     +       ' it will be assumed that:'/
     +        '  first effect is fixed like HYS'/
     +        '  the last effect is animal')
      if (var(nooffa-1).ne.0) print*,'  effect before the last is',
     +     ' permanent environment'
      write(nunout,'('' calculate approximate PEV (1=yes)?'')')
      read(nunin,*,end=500)ipev
      if (ipev.ne.1) goto 500
      icut=ilev(nooffa)-nr+1
      print '('' detected'',i8,'' unknown parent group(s)'')',icut
      write(*,'(/10x,'' Calculation of approximate PEV''/)')
      write(7,'(/10x,'' Calculation of approximate PEV''/)')
      write(*,'('' iteration no. convergence PEV(1)    PEV(2) ....'')')
      write(7,'('' iteration no. convergence PEV(1)    PEV(2) ....'')')
c initialize
      npmax=1
      npev=1
      var1=var(nooffa-1)
      var2=abs(var(nooffa))
      do 350 i=1,isize
        elm(i)=0
        pevn(i)=0
350     pevo(i)=0


c calculate and add effective lactations to new pev
c animal effects is assumed last on the list (nooffa)
c
        call iof(irew,ifact,0,4)
      do 430 ii=1,n
        call iof(iread,ifact,n1,4)
        do 440 j=2,nooffa-1
440        ifact(j)=ifact(j)+ilev1(j)
           i1=ifact(nooffa)
           elm(i1)=elm(i1)+(1-1/diag(ifact(1)))
430   continue
c
c If no permanent environment effect present, comment lines below
c
      if (var(nooffa-1).ne.0.0) then
       do 460 i=1,ilev(nooffa)-icut
           elm(i)=elm(i)*var1/(elm(i)+var1)
460      continue
      endif
c end of permanent environment section

            do 465 i=1,isize
465           pevo(i)=elm(i)

c start pev loop here
467    do 468 i=1,isize
468       pevn(i)=elm(i)
c calculate contributions due to relationships
      call iof1(irew,ifact,0,8)
      do 470 ii=1,nr-1
           call iof1(iread,ifact,4,8)
           ian=ifact(1)
           isr=ifact(2)
           idam=ifact(3)
           if (idam.le.ilev(nooffa)-icut .and. idam.ne.0) then
                          itot(1)=pevo(idam)
             else
               itot(1)=0
           endif
           if (isr.le.ilev(nooffa)-icut .and. isr.ne.0) then
               itot(2)=pevo(isr)
             else
               itot(2)=0
           endif
           itot(3)=pevo(ian)
           call deduct(itot,ifra,-var(nooffa))
c          print '(3i3,9f10.6)',idam,isr,ian,itot,ifra
           if (idam.le.ilev(nooffa)-icut .and. idam.ne.0)
     $                             pevn(idam)=pevn(idam)+ifra(1)
           if (isr.le.ilev(nooffa)-icut .and. isr.ne.0)
     $                             pevn(isr)=pevn(isr)+ifra(2)
           pevn(ian)=pevn(ian)+ifra(3)
470   continue

      adj=0
      adj1=0
      npev=npev+1
      do 480 i=1,ilev(nooffa)-icut
        adj=adj+abs(pevn(i)-pevo(i))
c480     pevo(i)=pevn(i)
480     pevo(i)=(pevn(i)+pevo(i))/2
      write(*,'(i13,8f10.4)')npev-1,adj,(pevn(i),i=1,5)
      write(7,'(i13,8f10.4)')npev-1,adj,(pevn(i),i=1,5)
      if (npev.le.npmax) goto 467
      print*,'more pev rounds?'
      read*,jj
      if (jj.ne.0) then
          npmax=npmax+jj
          goto 467
      endif

* write down solutions                                                  
500   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')
      if (ipev.eq.1) then
          do i=1,isize
            if (i.gt.ilev1(nooffa)) then
                xx=pevn(i-ilev1(nooffa))
              else
                 xx=0
            endif
            write(21,'(f11.4,f6.3)')sol(i),xx/(xx-var(nooffa))
          enddo
          print*,'solutions and reliabilities were stored in file sol'
        else
          do i=1,isize
            write(21,'(f11.4)')sol(i)
          enddo
          print*,'solutions were stored in file sol'
      endif
      end


      subroutine deduct(b,g,al)
c from total information on animal and its parents, calculates fraction
c due to this particular relationships.
c  b - on input contains total information for parent 1, parent 2, and animal
c  g - on output contains fraction f the total information due to this relationship
c  al - variance ratio
c
c unknown parent has information = 0.
*
      real b(3),g(3),al,adj,q(3),s(3),bg(3)
      real al2d,al2p,al2p4d,al3p,al3m2d,al2m
      integer i,j,n
*
      n=1
      al2p=al*al
      al2p4d=al2p/4
      al3p=al2p*al
      al3m2d=3*al/2
      al2m=2*al
      al2d=al/2
*
      do 10 i=1,3
10       q(i)=b(i)
*
15    bg(1)=al2d+q(1)-(al2p4d*(al2m+q(3))-al3p+al2p*(al3m2d+q(2)))/
     *  ((al3m2d+q(2))*(al2m+q(3))-al2p)

      bg(2)=al2d+q(2)-(al2p4d*(al2m+q(3))-al3p+al2p*(al3m2d+q(1)))/
     *  ((al3m2d+q(1))*(al2m+q(3))-al2p)

      bg(3)=al+q(3)-(al2p*(al3m2d+q(2))-al3p+al2p*(al3m2d+q(1)))/
     *  ((al3m2d+q(1))*(al3m2d+q(2))-al2p4d)
*
      adj=0
      do 20 j=1,3
        if (b(j).eq.0) then
            s(j)=0
          else
            s(j)=bg(j)-b(j)
        endif
        adj=adj+s(j)*s(j)
        q(j)=q(j)-s(j)
        if (q(j).lt.0) q(j)=0
        g(j)=bg(j)-q(j)
20    continue
*
      n=n+1
      if (n.lt.10 .and. adj.gt.10e-10) goto 15
*
      end



      subroutine iof(itype,vector,n,iun)
c I/O facility to read and write vectors of any size using
c large buffer, thus increasing unformatted i/o speed.
c
c itype  - operation type, as in the parameter statement below,
c vector - contains data to be stored or retrieved,
c      n - size of vector,
c    iun - unit number.
c 
c If the data is small enough to fit in the buffer, subsequent reads 
c are done entirely in memory with a RAM speed. Option "iopen" creates
c a temporary file. Definition of a named file should be done explicitly
c in a user program, followed by a call to this subroutine with "irew"
c option.
c
c For greater speed the subroutine should be used with one file only.
c More copies of it with different names (iof1,iof2...) can be used if
c more files are desired.
c
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      implicit integer (b,f,i)
      parameter(bufs=1024)
      integer x(bufs),vector(n)
      save x,fpos,bpos,iwr,lastit
      data fpos,bpos/2*-1/
***********************************************************************      
      if(itype.eq.iread) then
        iwr=0
        if (fpos.eq.0) then
          read(iun,end=50)x
          fpos=1  
          bpos=0
        endif
        n1=min(n,bufs-bpos)  
        do 30 i=1,n1
30        vector(i)=x(bpos+i)
        bpos=bpos+n1
        if (bpos.eq.bufs .and. n1.ne.n) then
          read(iun,end=50)x
          fpos=fpos+1
          bpos=n-n1
          do 40 i=1,bpos
40          vector(i+n1)=x(i)
          endif
* 
      else if(itype.eq.iwrite) then
        if (lastit.eq.irew) then
          fpos=0
          rewind iun
        endif  
        iwr=1
        n1=min(n,bufs-bpos)
        do 10 i=1,n1
10      x(bpos+i)=vector(i)
        bpos=bpos+n1
        if (bpos.eq.bufs) then
          write(iun)x
          fpos=fpos+1
          bpos=n-n1
          do 20 i=1,bpos
20          x(i)=vector(i+n1)                 
        endif
* 
      elseif (itype.eq.iopen) then
        bpos=0
        fpos=0
        open(iun,access='sequential',form='unformatted')
* 
      else if (itype.eq.irew) then
        if (bpos.ne.0 .and. iwr.eq.1) then
          write(iun)x
          fpos=fpos+1
        endif
        iwr=0
        bpos=0
        if (fpos.ne.1) then
          fpos=0
          rewind iun
        endif
* 
      else
        write(*,*)'unknown i/o code=',itype
        stop
      endif
*********************************************************************** 
      lastit=itype
      return
50    write(*,*)'end of file',iun
      write(*,*)'bpos,fpos,n',bpos,fpos,n
      stop
      end

      subroutine iof1(itype,vector,n,iun)
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      implicit integer (b,f,i)
      parameter(bufs=1024)
      integer x(bufs),vector(n)
      save x,fpos,bpos,iwr,lastit
      data fpos,bpos/2*-1/
***********************************************************************      
      if(itype.eq.iread) then
        iwr=0
        if (fpos.eq.0) then
          read(iun,end=50)x
          fpos=1  
          bpos=0
        endif
        n1=min(n,bufs-bpos)  
        do 30 i=1,n1
30        vector(i)=x(bpos+i)
        bpos=bpos+n1
        if (bpos.eq.bufs .and. n1.ne.n) then
          read(iun,end=50)x
          fpos=fpos+1
          bpos=n-n1
          do 40 i=1,bpos
40          vector(i+n1)=x(i)
          endif
* 
      else if(itype.eq.iwrite) then
        if (lastit.eq.irew) then
          fpos=0
          rewind iun
        endif  
        iwr=1
        n1=min(n,bufs-bpos)
        do 10 i=1,n1
10      x(bpos+i)=vector(i)
        bpos=bpos+n1
        if (bpos.eq.bufs) then
          write(iun)x
          fpos=fpos+1
          bpos=n-n1
          do 20 i=1,bpos
20          x(i)=vector(i+n1)                 
        endif
* 
      elseif (itype.eq.iopen) then
        bpos=0
        fpos=0
        open(iun,access='sequential',form='unformatted')
* 
      else if (itype.eq.irew) then
        if (bpos.ne.0 .and. iwr.eq.1) then
          write(iun)x
          fpos=fpos+1
        endif
        iwr=0
        bpos=0
        if (fpos.ne.1) then
          fpos=0
          rewind iun
        endif
* 
      else
        write(*,*)'unknown i/o code=',itype
        stop
      endif
*********************************************************************** 
      lastit=itype
      return
50    write(*,*)'end of file',iun
      write(*,*)'bpos,fpos,n',bpos,fpos,n
      stop
      end