c-----------------------------------------------------------------------      
c Program JAA
c
c Computes solutions of mixed model equations directly from the data.    
c Supports animal models. Contains fast I/O subroutines.
c
c  Copyright Ignacy Misztal & University of Illinois
c  Dec 31, 1988
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      
c-----------------------------------------------------------------------      
c
      parameter (mx=20000)                                               
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      real  acc(mx),diag(mx),sol(mx),dsol(mx) 
      dimension aver(20),var(20),buf(50)                                
      integer ifact(10),ilev(20),ilev1(20),ipos(20),ifil(9)      
      character*15 npar,ndat,nrel,nout,ncon

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,*)
      iposmx=itr                                                        
      isize=0                                                           
      irelat=0                                                          
      do 10 i=1,nooffa                                                  
        read(nunin,*)ipos(i),ilev(i),var(i)                                   
        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
27    read(3,*,end=26)(ifact(i),i=1,4)
      nr=nr+1  
      call iof1(iwrite,ifact,4,8)
      goto 27
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,ifact,n1,4)
          y=float(ifact(n1))                                        
          do 50 j=2,nooffa                                                  
50          ifact(j)=ifact(j)+ilev1(j)                                        
          j1=ifact(jj)                                                      
          diag(j1)=diag(j1)+1                                               
          acc(j1)=acc(j1)-y                                                 
          do 60 k=1,nooffa                                                  
            k1=ifact(k)                                                       
60          if (k1 .ne. j1) acc(j1)=acc(j1)+sol(k1)                           
                                                                        
* solve each factor separately, only fixed factors during the first     
* round                                                                 
40    aver(j)=0                                                         
      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)                                           
        aver(j)=aver(j)+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                                                             
*                                                                       
* write down solutions                                                  
      write(7,*)     
      write(7,*) n,' records',isize,' equations'                        
      do 150 i=1,nooffa                                                 
        write(7,*)                                                        
        write(7,*) 'solutions of factor',ipos(i),'  var. ratio',var(i)    
150     if (ilev(i).lt.100) write(7,170)(j,sol(j+ilev1(i)),j=1,ilev(i))
170   format(4(3x,i4,':',f11.4))                                        
      open(21,file='sol')
      do 180 i=1,isize
180     write(21,'(f11.4)')sol(i)                                    
      print*,'solutions were stored in file sol'
      stop                                                              
      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