******************************************************************** c Calculates variance components for multiple trait problems. c It uses the canonical transformation, is exact with 1 random c effect and approximate with more. c c This version accepts matrices not of full rank c c copyright (C) 1989-1996 Ignacy Misztal 9-1-1989 -- 10/01/96 ********************************************************************* * Description of parameters * * nf - maximum number of factors * ntrait - maximum number of traits * nrand - maximum number of random 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 i-j-a form, the first one is * reused for the sparse matrix solver ********************************************************************** integer nhash,mx,maxss,m1,nhashr,m1r,nf,ntrait,nrand,optim parameter (nhash=30000,mx=5000,m1=1.5*nhash,ntrait=14,nrand=4) parameter (nhashr=.8*nhash,m1r=2*nhashr, optim=1) 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 smpak--------------------------------- real*8 smprsp(1),tmp(mx) integer smpisp(1),smpesp c_______________ storage for variance components________________ integer icanon,rank real*8 w(ntrait,ntrait),winv(ntrait,ntrait), + uiAuj(ntrait,ntrait,nrand),red(ntrait,ntrait), + G(ntrait,ntrait,nrand),r(ntrait,ntrait), + scr1(ntrait,ntrait,nrand),scr2(ntrait,ntrait), + yy(ntrait,ntrait),y11y(ntrait),tracematrix real*8 tr(ntrait,nrand),var(nf),varw(ntrait,nrand),convr, + denom(ntrait,nrand),convg(nrand),sumr,sumg(nrand), + sigmae(ntrait,ntrait),sigma(ntrait,ntrait,nrand), + tracediag,dentot(ntrait) c_______________________________________________________________ c--------------- general storage-------------------------------- real*8 sol(mx,ntrait),rs(mx,ntrait) real second integer ifact(nf),ilev(nf),ilev1(nf),ipos(nf),itrno integer nooffa,itr,iposmx,nfix,irand,itmax,tt,jj, + i,ii,j,k,k1,l,iflag,n,nr,igroup,mg,isize,stat, + ranind(nrand),ibeg,iend,fratio 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 and verify parameters call parinp1(nooffa,itr,itrno,iposmx,nf,ipos,ilev,ilev1,ireg, +npar,var,isize,mx,ntrait,optim) c c initialize call zero(mx,ria) call zerof(3*nhashr,rx) call zerof(ntrait,y11y) call zerof(ntrait**2,yy) n=0 nr=0 igroup=0 itmax=1 ix=0 ixr=0 tmp(1)=second() c create the coefficient matrix in hash table + V call CoefV1(n,nf,iposmx,nooffa,ilev,ilev1,ireg,ipos,ntrait,yy, + y11y,itrno,itr,x,nhash,ix) * add variance ratio to diagonals of random factors do jj=1,nooffa if (var(jj).gt.0) then do k=ilev1(jj)+1,ilev(jj)+ilev1(jj) call hashmx(1d0,k,k,x,nhash,ix) enddo elseif (var(jj).lt.0) then * add relationship matrix if such exists 75 read(3,*,iostat=stat)(ifact(j),j=1,4) if (stat.eq.0) then nr=nr+1 call relhash(ifact,rx,nhashr,ixr,1.0d0,storfull) do j=1,3 ifact(j)=ifact(j)+ilev1(jj) enddo call relhash(ifact,x,nhash,ix,1d0,storhalf) goto 75 endif endif enddo * 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' if (n.eq.0) then print*,'empty data file' stop endif if (igroup.lt.0) then print*,'more animals in relationship file than declared' stop endif if (ixr.eq.0) then print*,'empty relationship file' stop endif c number of last equations to zero out (one group if exists) mg=min(1,igroup) isize=isize-mg c index and count random effects irand=0 nfix=0 do i=1,nooffa if (var(i).ne.0) then if (var(i).lt.0) then if (i.ne.nooffa) then print*,'animal effect must be last!' stop endif nfix=nfix+igroup-mg igroup=igroup-mg ilev(i)=ilev(i)-mg endif irand=irand+1 ranind(irand)=i else nfix=nfix+ilev(i) endif enddo if (irand.gt.nrand) then print*, 'increase parameter nrand' stop endif 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' call hashia(rx,nhashr,ilev(nooffa),ria,rja,ra,m1r,tmp) print*,'ra used ',ria(ilev(nooffa)+1), $ ' out of ',m1r,' words' print*,'before odrv',second() c verify that matrices have positive diagonals call checkdiag(isize,'coefficient',smpia,smpja,smpa) call checkdiag(ilev(nooffa),'relationship',ria,rja,ra) c call prspar(ra,ria,rja,ilev(nooffa)) c call prspar(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 fspak 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,irank) 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 c read g and r if exist, create otherwise open(31,file='gr'//npar) read(31,*,iostat=stat)(((g(i,j,k),i=1,itrno),j=1,itrno), + k=1,irand),((r(i,j),i=1,itrno),j=1,itrno) if (stat.eq.0) then print*,'initial (co)variances read from parameters'' file' else call zerof(ntrait**2,r) call zerof(ntrait**2*irand,g) do i=1,itrno r(i,i)=1 do k=1,irand g(i,i,k)=1/abs(var(ranind(k))) enddo enddo endif c calculate coefficients of the canonical transformation 152 icanon=0 250 icanon=icanon+1 call tcanon1(ntrait,itrno,irand,g,r,w,winv,varw) print*,'variance ratios after transformation in round ',icanon print '(6g12.4)',((varw(i,j),i=1,itrno),j=1,irand) if (icanon.eq.1) print*,'tcanon',second() rewind 7 write(7,'('' w'')') call prden(w,ntrait,itrno,7,'(6f13.5)') c maximum varw is 10e5 do i=1,itrno do j=1,irand if (abs(varw(i,j)).gt.1e5) varw(i,j)=1e5 enddo enddo 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 tt=1,itrno c---set up the system for the correct variance ratio rewind 10 read(10)isize,nsmp,(smpia(i),i=1,isize+1), + (smpa(i),i=1,nsmp),(smpja(i),i=1,nsmp) do j=1,irand jj=ranind(j) if (var(jj).gt.0) then do k=ilev1(jj)+1,ilev(jj)+ilev1(jj) call addspar(smpa,smpia,smpja,k,k,varw(tt,j)-1d0) enddo elseif (var(jj).lt.0) then call adjadd(smpa,smpia,smpja,ilev1(jj)+1,ilev(jj)+ilev1(jj), + ra,ria,rja,1,varw(tt,j)-1d0,storhalf,storfull) endif enddo c call prspar(smpa,smpia,smpja,isize) 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 (icanon.eq.1 .and. tt.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,tt)=rs(j,tt) enddo call fspak(50,isize,smpia,smpja,smpa,sol(1,tt),iflag,6,30,maxss, + smpesp,smpisp,j,j,j,j,j,j,rank) call fspakerr(iflag) c Invert rewind 10 read(10)isize,nsmp,(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; only for c speed optimization 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 (icanon.eq.1 .and. tt.eq.1) then print*,'sparse inversion time=',second() print '('' used'',f7.2,''% memory'')',smpesp*100./maxss endif call fspakerr(iflag) c calculate traces, and quadratic forms for non-animal effects do j=1,irand jj=ranind(j) if (var(jj).lt.0) then tr(tt,j)=tracematrix(isize,smpia,smpja,smpa,ilev(jj), $ ria,rja,ra,ilev1(jj),ilev1(jj),tmp) c u'Au c print*,'quadratic form',second() denom(tt,j)=ilev(jj)-igroup-tr(tt,j)*varw(tt,j) elseif (var(jj).gt.0) then ibeg=ilev1(jj)+1 iend=ilev1(jj)+ilev(jj) tr(tt,j)=tracediag(isize,smpia,smpja,smpa,ibeg,iend) denom(tt,j)=ilev(jj)-tr(tt,j)*varw(tt,j) endif enddo enddo c--------------------------------------------------------------------------- c transform back the R and G c--------------------------------------------------------------------------- c quadratic forms call zerof(ntrait**2*irand,uiauj) do i=1,irand ii=ranind(i) if (var(ii).gt.0) then do j=1,itrno do k=j,itrno do l=ilev1(ii)+1,ilev1(ii)+ilev(ii) uiauj(j,k,i)= uiauj(j,k,i)+sol(l,j)*sol(l,k) enddo uiauj(k,j,i)=uiauj(j,k,i) enddo enddo else call quadrmt(uiauj(1,1,i),ra,ria,rja,sol,ntrait,itrno,mx,1, + ilev(nooffa),ilev1(nooffa)) endif enddo c print*,'after uiauj',second() c reductions do i=1,itrno do j=1,itrno red(i,j)=0 do k=1,irand red(i,j)=red(i,j)+uiauj(i,j,k)*sqrt(varw(i,k)*varw(j,k)) enddo do k=1,isize c ********************** Is this right? ********************* red(i,j)=red(i,j)+sol(k,i)*rs(k,j) enddo enddo enddo call zerof(ntrait**2*irand,scr1) call zerof(ntrait**2,scr2) c contributions to sigmae denominator do i=1,itrno dentot(i)=0 do j=1,irand dentot(i)=dentot(i)+tr(i,j)*varw(i,j) enddo enddo do i=1,itrno do j=i,itrno do k=1,irand sigma(i,j,k)=uiAuj(i,j,k)/sqrt(denom(i,k)*denom(j,k)) sigma(j,i,k)=sigma(i,j,k) enddo sigmae(i,j)=(yy(i,j)-red(i,j))/ + (n-rank+sqrt(dentot(i)*dentot(j))) sigmae(j,i)=sigmae(i,j) enddo enddo print*,'computed variance ratios:' print '(6g12.4)',((sigmae(i,i)/sigma(i,i,j),i=1,itrno),j=1,irand) do i=1,itrno do j=1,itrno do k=1,itrno do l=1,itrno do k1=1,irand scr1(l,i,k1)=scr1(l,i,k1)+winv(l,k)*sigma(k,j,k1) + *winv(i,j) enddo scr2(l,i)=scr2(l,i)+winv(l,k)*sigmae(k,j)*winv(i,j) enddo enddo enddo enddo convr=0 sumr=0 do i=1,irand sumg(i)=0 convg(i)=0 enddo do i=1,itrno do j=1,itrno convr=convr+(scr2(i,j)-r(i,j))**2 r(i,j)=scr2(i,j) sumr=sumr+r(i,j)**2 do k=1,irand convg(k)=convg(k)+(scr1(i,j,k)-g(i,j,k))**2 g(i,j,k)=scr1(i,j,k) sumg(k)=sumg(k)+g(i,j,k)**2 enddo enddo enddo do i=1,irand convg(i)=convg(i)/sumg(i) enddo convr=convr/sumr write(7,'(''sigma & sigmae'')') do i=1,irand call prden(sigma(1,1,i),ntrait,itrno,7,'(8f13.5)') enddo call prden(sigmae,ntrait,itrno,7,'(8f13.5)') print '('' convR='',g9.2,'' convG='',5(g9.2,2x))',convr, + (convg(i),i=1,irand) write(7,*)'g,r, and heritabilities + correlations' rewind 31 do i=1,irand call prden(g(1,1,i),ntrait,itrno,7,'(6f13.5)') call prden(g(1,1,i),ntrait,itrno,31,'(6f13.5)') enddo call prden(r,ntrait,itrno,7,'(6f13.5)') call prden(r,ntrait,itrno,31,'(6f13.5)') call corr1(ntrait,itrno,irand,g,r,scr1,7) write(7,*)'sigma and sigmae' do i=1,irand call prden(sigma(1,1,i),ntrait,itrno,7,'(6f13.5)') enddo call prden(sigmae,ntrait,itrno,7,'(6f13.5)') c if not converged, repeat if (icanon.lt.itmax) goto 250 print*,'finish? (rounds-more?)' read*,j itmax=itmax+j if (icanon.lt.itmax .and. j.ne.0) goto 250 * print* Print*,' Heritability, repeatability and variances' print*,'trait h^2 r^2 residual effect' print '(24x,5i12)',(i,i=1,min(5,irand)) do i=1,itrno print '(i3,2f8.3,6g12.4)',i,scr1(i,i,1),scr1(i,i,2), + r(i,i),(g(i,i,j),j=1,min(5,irand)) enddo print*,'heritabilities on diagonal, genetic correlations above,', + 'phentotypic 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' stop end subroutine zerof(n,a) c zeroes real*8 matrix a of dimension n. integer n,i real*8 a(n) do i=1,n a(i)=0 enddo end subroutine parinp1(nooffa,itr,itrno,iposmx,nf,ipos,ilev,ilev1, + ireg,npar,var,isize,mx,ntrait,optim) c the parameters can be read either from the console or from a file character*15 npar,ndat,nrel,nout integer nooffa,itr,itrno,iposmx,nf,isize,ipos(nf),ilev(nf), + ilev1(nf),i,nunin,nunout,mx,ntrait,optim logical ireg(nf) real*8 var(nf) 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,6) 6 format(' no of factors, trait pos., no of traits ?') read(nunin,*)nooffa,itr,itrno write(nunout,*)' for each factor give position, no of levels,', + ' and initial variance ratio' write(nunout,*)' the last factor is assumed to be animal', * ' with groups' write(nunout,*) iposmx=itr+itrno-1 isize=0 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 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 (itrno.gt.ntrait) then print*,'too many traits, increase parameter ntrait' stop endif if (itrno.lt.1) then print*,'number of traits <=0' 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 end * * subroutine CoefV1(n,nf,iposmx,nooffa,ilev,ilev1,ireg,ipos,ntrait, + yy,y11y,itrno,itr,x,nhash,ix) c calculates coefficient matrix in hash table in memory integer n,nf,itr,iposmx,nooffa,ilev(nooffa),ilev1(nooffa) + ,ipos(nf),ifact(50),ntrait,itrno,nhash,ix,i,j,j1,k,k1 real*8 x(nhash,3),val(50),buf(50),yy(ntrait,ntrait),y11y(ntrait) logical ireg(nf) c create the coefficient matrix in hash table do 25 i=1,3 do 25 j=1,nhash 25 x(j,i)=0 c ******** while not eof(2) do ************* 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 do 55 i=1,itrno y11y(i)=y11y(i)+buf(itr+i-1) do 55 j=1,itrno 55 yy(i,j)=yy(i,j)+buf(itr+i-1)*buf(itr+j-1) do 60 j=1,nooffa j1=ifact(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 ****** 65 continue end subroutine posdef(g,itrno,ntrait,d,x,ineg) c checks for positive-definitness of G. If not, zeroes c negative eigevalues and restores G. integer itrno,ntrait,ineg,i,j,k real*8 g(ntrait,ntrait),x(ntrait,ntrait),d(itrno) c restore upper-diagonal elements of G destroyed by Jacobi call restor(g,itrno,ntrait) ineg=0 do 10 i=1,itrno if(d(i).lt.0) then ineg=ineg+1 do 20 j=1,itrno do 20 k=1,itrno 20 g(j,k)=g(j,k)-d(i)*x(i,k)*x(i,j) d(i)=1.0e-10 endif 10 continue end subroutine restor(g,itrno,ntrait) c restores symmetric g which upper diagonal elements have been destroyed integer ntrait,itrno,i,j real*8 g(ntrait,ntrait) do 10 i=1,itrno-1 do 10 j=i+1,itrno 10 g(i,j)=g(j,i) end subroutine trrhs(n,isize,iposmx,nooffa,ilev,ilev1,ipos,ireg + ,ntrait,itrno,itr,mx,rs,w,yy) c create transformed right-hand sides + sums of squares integer n,iposmx,nooffa,ilev(nooffa),ilev1(nooffa),ipos(nooffa), + itr,ntrait,itrno,ifact(50),i,j,j1,k,mx,isize real*8 rs(mx,1),buf1(50),val(50),buf(50) real*8 yy(ntrait,ntrait),w(ntrait,ntrait) logical ireg(nooffa) do 20 i=1,itrno do 10 j=1,itrno 10 yy(i,j)=0 do 20 j=1,isize 20 rs(j,i)=0 rewind 2 c ******** while not eof(2) do ************* c do 100 i=1,n 30 read(2,*)(buf(j),j=1,iposmx) do 50 j=1,nooffa ifact(j)=int(buf(ipos(j))) if (ireg(j)) then val(j)=ifact(j) ifact(j)=ilev1(j)+1 else if (ifact(j).le.ilev(j)) then val(j)=1 else val(j)=0 endif ifact(j)=ifact(j)+ilev1(j) endif 50 continue do 55 j=1,itrno buf1(j)=0 do 55 k=1,itrno 55 buf1(j)=buf1(j)+buf(itr+k-1)*w(j,k) do 57 j=1,itrno do 57 k=1,itrno 57 yy(j,k)=yy(j,k)+buf1(j)*buf1(k) do 60 j=1,nooffa j1=ifact(j) do 60 k=1,itrno 60 rs(j1,k)=rs(j1,k)+buf1(k)*val(j) 100 continue end subroutine tcanon1(ntrait,itrno,irand,G,R,W,WINV,varw) c diagonalizes R, and attempts to simultaneously diagonalize nrand c matrices in G integer ntrait,itrno,irand,i,j,k,l,m,nrot,nscr,flag,nsweep,msglev parameter (nscr=20) real*8 g(ntrait,ntrait,irand),r(ntrait,ntrait), + winv(ntrait,ntrait),w(ntrait,ntrait),g1(nscr,nscr,5), + r1(nscr,nscr),x(nscr,nscr),y(20,20,5),z(nscr,nscr),d(nscr), + varw(ntrait,irand),weight(5),scr(nscr),d1(nscr) data weight/5*1d0/,msglev/0/ c bypass procedures for single trait if (itrno.eq.1) then winv(1,1)=sqrt(r(1,1)) w(1,1)=1/winv(1,1) do i=1,irand varw(1,i)=r(1,1)/g(1,1,i) enddo return endif c numbers 20 and 5 are coded in subroutine FG if (nscr.lt.itrno .or. irand.gt.5 .or. nscr.ne.20) then print*,'nscr too low or too many random effects' stop endif c zero and copy matrices do i=1,itrno do j=1,itrno winv(j,i)=0 w(i,j)=0 r1(j,i)=r(j,i) do k=1,irand g1(i,j,k)=g(i,j,k) y(i,j,k)=0 enddo enddo enddo if (msglev.ge.2) then print*,'r and g' call prden(r,ntrait,itrno,6,'(20g10.2)') do i=1,irand call prden(g(1,1,i),ntrait,itrno,6,'(20g10.2)') enddo print*,'d' print '(5g10.2)',(d(i),i=1,itrno) endif c check for positive definitness and zero negative eigenvalues, if needed do i=1,irand c find e-values of G, to be able to check it is positive-definite call jacobi(g1(1,1,i),itrno,nscr,d,x,nrot) call posdef(g1(1,1,i),itrno,nscr,d,x,flag) if (flag.gt.0) then print*,'g(',i,') not positive definite' endif enddo c find e-values of R call jacobi(r1,itrno,nscr,d,x,nrot) c check for positive definitness and zero negative eigenvalues, if needed call posdef(r1,itrno,nscr,d,x,flag) if (flag.gt.0) then print*,'r not positive definite' endif c y()=sqrt(inv(d))*X'*G()*X*sqrt(inv(d)) do i=1,itrno do m=1,irand do l=i,itrno do j=1,itrno do k=1,itrno y(i,l,m)=y(i,l,m)+1/sqrt(d(i))*x(j,i)*g(j,k,m)*x(k,l)/ + sqrt(d(l)) enddo enddo y(l,i,m)=y(i,l,m) enddo enddo enddo if (msglev.ge.2) then print*,'x and y' call prden(x,nscr,itrno,6,'(20f8.2)') do i=1,irand call prden(y(1,1,i),nscr,itrno,6,'(20f8.2)') enddo endif c (semi) diagonalize all y() simulateneously using the FG algorithm by Flury c Constantine c (testing bypass for one random effect) c if (irand.eq.1) then c call jacobi(y,itrno,nscr,d1,z,nrot) c call posdef(y,itrno,nscr,d1,z,flag) c if (flag.gt.0) then c print*,'y(',i,') not positive definite' c endif c do i=1,itrno c y(i,i,1)=d1(i) c enddo c else do i=1,irand call jacobi(y(1,1,i),itrno,nscr,d1,z,nrot) call posdef(y(1,1,i),itrno,nscr,d1,z,flag) if (flag.gt.0) then print*,'y(',i,') not positive definite' endif enddo call fgalgm(y,weight,z,itrno,irand,flag,nsweep) if (flag.ne.0) then print*,'flag=',flag,' in multiple diagonalization' if (flag.le.3) stop endif c endif if (msglev.ge.2) then print*,'z and y after fg' call prden(z,nscr,itrno,6,'(20f8.2)') do i=1,irand call prden(y(1,1,i),nscr,itrno,6,'(20f8.2)') enddo endif c print weigthed off-diagonals for measure of diagonalization if (nscr.gt.irand) then do i=1,irand scr(i)=0 do j=1,itrno do k=j+1,itrno scr(i)=scr(i)+abs(y(j,k,i)**2/y(j,j,i)/y(k,k,i)) enddo enddo scr(i)=sqrt(+1e-20+scr(i))/(itrno*(itrno-1)/2) enddo print 10,(scr(i),i=1,irand) 10 format(' relative off-diagonals in G after diagonalization:', + 5g10.2) endif c calculate the transformation W=Z'D(-.5)X' c also winv = inv(W) = X D(-.5) Z do i=1,itrno do j=1,itrno do k=1,itrno winv(i,k)=winv(i,k)+x(i,j)*sqrt(d(j))*z(j,k) w(i,k)=w(i,k)+z(j,i)/sqrt(d(j))*x(k,j) enddo enddo enddo do i=1,itrno do j=1,irand varw(i,j)=1/y(i,i,j) enddo enddo if (msglev.ge.2) then print*,'w' call prden(w,ntrait,itrno,6,'(20f8.4)') endif end subroutine quadruv(uav,a,ia,ja,u,v,ibeg,iend,ivbeg) c calculates the quadratic form uAv of part of vectors q and u, and part of c sparse matrix a (a,ia,ja). Quadratic form is taken from row-column c ibeg to iend. Vectors u,v are used from ivbeg+1 location integer ia(1),ja(1),ibeg,iend,i,j,ivbeg real*8 uav real*8 a(1),u(1),v(1) uav=0 do 10 i=ibeg,iend do 10 j=ia(i),ia(i+1)-1 10 if (ja(j).ge.ibeg .and. ja(j).le.iend) $ uav=uav+a(j)*u(i+ivbeg)*v(ja(j)+ivbeg) end subroutine prden(x,n,m,un,form) c prints real*8 matrix x(n,n), first m rows and columns integer n,m,i,j,un character*(*) form real*8 x(n,n) do 10 i=1,m 10 write(un,form)(x(i,j),j=1,m) write(un,*) end SUBROUTINE JACOBI(A,N,NP,D,V,NROT) integer nmax,n,np,nrot,ip,iq,i,j PARAMETER (NMAX=100) real*8 A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX) real*8 sm,tresh,g,h,t,theta,c,s,tau DO 12 IP=1,N DO 11 IQ=1,N V(IP,IQ)=0. 11 CONTINUE V(IP,IP)=1. 12 CONTINUE DO 13 IP=1,N B(IP)=A(IP,IP) D(IP)=B(IP) Z(IP)=0. 13 CONTINUE NROT=0 DO 24 I=1,50 SM=0. DO 15 IP=1,N-1 DO 14 IQ=IP+1,N SM=SM+ABS(A(IP,IQ)) 14 CONTINUE 15 CONTINUE IF(SM.EQ.0.)RETURN IF(I.LT.4)THEN TRESH=0.2*SM/N**2 ELSE TRESH=0. ENDIF DO 22 IP=1,N-1 DO 21 IQ=IP+1,N G=100.*ABS(A(IP,IQ)) IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) * .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN A(IP,IQ)=0. ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN H=D(IQ)-D(IP) IF(ABS(H)+G.EQ.ABS(H))THEN T=A(IP,IQ)/H ELSE THETA=0.5*H/A(IP,IQ) T=1./(ABS(THETA)+SQRT(1.+THETA**2)) IF(THETA.LT.0.)T=-T ENDIF C=1./SQRT(1+T**2) S=T*C TAU=S/(1.+C) H=T*A(IP,IQ) Z(IP)=Z(IP)-H Z(IQ)=Z(IQ)+H D(IP)=D(IP)-H D(IQ)=D(IQ)+H A(IP,IQ)=0. DO 16 J=1,IP-1 G=A(J,IP) H=A(J,IQ) A(J,IP)=G-S*(H+G*TAU) A(J,IQ)=H+S*(G-H*TAU) 16 CONTINUE DO 17 J=IP+1,IQ-1 G=A(IP,J) H=A(J,IQ) A(IP,J)=G-S*(H+G*TAU) A(J,IQ)=H+S*(G-H*TAU) 17 CONTINUE DO 18 J=IQ+1,N G=A(IP,J) H=A(IQ,J) A(IP,J)=G-S*(H+G*TAU) A(IQ,J)=H+S*(G-H*TAU) 18 CONTINUE DO 19 J=1,N G=V(J,IP) H=V(J,IQ) V(J,IP)=G-S*(H+G*TAU) V(J,IQ)=H+S*(G-H*TAU) 19 CONTINUE NROT=NROT+1 ENDIF 21 CONTINUE 22 CONTINUE DO 23 IP=1,N B(IP)=B(IP)+Z(IP) D(IP)=B(IP) Z(IP)=0. 23 CONTINUE 24 CONTINUE PAUSE '50 iterations should never happen' RETURN END subroutine corr1(n,m,irand,g,r,cor,un) c calculates heritabilities on the diagonal, phenotypic correlations c below, and genetic correlations above; assumes the animal effect is last c Also in the highest dimension, calculates repeatabilities on the diagonal c (only if irand>1) integer n,m,i,j,k,un,irand real*8 g(n,n,irand),r(n,n),cor(n,n,2),nom,denomi,denomj do i=1,m cor(i,i,1)=0 do j=1,irand cor(i,i,1)=cor(i,i,1)+g(i,i,j) enddo cor(i,i,2)=cor(i,i,1)/(cor(i,i,1)+r(i,i)) cor(i,i,1)=g(i,i,irand)/(cor(i,i,1)+r(i,i)) do j=i+1,m cor(i,j,1)=g(i,j,irand)/sqrt(g(i,i,irand)*g(j,j,irand)) nom=r(i,j) denomi=r(i,i) denomj=r(j,j) do k=1,irand nom=nom+g(i,j,k) denomi=denomi+g(i,i,k) denomj=denomj+g(j,j,k) enddo cor(j,i,1)=nom/sqrt(denomi*denomj) enddo enddo do 20 i=1,m 20 write(un,'(10f8.3)')(cor(i,j,1),j=1,m) write(un,*) end c this file includes asr 71 and asr 74 which provide replacements for c subroutines falg and galg respectively. the replacements have been c placed at the end. c subroutine fgalgm(a, rn, b, ip, k, ifault, nf) c c algorithm as211 appl. statist. (1985) vol. 34, no. 2 c f-g diagonalization algorithm c c modified by I. Misztal for different index order c real*8 a(20, 20, 5), h(5, 20, 20), rn(5), b(20, 20) real*8 zero, one, p8, p6, g1, g2 data zero /0.0/, one /1.0/, p8 /0.8/, p6 /0.6/ c c check input parameters c ifault = 1 if (k .lt. 1 .or. k .gt. 5) return ifault = 2 if (ip .lt. 2 .or. ip .gt. 20) return ifault = 3 do 3 i = 1, k if (rn(i) .le. zero) return 3 continue do 4 i = 1, k ifault = -i do 4 l = 2, ip jend = l - 1 do 4 j = 1, jend if (a(l, j, i) .ne. a(j, l, i)) return 4 continue ifault = 0 c do 5 l = 1, ip do 5 j = 1, ip b(l, j) = zero if (l .eq. j) b(l, j) = one 5 continue c do 7 i = 1, k do 7 l = 1, ip do 7 j = 1, ip 7 h(i, l, j) = a(l, j, i) call falg(a, rn, b, ip, k, ifault, nf) if (ifault .ne. 0 .or. nf .gt. 1) return c c if f-algorithm stopped after only 1 iteration, try c another initial approximation to the orthogonal matrix b. c do 8 i = 1, k do 8 l = 1, ip do 8 j = 1, ip 8 a(l, j, i) = h(i, l, j) ip1 = ip - 1 do 6 l = 1, ip1 do 6 j = 1, ip g1 = p8 * b(l, j) + p6 * b(l + 1, j) g2 = p8 * b(l + 1, j) - p6 * b(l, j) b(l, j) = g1 b(l + 1, j) = g2 6 continue call falg(a, rn, b, ip, k, ifault, nf) return end c c subroutine eigvec(u, q) real*8 u(2, 2), q(2, 2) real*8 zero, one, two, root, denom, ep data zero /0.0/, one /1.0/, two /2.0/, ep /1d-10/ c root = (u(1, 1) + u(2, 2)) / two root = root + sqrt(((u(1, 1) - u(2, 2)) / two) ** 2 + u(1, 2) * ** 2) denom = sqrt((root - u(1, 1)) ** 2 + u(1, 2) ** 2) q(1, 1) = one q(2, 1) = zero if (denom .lt. ep) goto 1 q(1, 1) = u(1, 2) / denom q(2, 1) = (root - u(1, 1)) / denom 1 q(1, 2) = -q(2, 1) q(2, 2) = q(1, 1) return end c subroutine mult(a, b, ip, iq) real*8 a(20, 20), b(20, 20), h(20, 20) real*8 zero data zero /0.0/ c do 1 i = 1, iq do 1 j = 1, ip h(i, j) = zero do 1 k = 1, ip h(i, j) = h(i, j) + b(k, i) * a(k, j) 1 continue do 2 i = 1, iq do 2 j = 1, i a(i, j) = zero do 3 k = 1, ip 3 a(i, j) = a(i, j) + h(i, k) * b(k, j) a(j, i) = a(i, j) 2 continue return end c c--------------------------------------------------------------------- c subroutine falg(a, rn, b, ip, k, ifault, nf) c c asr 71 (remark on as 211) appl. statist. (1988) vol. 37, no. 1 c integer i, ifault, ip, ip1, j, jstart, k, l, m, nf, maxf real*8 a(20, 20, 5), aux(20, 20), b(20, 20), b1, b2, bold(20, 20) real*8 c, c2, diff, eps, epsf, q(2, 2), r1, rn(5), s, t(5, 2, 2) real*8 t1, t2, zero, one, two c STORAGE FOR REORTHOGONALIZATION REAL*8 SCR1(20,20),SCR2(20,20) data zero /0.0/, one /1.0/, two /2.0/, epsf /1d-10/, maxf /15/ c nf = 0 ifault = 0 ip1 = ip - 1 c c initial multiplication c do 2 i = 1, k do 1 l = 1, ip do 1 j = 1, ip 1 aux(l, j) = a(l, j, i) call mult(aux, b, ip, ip) do 2 l = 1, ip do 2 j = 1, ip 2 a(l, j, i) = aux(l, j) c c start iteration step of f-algorithm c 4 nf = nf + 1 do 5 l = 1, ip do 5 j = 1, ip 5 bold(l, j) = b(l, j) do 8 l = 1, ip1 jstart = l + 1 do 8 j = jstart, ip do 6 i = 1, k t(i, 1, 1) = a(l, l, i) t(i, 1, 2) = a(l, j, i) t(i, 2, 1) = a(j, l, i) t(i, 2, 2) = a(j, j, i) 6 continue c c get rotation sine and cosine c call galg(t, q, rn, k, ifault) if (ifault .ne. 0) return c c rotate b c c = q(1, 1) s = q(2, 1) r1 = s / (one + c) do 7 m = 1, ip b1 = b(m, l) b2 = b(m, j) b(m, l) = b1 + s * (b2 - r1 * b1) b(m, j) = b2 - s * (b1 + r1 * b2) 7 continue c c update the a matrices c t1 = s / c t2 = t1 * t1 c2 = c * c do 8 i = 1, k a(l, l, i) = c2 * (t(i, 1, 1) + two * t1 * t(i, 1, 2) + t2 * * t(i, 2, 2)) a(j, j, i) = c2 * (t2 * t(i, 1, 1) - two * t1 * t(i, 1, 2) * + t(i, 2, 2)) a(l, j, i) = c2 * (t1 * (t(i, 2, 2) - t(i, 1, 1)) + (one - t2) * * t(i, 1, 2)) a(j, l, i) = a(l, j, i) do 8 m = 1, ip if (m .eq. j .or. m .eq. l) goto 8 b1 = a(m, l, i) b2 = a(m, j, i) a(m, l, i) = b1 + s * (b2 - r1 * b1) a(m, j, i) = b2 - s * (b1 + r1 * b2) a(l, m, i) = a(m, l, i) a(j, m, i) = a(m, j, i) 8 continue c c check for convergence c C REORTHOGONALIZE, AS SUGGESTED IN THE BOOK FOR LARGE N call qrfact(b,scr1,scr2,ip,20) do i=1,ip do j=1,ip b(i,j)=scr1(i,j) enddo enddo eps = zero do 9 l = 1, ip do 9 j = 1, ip diff = abs(b(l, j) - bold(l, j)) if (diff .gt. eps) eps = diff 9 continue if (eps .ge. epsf .and. nf .lt. maxf) goto 4 if (eps .ge. epsf) ifault = 4 return end c c--------------------------------------------------------------------- c c this is asr 74 c subroutine galg (t, q, rn, k, ifault) real*8 c, co, eight, epsg, fmg, four, h, pi, q(2,2), rn(5), & si, t(5,2,2), theta, two, u, wt, zero c data epsg/1.0d-14/, zero/0.0d0/, two/2.0d0/, four/4.0d0/, & eight/8.0/, pi/3.141592653589793238/ c fmg = zero h = zero c compute iteration constants do 1 i=1, k u = (t(i,2,2)-t(i,1,1))/two c = t(i,1,2) wt = rn(i) fmg = fmg + wt*(u*u-c*c) 1 h = h + wt*u*c c check for degenerate cases if (abs(h) .lt. epsg .or. abs(fmg) .lt. epsg) goto 2 c c nodegenerate cases (no. 1 to 4) c theta = atan(-two*h/fmg)/four c cases 2 and 4 if (h .gt. zero .and. fmg .lt. zero) theta = theta - pi/four if (h .lt. zero .and. fmg .lt. zero) theta = theta + pi/four go to 3 c c degenerate cases (no. 5 to 9) c 2 theta = zero if (fmg .le. -epsg) theta = -pi/four if (h .ge. epsg) theta = -pi/eight if (h .le. -epsg) theta = pi/eight c get cosines and sines from angle 3 si = sin(theta) co = cos(theta) q(1,1) = co q(1,2) = -si q(2,1) = si q(2,2) = co return end subroutine qrfact(a,q,r,n,np) c Modified Gram-Schmidt qr factorization: a=qr for real matrices based c on algorithm 5.2.5 in "matrix computations" by Golub and Van Loan. c Matrices a, q and r are declared as np x np and have size n integer n,np,j,k,l real*8 a(np,np),q(np,np),r(np,np) c for k=1:n c do k=1,n c c r(k,k)=||a(1:m,k)/r(k,k) c r(k,k)=0 do l=1,n r(k,k)=r(k,k)+a(l,k)**2 enddo r(k,k)=sqrt(r(k,k)) c c q(1:m,k)=a(1:m,k)/r(k,k) c do l=1,n q(l,k)=a(l,k)/r(k,k) enddo c c for j=k+1:n c do j=k+1,n c c r(k,j)=q(1:m,k)'a(1:m,j) c r(k,j)=0 do l=1,n r(k,j)=r(k,j)+q(l,k)*a(l,j) enddo c c a(1:m,j)=a(1:m,j)-q(1:m,k)r(k,j) c do l=1,n a(l,j)=a(l,j)-q(l,k)*r(k,j) enddo enddo enddo end