c computes solutions of sire models with sire relationships
c either for linear or threshold model.
c matrix version, uses absorption. has fast iof i/o
c Vesrion 2.0. Refurbished 05/20/1991
c Modified for the maternal effects 4-5-1990.
c
c copyright 1989-1991 Ignacy Misztal and the University of Illinois
c Use of this program should be acknowledged in publications.
c Free use granted for noncommercial applications.
c
      integer mx,mxc,mxabs,irew,iopen,iwrite,iread
      parameter (mx=30,mxc=5,mxabs=100)
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      real*8 sol(0:mx+mxabs),rs(mx),absp(mx)
      real*8 x(mx,mx),work(mx),xlhs(mxc,mxc),xrhs(mxc)
      real*8 var(0:20),buf(30),r(-10:20)
      real*8 uu,um,mm,truu,trum,trmm,siguu,sigmm,sigum,varum,beta,
     +     rel,vabsp,dabsp,yabsp,yy,yab,xx,ww,y,ss,sum,sige,
     +     tr,denom1,denom2,denom3,denom,sigu
      logical mgs
      integer ifact(-10:20),ilev(0:20),ilev1(0:20),ipos(0:20),ifil(30)
      character*15 npar,ndat,nrel,nout
      integer nunin,nunout,icat,ith,ith1,ibeg,nooffa,itr,imax,iposmx,
     +        isize,irand,ifcat,iabsp,ivarc,ivarc1,ivarc2,n,n1,nr,
     +        j,ipoint,k,nherd,ii,iherd,jj,j1,k1,iivarc,j2,iii,i,
     +        irank,iv,ij

c
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 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,*)'number of categories?'
      read(nunin,*)icat
      ith=icat-1
      ith1=1-ith
      ibeg=1-min(1,ith)
      ipos(0)=0
      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 of fixed'
      write(nunout,*)
      iposmx=itr
      isize=ith
      ilev(0)=ith
      ilev1(0)=0
      var(0)=0
      do 7 i=1,ith
7       ifact(i-ith)=i
 
      irand=0
      do 10 i=1,nooffa
        read(nunin,*)ipos(i),ilev(i),var(i)
        if (var(i).ne.0) irand=irand+ilev(i)
        iposmx=max(iposmx,ipos(i))
        ilev1(i)=isize
10      isize=isize+ilev(i)
c
      if (icat.gt.1) then
         print*,' zero first threshold (1),'
         print*,'  or first level of second and later fixed effects (2)'
         print*,'  or first level of all fixed effects (3),'
         print*,' 0 if none'
         read*,ifcat
      endif
c      
      write(*,*)'position of factor to be absorbed and its variance rati
     *o, 0 0 if none'
      read(*,*)iabsp,vabsp
      iposmx=max(iposmx,iabsp)
      ilev1(nooffa+1)=isize
 
      write(*,*)'how many rounds of reml, 0 if not used'
      read(*,*)ivarc
      if (ivarc.gt.0 .and. icat.gt.1) then
          write(*,*)'how many rounds of newton-raphson before reml'
     *    ,' starts?'
          read(*,*)ivarc1
          write(*,*)'how many rounds of newton-raphson between'
     *    ,' subsequent reml estimations?'
          read(*,*)ivarc2
      endif
      if (isize.gt.mx) then
           print*,'too small matrix, increase mx parameter'
           stop
      endif
      if (icat.gt.mxc) then
         print*,'too many categories, increase mxc parameter in main'
     *   ,'program and subroutine cat'
      endif
c 
c maternal effect model if: two last effecsts are random, 
c have the same number of levels, and the first has a negative variance ratio
      mgs=var(nooffa-1).lt.0 .and. ilev(nooffa-1).eq.ilev(nooffa) 
     +    .and. var(nooffa).gt.0
      if (mgs) varum=-.1*(abs(var(nooffa-1))+var(nooffa))
      print*, ' beta=?'         
      read*,beta
c      
      do 20 i=1,nooffa
20     r(i)=1
c
c 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+min(iabsp,1)
      if (iabsp.ne.0) ifil(n1)=buf(iabsp)
      call iof(iwrite,ifil,n1,4)
      goto 22
c
24    nr=1
      if (nrel.eq.'scratch') goto 26
27    read(3,*,end=26)(ifil(i),i=1,3)
      nr=nr+1
      call iof1(iwrite,ifil,3,8)
      goto 27
 
c open file for absorption variables so that backsolution is possible
26    open(9,form='unformatted')
c open file for coefficient matrix to be used for reml
      open(10,form='unformatted')
 
      do 28 i=0,mx+mxabs
28      sol(i)=0
 
      do 29 i=1,mxc
        xrhs(i)=0
          do 29 j=1,mxc
29          xlhs(i,j)=0
      ifact(nooffa+2)=0
      siguu=0
      sigum=0
      sigmm=0
      i=1
c
c initialize pointer to first solution of absorbed factor
25    ipoint=ilev1(nooffa+1)+1
      do 30 j=1,isize
        rs(j)=0
        absp(j)=0
        do 30 k=1,isize
30        x(j,k)=0
      dabsp=0
      yabsp=0
      yy=0
      nherd=0
 
c do one round of newton-raphson
      call iof(irew,ifil,0,4)
      rewind 9
      do 60 ii=1,n
        call iof(iread,ifil,n1,4)
        do 62 j=1,n1
62        ifact(j)=ifil(j)
 
        if (iabsp.eq.0) then
            yab=0
          elseif (ii.eq.1) then
c absorption
c beginning of data
            iherd=ifact(n1)
            nherd=1
            if (icat.gt.1) yab=sol(ipoint)
          elseif (iherd.ne.ifact(n1)) then
c end of herd
             iherd=ifact(n1)
             nherd=nherd+1
             dabsp=dabsp+vabsp
             do 55 j=1,isize
               if (absp(j).eq.0) goto 55
               xx=absp(j)/dabsp
               rs(j)=rs(j)-xx*yabsp
               do 54 k=1,isize
54               x(j,k)=x(j,k)-xx*absp(k)
55             continue
             ipoint=ipoint+1
             if (icat.gt.1) yab=sol(ipoint)
             write(9)(absp(j),j=1,isize),dabsp,yabsp
             do 56 j=1,isize
56             absp(j)=0
             dabsp=0
             yabsp=0
         endif
         ww=1
c
        do 50 j=1,nooffa
50        if (ifact(j).ne.0) ifact(j)=ifact(j)+ilev1(j)
 
          if (icat.eq.1) then
c linear model
              y=float(ifact(nooffa+1))
            else
c threshold model
c  recalculate weights, pseudo-yields and -regressions
              call categ(ifact,icat,nooffa,sol,ww,r,y,yab,xlhs,xrhs,mx)
          endif
          yy=yy+y*y
          if (iabsp.ne.0) then
                yabsp=yabsp+y*ww
                dabsp=dabsp+ww
          endif
          do 60 jj=ith1,nooffa
            j1=ifact(jj)
            if (j1.eq.0) goto 60
            rs(j1)=rs(j1)+y*ww*r(jj)
            if (iabsp.ne.0) absp(j1)=absp(j1)+ww*r(jj)
            do 61 k=ith1,nooffa
              k1=ifact(k)
              if (k1.eq.0) goto 61
              x(j1,k1)=x(j1,k1)+ww*r(jj)*r(k)
61            continue
60          continue
 
c add extra matrix and vectors to left- and right-hand side of threshold
c sub-equation system
      if (icat.ge.2) then
          do 65 j=1,ith
            rs(j)=rs(j)+xrhs(j)
            xrhs(j)=0
            do 65 k=1,ith
              x(j,k)=x(j,k)+xlhs(j,k)
              xlhs(j,k)=0
65            continue
      endif
c
c absorb last level
      if (nherd.gt.mxabs) then
          print*,'too many levels of absorbed factor, increase '
     *         ,'parameter mxabs'
          stop
      endif
      if (iabsp.ne.0) then
          dabsp=dabsp+vabsp
          do 63 j=1,isize
            rs(j)=rs(j)-absp(j)*yabsp/dabsp
            do 63 k=1,isize
63            x(j,k)=x(j,k)-absp(j)*absp(k)/dabsp
          write(9)(absp(j),j=1,isize),dabsp,yabsp
      endif
 
      iivarc=0
      if (ivarc.ge.1) then
          rewind 10
          write(10)((x(j,k),j=1,isize),k=1,isize)
      endif
67    continue
c beginning of a variance component loop
 
      if (iivarc.gt.0) then
          rewind 10
          read(10)((x(j,k),j=1,isize),k=1,isize)
      endif
      do 110 jj=1,nooffa
        j1=ilev1(jj)+1
        j2=ilev(jj)+ilev1(jj)
 
c add variance ratio to diagonals of random factors
        if (var(jj).gt.0) then
              do 70 k=j1,j2
70              x(k,k)=x(k,k)+var(jj)
           elseif (var(jj).lt.0) then
c ---------- add sire - maternal grandsire matrix
              if (mgs) then
c -------------- maternal effect              
            call addmgs(x,mx,nr,ilev1(jj),ilev1(jj),ilev(jj),-var(jj))
            call addmgs(x,mx,nr,ilev1(jj+1),ilev1(jj),ilev(jj),varum)
            call addmgs(x,mx,nr,ilev1(jj),ilev1(jj+1),ilev(jj),varum)
            call addmgs
     +            (x,mx,nr,ilev1(jj+1),ilev1(jj+1),ilev(jj),var(jj+1))
            goto 80
            else
c ------------- no maternal effects            
            call addmgs(x,mx,nr,ilev1(jj),ilev1(jj),ilev(jj),-var(jj))
            endif
        endif    
110   continue
c
c solve,  add constraints if applicable
80    if (ifcat.ge.2) then
          if (ifcat.eq.2) then
              iii=ibeg+1
            else
              iii=ibeg
          endif     
          do 90 ii=iii,nooffa
            if (var(ii).eq.0) then
               ij=ilev1(ii)+1
               do 92 jj=1,isize
                 x(jj,ij)=0
92               x(ij,jj)=0
            endif
90        continue
        elseif (ifcat.eq.1) then
         do 95 ii=1,isize
           x(ii,1)=0
95         x(1,ii)=0
      endif
      call inv(mx,isize,x,sol,rs,work,irank)
c
        ss=0
c obtain backsolution for
        rewind 9
        do 120 j=1,nherd
          read(9)(absp(k),k=1,isize),dabsp,yabsp
          sum=0
          do 125 k=1,isize
            ss=ss+sol(k)*absp(k)*yabsp/dabsp
125         sum=sum+absp(k)*sol(k)
          sol(ilev1(nooffa+1)+j)=(yabsp-sum)/dabsp
120       ss=ss+sol(ilev1(nooffa+1)+j)*yabsp
 
      if (ivarc.le.iivarc) goto 500
      if (icat.gt.1) then
         if (i.le. ivarc1 .or.mod(i-1-ivarc1,ivarc2).ne.0) goto 500
      endif
c calculate reduction due to non-absorbed factors
      do 130 j=1,isize
130     ss=ss+sol(j)*rs(j)
 
c variance component estimation
      if (icat.gt.1) then
          sige=1
        else
          if (vabsp.eq.0) then
              sige=(yy-ss)/(n-nherd-irank+irand)
            else
              sige=(yy-ss)/(n-irank+irand)
          endif
      endif
      write(*,135)iivarc+1,sige
      write(7,135)iivarc+1,sige
135   format(' vcround=',i2,'  sige=',f10.4)
c
      do 200 iv=1,nooffa
        if (var(iv).eq.0) goto 200
        tr=0
        uu=0
        j1=ilev1(iv)+1
        j2=ilev1(iv)+ilev(iv)
        if (var(iv).gt.0) then
            do 140 j=j1,j2
140            tr=tr+x(j,j)
            do 150 j=j1,j2
150            uu=uu+sol(j)*sol(j)
          else
c           calculate uAu and tr(A-1*X)
            if (mgs) then 
c ---------- maternal effects
                 call uautr(uu,truu,x,sol,mx,nr,ilev1(iv),ilev1(iv))
                 call uautr(um,trum,x,sol,mx,nr,ilev1(iv+1),ilev1(iv))
                 call uautr(mm,trmm,x,sol,mx,nr,ilev1(iv+1),ilev1(iv+1))
c------------------first round EM_REML because sig* are all 0
                 if (siguu.eq.0) then                 
                     siguu=(uu+truu)/ilev(iv)
                     sigum=(um+trum)/ilev(iv)
                     sigmm=(mm+trmm)/ilev(iv)
                   else  
                     denom1=ilev(iv)-truu*sige/siguu
                     siguu=(uu-siguu*beta*denom1)/denom1/(1-beta)
                     denom2=ilev(iv)-trum*sige/sigum
                     sigum=(um-sigum*beta*denom2)/denom2/(1-beta)
                     denom3=(ilev(iv)-trmm*sige/sigmm)
                     sigmm=(mm-sigmm*beta*denom3)/denom3/(1-beta)
                 endif    
                 denom=siguu*sigmm-sigum**2
                 var(iv)=-sige*sigmm/denom
                 varum=sige*(-sigum)/denom
                 var(iv+1)=sige*siguu/denom
                 write(*,170)iv,siguu,iv,var(iv)
                 write(*,160)sigum,varum
                 write(*,170)iv+1,sigmm,iv+1,var(iv+1)
                 write(7,170)iv,siguu,iv,var(iv)
                 write(7,160)sigum,varum
                 write(7,170)iv+1,sigmm,iv+1,var(iv+1)
160              format('  sigcov =',f10.4,'  alfacov =',f8.4)
                 goto 210                 
              else
                 call uautr(uu,tr,x,sol,mx,nr,ilev1(iv),ilev1(iv))
          endif       
        endif
        denom=ilev(iv)-tr*abs(var(iv))
        sigu=(uu- (sige/var(iv)) *beta*denom)/denom/(1-beta)
        var(iv)=sige/sigu*dsign(1.d0,var(iv))
        write(*,170)iv,sigu,iv,sige/sigu
        write(7,170)iv,sigu,iv,sige/sigu
170     format('  sigu',i2,' =',f10.4,'  alfa',i2,' =',f8.4)
200     continue
210   iivarc=iivarc+1
      if (iivarc.lt.ivarc) goto 67
      if (icat.eq.1) then
         print*,' how many vc iterations more, 0-none?'
         read*,j
         if (j.gt.0) then
           ivarc=ivarc+j
           goto 67
         endif
      endif      
c
c
500   write(*,510)i,sol(2),isize,sol(isize)
      write(7,510)i,sol(2),isize,sol(isize)
      if (icat.eq.1) goto 537
510   format(' round=',i3,'   sol(2)=',f10.4,'   sol(',i5,')=',f10.4)
      i=i+1
      if (i.le.imax ) goto 25
      write(*,535)
535   format(' how many more iterations, 0-none?')
      read(*,*)j
      if (j.gt.0) then
      imax=imax+j
      goto 25
      endif
c
c write down solutions
537   write(7,*)
      write(7,*) n,' records',isize,' equations   rank',irank
      if (nherd.gt.0) write(7,540)iabsp,nherd,vabsp
540   format(' absorbed factor',i3,i8,' levels, var. ratio=',f6.2)
      do 550 i=ibeg,nooffa
        write(7,*)
        write(7,560)ipos(i),var(i),ilev(i)
560     format('solutions of factor',i2,'  var. ratio:',f5.2
     *  ,'  levels:',i6)
        if (ilev(i).lt.100) write(7,570)(j,sol(j+ilev1(i)),j=1,ilev(i))
550     continue
570   format(4(3x,i4,':',f11.4))
      open(21,file='sol')
      do 180 i=1,isize
180     write(21,'(2f15.8)')sol(i),x(i,i)
      print*,'solutions and inv. diagonals were stored in file sol'
      close(4,status='delete')
      close(8,status='delete')
      close(9,status='delete')
      close(10,status='delete')
      stop
      end

      subroutine categ(ifact,icat,nvar,sol,ww,r,y,yab,xlhs,xrhs,mx)
c calculates threshold pseudo-values. modified for correct mgs
c model 9/2/87.
      integer nd,mxc,mx
      real*8 xpmin,ww,y
      parameter (nd=2000,mxc=5,xpmin=1e-10)
      real*8 xp(mxc+3),l(mxc+3),p(mxc+3)
      real*8 dist(nd),dint(nd),fun,xlhs(mxc,mxc),xrhs(mxc)
      real*8 t(mxc+3,mxc+3),eta(mxc+3),xint(mxc+3),xdis(mxc+3)
      real*8 sol(0:mx),r(-10:20)
      integer ifact(-10:20),icat,nvar,ifirst,i,j,k,itt
      real*8 xmi,wwsq,w,xl1,xl2,v,yab
      save ifirst,dist,dint

      data ifirst/0/
      
      if (ifirst.eq.0) then
          call setdst(nd,dist,dint)
          ifirst=1
          do 30 i=1,icat-1
30          sol(i)=(i-1)/(icat-1.)
          xint(1)=0
          xint(icat+1)=1
          xdis(1)=0
          xdis(icat+1)=0
      endif
 
      xmi=yab
      wwsq=sqrt(ww)
      do 40 i=1,nvar
40      xmi=xmi+sol(ifact(i))*r(i)
      do 53 j=2,icat
        eta(j)=wwsq*(sol(j-1)-xmi)
        xint(j)=fun(eta(j),nd,dint)
53      xdis(j)=fun(eta(j),nd,dist)
 
      itt=ifact(nvar+1)

      do 55 j=1,icat
        p(j)=0
        xp(j+1)=xint(j+1)-xint(j)
        xp(j+1)=max(xp(j+1),xpmin)
55      if (abs(xp(j+1)).lt.xpmin) xp(j+1)=sign(xpmin,xp(j+1))
      v=(xdis(itt)-xdis(itt+1))/xp(itt+1)*wwsq
      if (itt.eq.1) then
           p(1)=xdis(2)/xp(2)*wwsq
         elseif (itt.eq.icat) then
           p(icat-1)=-xdis(icat)/xp(icat+1)*wwsq
         else
          p(itt-1)=-xdis(itt)/xp(itt+1)*wwsq
          p(itt)=xdis(itt+1)/xp(itt+1)*wwsq
      endif

      w=0
      do 57 j=2,icat+1
57    w=w+(xdis(j-1)-xdis(j))**2/xp(j)
      do 58 j=2,icat
      xl1=(xdis(j)-xdis(j-1))/xp(j)
      xl2=(xdis(j+1)-xdis(j))/xp(j+1)
58    l(j)=-xdis(j)*(xl1-xl2)
      do 59 j=2,icat
      t(j,j)=(xp(j)+xp(j+1))/(xp(j)*xp(j+1))*xdis(j)**2*ww
      k=j+1
      if (k.gt.icat) goto 59
      t(j,k)=-xdis(j)*xdis(k)/xp(k)*ww
      t(k,j)=t(j,k)
59    continue
 
      y=0
      do 90 i=2-icat,0
90      r(i)=l(i+icat)/w
      ww=ww*w
      do 100 i=2-icat,nvar
100     y=y+r(i)*sol(ifact(i))
      y=y+v/w+yab
      do 120 i=1,icat-1
        xrhs(i)=xrhs(i)+p(i)-r(i-icat+1)*ww*(y-xmi)
          do 120 j=1,icat-1
            xrhs(i)=xrhs(i)+t(i+1,j+1)*sol(j)
120         xlhs(i,j)=xlhs(i,j)+t(i+1,j+1)-r(i-icat+1)*r(j-icat+1)*ww
      return
      end
      subroutine inv(n,m,x,sol,rs,work,irank)
c equation system x*sol=rs is solved using generalized inverse, in
c memory. matrix x is of order n*n, sol and rs are of order n, only m
c initial equations are solved. w is a work vector. irank contains the
c rank of inverted matrix.
      real*8 sol(0:n),rs(n)
      real*8 x(n,n),work(n),acc
      integer i,j
      call ginv1(x,m,n,1.0d-8,irank,work)
      do 10 i=1,m
        acc=0
        do 20 j=1,m
20        acc=acc+x(i,j)*rs(j)
10      sol(i)=acc
      return
      end
 
      subroutine ginv1(a,n,m,tol,irank,w)
c returns generalized inverse of matrix x of size n x n declared
c as m x m. tol is working zero (real*8) and irank returns the rank of
c the matrix. w is a work vector of size m,
c by rohan fernando, slightly structured by i. misztal 05/05/87
      real*8 a(m,m),w(m),re,sum,tol
      integer i,j,ii,iim1,iii
      irank=n
      do 10 i=1,n
         do 20 j=1,i-1
              re=a(i,j)
              do 20 ii=i,n
20                 a(ii,i)=a(ii,i)-re*a(ii,j)
         if (a(i,i).lt.tol) then
              a(i,i)=0.0
              do 45 ii=i+1,n
45                 a(ii,i)=0.
           irank=irank-1
           else
              a(i,i)=sqrt(a(i,i))
              do 40 ii=i+1,n
40                a(ii,i)=a(ii,i)/a(i,i)
         endif
10    continue
 
      do 100 i=1,n
         if (a(i,i).eq.0.) then
              do 150 ii=i+1,n
150                a(ii,i)=0.
           else
              a(i,i)=1.0/ a(i,i)
              do 200 ii=i+1,n
200               w(ii)=0.0
              do 300 ii=i+1,n
                  iim1=ii-1
                  re=a(iim1,i)
                  do 400 iii=ii,n
400                   w(iii)=w(iii)-a(iii,iim1)*re
                  if (a(ii,ii).eq.0.) then
                      a(ii,i)=0.
                    else
                      a(ii,i)=w(ii)/a(ii,ii)
                  endif
300           continue
          endif
100     continue
 
      do 110 j=1,n
         do 110 i=j,n
              sum=0
              do 130 ii=i,n
130                sum=sum+a(ii,j)*a(ii,i)
110           a(i,j)=sum
      do 600 i=1,n
          do 600 j=i,n
600           a(i,j)=a(j,i)
      return
      end
      subroutine setdst(n,dist,dint)
c sets n-point look-up tables for the normal distribution and the
c normal integral
      real*8 dist(n),dint(n),sum,d,dd,x,delta
      integer i
      sum=0
      d=1/sqrt(2*3.14159265358979d0)
      delta=24./n
      x=-12
      do 10 i=1,n
        dist(i)=d*dexp(-x*x/2)
        dint(i)=sum
        dd=d*dexp(-(x+delta/2)**2/2)
        x=x+delta
10      sum=sum+delta*dd
      return
      end
 
      function fun(x,n,table)
c return the value of a function in point x stored in table(n) look-up
c table
      integer n,nr
      real*8 x,table(n),anr,fun

      anr=1+(x+12)/24.*n
      nr=int(anr)
      if (nr.le.1) then
        fun=table(1)
       else
        if (nr.gt.n-1) then
       fun=table(n)
       else
        fun=table(nr)+(anr-nr)*(table(nr+1)-table(nr))
      endif
      endif
      return
      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
      integer irew,iopen,iwrite,iread,bufs
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      parameter(bufs=1024)
      integer x(bufs),vector(n),fpos,bpos,iwr,n1,i,lastit
      save x,fpos,bpos,iwr,lastit
      data fpos,bpos/2*-1/
c**********************************************************************
      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
c
      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
c
      elseif (itype.eq.iopen) then
        bpos=0
        fpos=0
        open(iun,access='sequential',form='unformatted')
c
      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
c
      else
        write(*,*)'unknown i/o code=',itype
        stop
      endif
c**********************************************************************
      lastit=itype
      return
50    write(*,*)'end of file',iun
      write(*,*)'bpos,fpos,n',bpos,fpos,n
      stop
      end
 
      subroutine iof1(itype,vector,n,iun)
      integer irew,iopen,iwrite,iread,bufs
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      parameter(bufs=1024)
      integer x(bufs),vector(n),fpos,bpos,iwr,n1,i,lastit
      save x,fpos,bpos,iwr,lastit
      data fpos,bpos/2*-1/
c**********************************************************************
      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
c
      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
c
      elseif (itype.eq.iopen) then
        bpos=0
        fpos=0
        open(iun,access='sequential',form='unformatted')
c
      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
c
      else
        write(*,*)'unknown i/o code=',itype
        stop
      endif
c**********************************************************************
      lastit=itype
      return
50    write(*,*)'end of file',iun
      write(*,*)'bpos,fpos,n',bpos,fpos,n
      stop
      end
 
        subroutine addmgs(x,mx,nr,irow,icol,level,var)
c add the sire-maternal grandsire relationship matrix to x beginning with
c row irow+1 and column icol+1. The sire-mgs list is stored by the iof1 routine
c and has nr entries. nr should equal lev.
      integer irew,iopen,iwrite,iread
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      real*8 x(mx,mx)
      real*8 var,coef(3),weight
      integer i,j,ifil(3),ib,isr,mgs,k
      equivalence (ifil(1),ib),(ifil(2),isr),(ifil(3),mgs)
      data coef/1.,-.5,-.25/
c
        call iof1(irew,ifil,0,8)
c        
        do 10 k=1,nr-1
          call iof1(iread,ifil,3,8)
          if (max(ib,isr,mgs).gt.level .or.min(ib,isr,mgs).lt.0) 
     +                                                         then
              print*,'relationship record no. ',nr,'incorrect'
              print*,'bull,sire,mgs=  ',ib,isr,mgs
              stop
          endif    
c
c calculate the weight factor in the relationships
        weight=1
        if (isr.ne.0) weight=.75
        if (mgs.ne.0) weight=weight-1./16
        weight=1/weight
        do 20 i=1,3
          if (ifil(i).ne.0) then
            do 30 j=1,3
              if (ifil(j).ne.0) then
                 x(ifil(i)+irow,ifil(j)+icol)=
     +              x(ifil(i)+irow,ifil(j)+icol)+coef(i)*coef(j)*weight
     +              *var 
                endif
30           continue
          endif
20      continue                        
10    continue
      end
c
c
      subroutine uautr(uau,tr,x,sol,mx,nr,irow,icol)
c calculates uiAuj and tr(A-1 Xij) ,having solutions sol, matrix x and a 
c list of relationships through the iof1 subroutine.       
c subvectors of u and submatrix of x are defined by irow and icol.
      integer irew,iopen,iwrite,iread,mx,nr,irow,icol
      parameter(irew=1,iopen=2,iwrite=3,iread=4)
      real*8 uau,tr,x(mx,mx)
      real*8 sol(0:mx),coef(3),weight
      integer ifil(3),ib,isr,mgs,i,j,k
      equivalence (ifil(1),ib),(ifil(2),isr),(ifil(3),mgs)
      data coef/1.,-.5,-.25/
c
      uau=0
      tr=0
      call iof1(irew,ifil,0,8)
      do 10 k=1,nr-1
        call iof1(iread,ifil,3,8)
c calculate the weight factor in the relationships
        weight=1
        if (isr.ne.0) weight=.75
        if (mgs.ne.0) weight=weight-1./16
        weight=1/weight
        do 20 i=1,3
          if (ifil(i).ne.0) then
             do 30 j=1,3
                if (ifil(j).ne.0) then
                   uau=uau+coef(i)*coef(j)*weight*sol(ifil(i)+irow)*
     +                 sol(ifil(j)+icol)
                   tr=tr+coef(i)*coef(j)*weight*
     +                x(ifil(i)+irow,ifil(j)+icol)
                endif
30           continue
          endif
20      continue                        
10    continue
      end