subroutine parinp(nooffa,itr,itrno,iposmx,nf,ipos,ilev,ilev1,ireg, +isize,npar,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) 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,' 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) 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 CoefV(n,nf,iposmx,nooffa,ilev,ilev1,ireg,ipos,ntrait, + yy,y11y,V,itrno,itr,x,nhash,ix) 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) real*8 yy(ntrait,ntrait),y11y(ntrait),v(ntrait,ntrait) logical ireg(nf) c create the coefficient matrix in hash table + V do 25 i=1,3 do 25 j=1,nhash 25 x(j,i)=0 do 27 i=1,itrno do 27 j=1,itrno 27 v(j,i)=0 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 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 c calculate v do 100 i=1,itrno do 100 j=1,itrno 100 v(i,j)=v(i,j)+(yy(i,j)-y11y(i)*y11y(j)/n)/(n-1) end subroutine trinp(tab,ntab,ipnts,n1,inorec,alfa,tr) c interploates and extrapolates trace elements based on a c 2-step function: c trace =n[x/(d+alfa)+y/(e+alfa)],x=.15, y=.85 c Tab contains tabulated ipnts points, inorec is number of points that c contribute to trace n/alfa. integer ntab,ipnts,n,inorec real*8 tab(ntab,3),alfa,tr real*8 t1,t2,alfa1,alfa2,tau1,tau2,d,e,delta,p1,p2,q1,q2,x,y,a,b,c integer k,n1 data x/.15/,y/.85/ * n=n1 c find the right interval k=2 10 if (tab(k,1).lt.alfa .and. k.lt.ipnts) then k=k+1 goto 10 endif alfa1 = tab(k-1,1) alfa2 = tab(k,1) t1 = tab(k-1,2)-inorec/alfa1 t2 = tab(k,2)-inorec/alfa2 n = n-inorec tau1 = t1 / n tau2 = t2 / n p1 = x - tau1 * alfa1 p2 = x - tau2 * alfa2 q1 = x + y - tau1 * alfa1 q2 = x + y - tau2 * alfa2 a = p1 * tau2 - p2 * tau1 b = p1 * (x - q2) + tau2 * alfa1 * q1 - + p2 * (x - q1) - tau1 * alfa2 * q2 c = alfa1 * q1 * (x - q2) - alfa2 * q2 * (x - q1) delta = SQRT(b**2 - 4 * a * c) e = (-b - delta) / (2 * a) d = (e * p1 + alfa1 * q1) / (tau1 * e + tau1 * alfa1 - y) tr= n * (x / (d + alfa) + y / (e + alfa)) + inorec/alfa end subroutine posdef(g,itrno,ntrait,d,x,ineg) c checks for positive-definitness of G. If not, zeroes c negative eigevalues and eigenvectors, 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 tcanon(ntrait,itrno,G,R,W,WINV,x,y,z,varw) c calculates the coefficients of the canonical transformation integer ntrait,itrno,i,j,k,l,nrot,nscr parameter (nscr=20) real*8 g(ntrait,ntrait),r(ntrait,ntrait),winv(ntrait,ntrait), + w(ntrait,ntrait),x(ntrait,ntrait),y(ntrait,ntrait), + z(ntrait,ntrait),d(nscr),d1(nscr) real*8 varw(ntrait) if (nscr.lt.itrno) then print*,'nscr too low' stop endif c zero matrices do 5 i=1,itrno do 5 j=1,itrno y(i,j)=0 winv(i,j)=0 5 w(i,j)=0 c find e-values of G, to be able to check it is positive-definite call jacobi(g,itrno,ntrait,d,x,nrot) c check for positive definitness call posdef(g,itrno,ntrait,d,x,i) if (i.gt.0) then print*,'g not positive definite' endif c find e-values of R call jacobi(r,itrno,ntrait,d,x,nrot) c check for positive definitness call posdef(r,itrno,ntrait,d,x,i) if (i.gt.0) then print*,'r not positive definite' endif c y=sqrt(inv(d))*X'*G*X*sqrt(inv(d)) do 10 i=1,itrno do 10 j=1,itrno do 10 k=1,itrno do 10 l=1,itrno 10 y(i,l)=y(i,l)+1/sqrt(d(i))*x(j,i)*g(j,k)*x(k,l)/ + sqrt(d(l)) c find e-values of y call jacobi(y,itrno,ntrait,d1,z,nrot) c check for positive definitness call posdef(y,itrno,ntrait,d1,z,i) if (i.gt.0) then print*,'y not positive definite' endif c calculate the transformation W=Z'D(-.5)X' c also winv = inv(W) = X D(-.5) Z do 20 i=1,itrno do 20 j=1,itrno do 20 k=1,itrno winv(i,k)=winv(i,k)+x(i,j)*sqrt(d(j))*z(j,k) 20 w(i,k)=w(i,k)+z(j,i)/sqrt(d(j))*x(k,j) do 30 i=1,itrno 30 varw(i)=1/d1(i) 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 corr(n,m,g,r,cor,un) c calculates heritabilities on the diagonal, residual correlations c below, and genetic correlations above integer n,m,i,j,un real*8 g(n,n),r(n,n),cor(n,m) do 10 i=1,m cor(i,i)=g(i,i)/(r(i,i)+g(i,i)) do 10 j=i+1,m cor(i,j)=g(i,j)/sqrt(g(i,i)*g(j,j)) 10 cor(j,i)=r(j,i)/sqrt(r(i,i)*r(j,j)) do 20 i=1,m 20 write(un,'(8f10.3)')(cor(i,j),j=1,m) write(un,*) end subroutine stderr(r,g,ntrait,itrno,trtab,ntab,ipnts,inorec, + n,na,nfix,un) c calculates standard errors of residual+genetic variances+heritabilities c follows algorithm by Paul VanRaden integer ntrait,itrno,ntab,ipnts,inorec,n,na,nfix,un,i real*8 r(ntrait,ntrait),g(ntrait,ntrait),trtab(ntab,2), + alfa,ve,vg,vh,ne,tvar,tr,h2 c print 15 write(un,15) do 10 i=1,itrno tvar=g(i,i)+r(i,i) alfa=r(i,i)/g(i,i) ve=r(i,i)*sqrt(2./(n - nfix)) call trinp + (trtab,ntab,ipnts,na,inorec,alfa,tr) ne=alfa/(na/(na-tr*alfa) - na/(na-1.)) vg=2.*(na-1.+(na-1.)**2/(n-nfix))*(r(i,i)/ne/na)**2 vg=vg + 4.*r(i,i)*g(i,i)/ne/na vg=vg + 2.*g(i,i)**2/(na-1.) vh=vg/g(i,i)**2 + (ve**2 + vg)/tvar**2 h2=g(i,i)/tvar vh=sqrt(vh - 2.*vg/g(i,i)/tvar)*h2 vg=sqrt(vg) print 30, i,r(i,i),ve,g(i,i),vg,g(i,i)/(g(i,i)+r(i,i)),vh 10 write(un,30)i,r(i,i),ve,g(i,i),vg,g(i,i)/(g(i,i)+r(i,i)),vh c 15 format(' estimates of variance components (appr. standard errors)' + /' trait',3x,'residual',18x,'genetic',18x, + ' heritability') 30 format(i2,4x,2(g12.4,'(',g12.4,') '),f6.3,'(',f6.3,')') end subroutine quadrmt(uav,a,ia,ja,u,ntrait,itrno,mx,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 c------multiple-trait version integer ia(1),ja(1),ibeg,iend,i,j,k,l ,ntrait,mx,itrno,ivbeg,n real*8 uav(ntrait,ntrait),u(mx,ntrait),a(1) do 5 i=1,itrno do 5 j=1,itrno 5 uav(i,j)=0 do 10 i=ibeg,iend do 10 j=ia(i),ia(i+1)-1 if (ja(j).ge.ibeg .and. ja(j).le.iend) then do 20 k=1,itrno do 20 l=1,itrno 20 uav(k,l)=uav(k,l)+a(j)*u(i+ivbeg,k)*u(ja(j)+ivbeg,l) endif 10 continue end subroutine sortqr(a,n,m,ix,iy,col,icol) c sorts matrix of reals by icol columns given by col, in ascending c order. matrix a is n x m, only first ix rows and iy columns are sorted integer n,m,ix,iy,icol integer i,j,j1,j2,k,l,r,s,col(icol) real*8 a(n,m),x(20),a1,a2 integer stack(50,2) s=1 stack(1,1)=1 stack(1,2)=ix 10 l=stack(s,1) r=stack(s,2) s=s-1 20 i=l j=r j1=(l+r)/2 do 25 k=1,icol 25 x(k)=a(j1,col(k)) 30 do 35 k=1,icol a1=a(i,col(k)) a2=x(k) if (a1.lt.a2) then i=i+1 goto 30 else if(a1.gt.a2) then goto 40 endif 35 continue 40 do 45 k=1,icol a1=a(j,col(k)) a2=x(k) if (a1.gt.a2) then j=j-1 goto 40 else if(a1.lt.a2) then goto 47 endif 45 continue 47 if (i.le.j) then do 50 j2=1,iy a1=a(i,j2) a(i,j2)=a(j,j2) 50 a(j,j2)=a1 i=i+1 j=j-1 endif if (i.le.j) goto 30 if (i.lt.r) then s=s+1 stack(s,1)=i stack(s,2)=r endif r=j if (l.lt.r) goto 20 if (s.ne.0) goto 10 return end subroutine hashmx(y,j1,k1,ind,m,nr) c stores or acculuates a spare matrix element y(j1,k1) in a hash c table ind, which rank is m x 3. nr returns the number of stored elemen c integer j1,k1,m,nr real*8 y,ind(m,3) integer*4 iaddr integer ie,ieq,izer,iaddress,k,j data ie/641/ c iaddr=433*j1+53*k1 iaddress=mod(iabs(iaddr),m)+1 do 10 k=1,200 j=iaddress if (ind(iaddress,1).ne.j1 .or. ind(iaddress,2).ne.k1) then ieq=1 else ieq=0 endif if (ind(iaddress,1).ne.0) then izer=1 else izer=0 endif if (izer.eq.0 .or. ieq.eq.0) then if (izer.eq.0) then ind(iaddress,1)=j1 ind(iaddress,2)=k1 ind(iaddress,3)=y nr=nr+1 else ind(iaddress,3)=ind(iaddress,3)+y endif c if (ind(iaddress,3).eq.0) then c print*,'hashmx: zero accumulation',j1,k1,ind(iaddress,3) c endif return endif iaddress=mod(iaddress+ie-1,m)+1 10 continue write(*,60)nr,m 60 format(' hash matrix too small,filled',i8,' out of',i8) stop end subroutine relhash(ifact,x,nhash,ix,var,stor) integer ifact(4),nhash,ix,ian,isr,idam, itype,stor,storfull, $storhalf real*8 x(nhash,3),type,var parameter (storfull=1,storhalf=2) logical st st=stor.eq.storfull ian=ifact(1) isr=ifact(2) idam=ifact(3) itype=ifact(4) if (itype.eq.1) type=var if (itype.eq.2) type=var*2./3 if (itype.eq.3) type=var*.5 call hashmx(2*type,ian,ian,x,nhash,ix) call hashmx(type/2,idam,idam,x,nhash,ix) call hashmx(type/2,isr,isr,x,nhash,ix) if (st .or. ian.le.isr) call hashmx(-type,ian,isr,x,nhash,ix) if (st .or. ian.le.idam) call hashmx(-type,ian,idam,x,nhash,ix) if (st .or. isr.le.ian) call hashmx(-type,isr,ian,x,nhash,ix) if (st .or. isr.le.idam) call hashmx(type/2,isr,idam,x,nhash,ix) if (st .or. idam.le.ian) call hashmx(-type,idam,ian,x,nhash,ix) if (st .or. idam.le.isr) call hashmx(type/2,idam,isr,x,nhash,ix) return end subroutine hashia(x,nhash,n,ia,ja,a,m,tmp) c copies data form hash to ia-ja-a form, tmp is a temporary array of size n c Entries are sorted within rows integer nhash,n,m,ia(n+1),ja(m),tmp(n) real*8 x(nhash,3),a(m) integer i,imin,imax,j,k,maxcol,maxrow,size c c count the number of entries in each column of a call zero(n,tmp) maxrow=0 maxcol=0 do i=1,nhash j=x(i,1) if (j.ne.0) then if (j.gt.n) then maxrow=max(maxrow,int(x(i,1))) x(i,1)=0 elseif (x(i,2).gt.n) then maxcol=max(maxcol,int(x(i,2))) x(i,1)=0 else tmp(j)=tmp(j)+1 endif endif enddo c create the row count for b ia(1)=1 do i=1,n ia(i+1)=ia(i)+tmp(i) enddo if (ia(n+1)-1.gt.m) then print*,'Too small parameter m in hashia:should be > ', ia(n+1) stop endif c load a into b call zero(n,tmp) do i=1,nhash j=x(i,1) if (j.ne.0) then k=ia(j)+tmp(j) ja(k)=x(i,2) a(k)=x(i,3) tmp(j)=tmp(j)+1 endif enddo c if (maxrow*maxcol.ne.0) then c print*,'HASHIA:declared size=',n,' found columns=',maxcol, c + ' and rows =',maxrow c endif c sort columns do i=1,n imin=ia(i) imax=ia(i+1)-1 size=imax-imin+1 if (size.gt.1) then call sortjsp(ja(imin),a(imin),size) endif enddo end subroutine sortjsp(ja,xa,n) c sorts intger vector ja in ascending order, and orders the real*8 vector a c accordingly integer i,j,s,l,r,j1,a1,n,x integer ja(n),stack(50,2) real*8 xa(n),xb s=1 stack(1,1)=1 stack(1,2)=n 10 l=stack(s,1) r=stack(s,2) s=s-1 20 i=l j=r j1=(l+r)/2 x=ja(j1) 30 if (ja(i).lt.x) then i=i+1 goto 30 endif 40 if (ja(j).gt.x) then j=j-1 goto 40 endif c print*,i,j,ja(i),ja(j) if (i.le.j) then a1=ja(i) ja(i)=ja(j) ja(j)=a1 xb=xa(i) xa(i)=xa(j) xa(j)=xb i=i+1 j=j-1 endif if (i.le.j) goto 30 if (i.lt.r) then s=s+1 stack(s,1)=i stack(s,2)=r endif r=j if (l.lt.r) goto 20 if (s.ne.0) goto 10 end subroutine quadrf(qaq,a,ia,ja,q,ibeg,iend,ivbeg) c calculates the quadratic form qAq of part of vector q and part of c sparse matrix a (a,ia,ja). Quadratic form is taken from row-column c ibeg to iend. Vector q is used from ivbeg+1 location integer ia(1),ja(1),ibeg,iend,ivbeg,i,j real*8 qaq real*8a(1),q(1) qaq=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) $ qaq=qaq+a(j)*q(i+ivbeg)*q(ja(j)+ivbeg) return end subroutine adjadd(a,ia,ja,i1,i2,b,ib,jb,j1,mult,stora,storb) c to sparse matrix a (a,ia,ja), rows j1:j2 add sparse matrix b (b,ib,jb) c beginning from row j1, multiplied by mult.It is assumed that elements c of a that are being added already exist. integer storhalf,storfull parameter(storfull=1,storhalf=2) integer ia(1),ja(1),i1,i2,ib(1),jb(1),j1,stora,storb,i,j,rowb, $ rowb1 real*8 a(1),b(1),mult if (stora.ne.storhalf.or.storb.ne.storfull) then print*,'option not implemented' stop endif do 10 i=i1,i2 rowb=ib(i+j1-i1) rowb1=ib(i+j1-i1+1) do 10 j=ia(i),ia(i+1)-1 20 if (rowb.lt.rowb1) then if (jb(rowb)+i1-j1.eq.ja(j)) then a(j)=a(j)+mult*b(rowb) rowb=rowb+1 else rowb=rowb+1 goto 20 endif endif 10 continue return end subroutine addspar(a,ia,ja,k,l,x) c adds x to (j,k) element of sparse matrix (a,ia,ja). The sparse matrix c must already contain (k,l) location real*8a(1),x integer ia(1),ja(1),k,l,i do 10 i=ia(k),ia(k+1)-1 if (ja(i).eq.l) then a(i)=a(i)+x goto 20 endif 10 continue print*,'addspar:indices not found>>',k,l 20 return end subroutine prhash(h,id,nhash) c prints matrix in the hash table, must be <13 integer nhash,i,j,isize,id real*8 h(id,3),x(12,12) isize=0 do 10 i=1,12 do 10 j=1,12 10 x(i,j)=0 do 15 i=1,nhash if (h(i,1).ne.0) then if (h(i,1).le.12 .and. h(i,2).le.12) then isize=max(int(h(i,1)),int(h(i,2)),isize) x(int(h(i,1)),int(h(i,2)))=h(i,3) else print*,'prhash:elements too large>',(x(i,j),j=1,3) endif endif 15 continue do 30 i=1,isize 30 print '(12f7.2)',(x(i,j),j=1,isize) return end subroutine prset(h,nh,isize) c prints a small sparse matrix of rank isize<=12 stored in hash c table h(nh,3) integer nh,isize real*8 x(12,12),h(nh,3) integer i,j if (isize.gt.12) then print*,'matrix too large to print' stop endif do 10 i=1,isize do 10 j=1,isize 10 x(i,j)=0 do 20 i=1,isize do 20 j=h(i,1),h(i+1,1)-1 if (h(j,2).gt.isize) then print*,'element outside bound',i,h(j,2),h(j,3) else x(i,int(h(j,2)))=h(j,3) endif 20 continue do 30 i=1,isize 30 print '(12f7.2)',(x(i,j),j=1,isize) return end subroutine prspar(a,ia,ja,isize) c prints a small sparse matrix (a,ia,ja) of rank isize<=12 real*8 a(1),x(12,12) integer ia(1),ja(1),isize,i,j if (isize.gt.12) then print*,'matrix too large to print' stop endif do 10 i=1,isize do 10 j=1,isize 10 x(i,j)=0 do 20 i=1,isize do 20 j=ia(i),ia(i+1)-1 if (ja(j).gt.isize) then print*,'element outside bound',i,ja(j),a(j) else if (ja(j).ne.0) x(i,ja(j))=a(j) endif 20 continue do 30 i=1,isize 30 print '(12f7.2)',(x(i,j),j=1,isize) return end subroutine prsparu(a,ia,ju,iju,isize) c prints a small sparse matrix (a,ia,ja) of rank isize<=12 real*8 a(1),x(12,12) integer ia(1),ju(1),iju(1),isize,i,j if (isize.gt.12) then print*,'matrix too large to print' stop endif do 10 i=1,isize do 10 j=1,isize 10 x(i,j)=0 do 20 i=1,isize do 20 j=ia(i),ia(i+1)-1 if (ju(iju(i)-ia(i)+j).gt.isize) then print*,'element outside bound',i,ju(iju(i)-ia(i)+j),a(j) else if (ju(iju(i)-ia(i)+j).ne.0) x(i,ju(iju(i)-ia(i)+j))=a(j) endif 20 continue do 30 i=1,isize 30 print '(12f7.2)',(x(i,j),j=1,isize) return end subroutine prindspar(ia,ja,isize) c prints indices ofs a small sparse matrix (a,ia,ja) of rank isize<=12 integer ia(1),ja(1),isize,i,j,x(12,12) if (isize.gt.12) then print*,'matrix too large to print' stop endif do 10 i=1,isize do 10 j=1,isize 10 x(i,j)=0 do 20 i=1,isize do 20 j=ia(i),ia(i+1)-1 if (ja(j).gt.isize) then print*,'element outside bound',i,ja(j) else if (ja(j).ne.0) x(i,ja(j))=x(i,ja(j))+1 endif 20 continue do 30 i=1,isize 30 print '(12i4)',(x(i,j),j=1,isize) return end subroutine checkdiag(n,text,ia,ja,a) c quits if sparse matrix ia-ja-a of rank n has missing or <=0 diagonals integer n,ia(1),ja(1),i,j,inddiag real*8 a(1) character text*(*) c do i=1,n inddiag=0 if (ia(i+1)-ia(i).le.0) then print '('' equation '',i6,'' in '',a,'' matrix empty'')', + i,text stop endif do j=ia(i), ia(i+1)-1 if (ja(j).eq.i) inddiag=j enddo if (inddiag.eq.0) then print '('' missing diagonal '',i6,'' in '',a,'' matrix'')',i stop elseif(a(inddiag).le.0) then print '('' <=0 diagonal '',i6,'' in '',a,'' matrix'')',i stop endif enddo end subroutine fspakerr(iflag) c error display routine for fspak integer iflag if (iflag.eq.0) then return elseif (iflag.gt.100 .and. iflag.lt.100000) then print*,'Insufficient memory in FSPAK for option:',iflag/100 elseif (iflag.gt.4000000) then print*,'not positive-definite matrix' else print*,'undetermined problem in FSPAK, flag=',iflag endif stop end function tracediag(n,ia,ja,a,l1,l2) c computes the sum of diagonal elements l1-l2 in sparse matrix ia-ja-a integer n,ia(1),ja(1),l1,l2,i,j,msglev real*8 tracediag,a(1) data msglev/0/ tracediag=0 if (l1.le.0 .or. l2.gt.n) then print*,' l1 or l2 out of range. l1.l2.n:',l1,l2,n stop endif do i=l1,l2 do j=ia(i),ia(i+1)-1 if (ja(j).eq.i) then tracediag=tracediag+a(j) goto 10 endif enddo print*,'tracediag: missing diagonal element',i stop 10 continue enddo if (msglev.ge.1) then print*,'TRACEDIAG: input matrix' call prspar(a,ia,ja,n) print*,'tracediag,l1,l2',tracediag,l1,l2 endif end function tracematrix(n,ia,ja,a,n1,ib,jb,b,ix,iy,tmp) c calculates trace of sparse matrices a-ia-ja of order c n and ib-jb-b of order n1. Columns of ib-jb-b over n1 are ignored. c ix and iy specify offsets for the second matrix. c tmp is a real(*8) matrix of order n. c c The matrices may be full or half stored. In the later case, c the trace is computed as a sum of trace due to diagonal elements c plus 2 * sum of trace due to off-diagonals c integer n,ia(1),ja(1),n1,ib(1),jb(1),ix,iy,i,j,msglev real*8 tracematrix,a(1),b(1),tmp(1),tl,td,tu data msglev/0/ do i=1,n tmp(i)=0 enddo tl=0 td=0 tu=0 c scatter row i of b do i=1,n1 do j=ib(i),ib(i+1)-1 if (jb(j).le.n1) then tmp(jb(j)+iy)=b(j) endif enddo c actual trace computation do j=ia(i+ix),ia(i+ix+1)-1 if (ja(j).gt.i+ix) then tl=tl+a(j)*tmp(ja(j)) elseif (ja(j).eq.i+ix) then td=td+a(j)*tmp(ja(j)) else tu=tu+a(j)*tmp(ja(j)) endif enddo c zero row in tmp do j=ib(i+ix),ib(i+ix+1)-1 if (jb(j).le.n1) then tmp(jb(j)+iy)=0 endif enddo enddo c Matrices full stored if upper and lower traces equal if (tu.eq.0) then tracematrix=td+2*tl elseif (tl.eq.0) then tracematrix=td+2*tu else if (abs((tl-tu)/(tl+tu)).lt.1e-6) then tracematrix=td+2*tu else print*,'tracematrix:tl,td,tu=',tl,td,tu stop endif if (msglev.ge.1) then print*,'tracematrix:tl,td,tu=',tl,td,tu print*,'TRACEMATRIX: matrix A' call prspar(a,ia,ja,n) print*,'matrix B' call prspar(b,ib,jb,n1) print*,'tracematrix,ix,iy',tracematrix,ix,iy endif end C 24 OCT 1993 / IBM c ---------------- subroutine fspak c ---------------- + (option,n,ia,ja,a,b,flag,iout,ioor, + available,needed,is,fnode,bnode, + fnseg,bnseg,fx,feqb) c ---------------------------------------------------- c Public Sparse Matrix Package for Symmetric Matrices c Authors: Miguel Perez-Enciso, Ignacy Misztal c Mon Aug 24, 1992 12:39:46 - Oct 24, 1993 c ---------------------------------------------------- implicit none integer zero,order,symfac,numfac,solve,sinit,ssolve + ,det,ldet,check,spin,spin1,restart,stat,equal,mx_st parameter (order=10, symfac=20, numfac=40, solve=50, + sinit=51, ssolve=52, det=54, ldet =55, + check=56, spin=60, spin1=61, restart=70, + stat=80, equal=1, mx_st=50 ) integer available,r_avail + ,delta,flag,ioor,iout,n,ia(1),ja(1),is(available) + ,fnode(1),bnode(1),fnseg,bnseg,feqb,ratio,fratio + ,path,found,fpath,fstart,fpnode,bpath,bstart,bpnode + ,maxju,maxu,needed,maxneed,i,adjncy,d,invp,iu,iju,ju,u + ,maxint,nza,nzad,option,perm,iap,jap,ap,xadj,first,ios + ,zja,za,zd,utia,utja,work1,work2,work3,work4,work5,work6 real*8 a(1),b(1),fx(1),st(mx_st),t_t,t0,second,tol equivalence (st(1),maxju),(st(2),maxu),(st(3),maxneed),(st(4),nza) c this is to avoid that the value of ratio is lost c with dynamic allocation memory save ratio data st/mx_st*0./,first/0/ c.. tolerance for accuracy in solution (in option check) + ,tol/1.d-6/ flag=0 t_t=second() c ----------------------------------------------- c first time it checks whether ratio equal 1 or 2 c ----------------------------------------------- if (first.eq.0) then ratio=fratio() first=1 endif c ------------------------------------------- c no. of off-diagonal nze in upper triangular c ------------------------------------------- if(option.le.numfac) nza=(ia(n+1)-1-n) c -------------------- c no. of nze in adjncy c -------------------- nzad=2*nza c ---------------- c allocate storage c ---------------- perm = 1 invp = perm + n c ------------------------- if (option.eq.order) then c ------------------------- xadj = invp + n adjncy = xadj + n+1 work1 = adjncy + nzad work2 = work1 + n work3 = work2 + n work4 = work3 + n needed = work4 + n call checkst (option,needed,maxneed,available,iout,flag) if (flag.ne.0) return c -------------------------------------------------- c max integer (used as flag in the ordering routine) c -------------------------------------------------- maxint=9999999 c --------------------------------------------------------- c parameter for type of ordering -1=MD, 0=MMD, >0=slack MMD c --------------------------------------------------------- delta=0 call unfold + (n,ia,ja,is(xadj),is(adjncy),is(work1)) t0=second() c ------------------------------- c if mmd ordering routines by Liu c ------------------------------- call genmmd + (n,is(xadj),is(adjncy),is(invp),is(perm),delta, + is(work1),is(work2),is(work3),is(work4),maxint, + maxju) c --------------------------- c elseif md sparspak routines c --------------------------- c work5 = work4 + n c work6 = work5 + n c needed= work6 + n c call checkst (option,needed,maxneed,available,iout,flag) c if(flag.ne.0) return c call genqmd c + (n,is(xadj),is(adjncy),is(perm),is(invp), c + is(work1),is(work2),is(work3),is(work4), c + is(work5),is(work6),maxju) st(21)=st(21)+second()-t0 c ------------------------------ c write down permutation vectors c ------------------------------ rewind ioor write(ioor,*) n,maxju write(ioor,'(10i7)') (is(i),i=perm,perm+n-1), + (is(i),i=invp,invp+n-1) c ------------------------------- elseif (option.eq.restart) then c ------------------------------- rewind ioor i=0 read(ioor,*,end=1) i,maxju 1 if(i.eq.n) then read(ioor,'(10i7)',end=2) (is(i),i=perm,perm+n-1), + (is(i),i=invp,invp+n-1) 2 inquire(ioor,iostat=ios) if (ios.ne.0) then write(iout,*) 'error in unit IOOR ' flag=-2 endif else if (i.ne.0) then write(iout,*) 'incorrect size of permutation vector' write(iout,*) 'size ',i,' found; size ',n, ' required' endif flag=-1 endif c ------------------------------ elseif (option.eq.symfac) then c ------------------------------ iu = invp + n iju = iu + n+1 ju = iju + n+1 xadj = ju + maxju adjncy = xadj + n+1 work1 = adjncy + nzad work2 = work1 + n work3 = work2 + n needed = work3 + n call checkst (option,needed,maxneed,available,iout,flag) if (flag.ne.0) return t0=second() call unfold + (n,ia,ja,is(xadj),is(adjncy),is(work1)) call smbfct + (n,is(xadj),is(adjncy),is(perm),is(invp),is(iu), + maxu,is(iju),is(ju),maxju,is(work1),is(work2), + is(work3),flag) st(22)=st(22)+second()-t0 else iu = invp + n iju = iu + n+1 ju = iju + n+1 d = ju + maxju u = d + n*ratio c -------------------------- if (option.eq.numfac) then c -------------------------- iap = u + maxu*ratio jap = iap + n+1 ap = jap + nza work1 = ap + nza*ratio needed = work1 + n*ratio call checkst (option,needed,maxneed,available,iout,flag) if (flag.ne.0) return t0=second() c --------------------------------------- c first load numerical values of a into u c --------------------------------------- call loadap + (n,ia,ja,a,is(iap),is(jap),is(ap),is(d),is(invp), + is(work1)) call loadu + (n,is(iap),is(jap),is(ap),is(iu),is(ju),is(iju), + is(u),is(work1)) work1 = u + maxu*ratio work2 = work1 + n work3 = work2 + n needed = work3 + n*ratio call checkst (option,needed,maxneed,available,iout,flag) if (flag.ne.0) return c ----------------------- c numerical factorization (slightly modified sparspak routine) c ---------------------------------- call fgsfct + (n,is(iu),is(u),is(iju),is(ju),is(d),is(work1), + is(work2),is(work3),flag) call merlin (n,is(iu),is(u),is(d)) st(24)=st(24)+second()-t0 st(11)=st(11)+1 c ---------- c zero pivot c ---------- if(flag.ne.0) flag=numfac*1000000 + flag c ----------------------------- elseif (option.eq.solve) then c ----------------------------- work1 = u + maxu*ratio needed = work1 + n*ratio t0=second() call permute (n,b,is(work1),is(invp)) call fullfb + (n,is(iu),is(ju),is(iju),is(u),is(d),is(work1)) call permute (n,is(work1),b,is(perm)) st(25)=st(25)+second()-t0 st(12)=st(12)+1 c -------------------------------------------------- elseif (option.ge.sinit.and.option.le.ssolve) then c -------------------------------------------------- path = u + maxu*ratio found = path + n work1 = found + n fstart= work1 + n*ratio fpnode= fstart+ fnseg+1 needed= fpnode+ fnseg call checkst c ------------------------------------------ c allows a safety margin for fpath and bpath c ------------------------------------------ + (option,needed+200,maxneed,available,iout, + flag) if(flag.ne.0) return c ------------------------- if (option.eq.sinit) then c ------------------------- c initialize and computes factorization table c ------------------------------------------- call rzero (n,is(work1)) call lzero (n,is(found)) call makepath (n,is(path),is(iu),is(ju),is(iju)) c ------------------------------ elseif (option.eq.ssolve) then c ------------------------------ t0=second() fpath = fpnode + fnseg c ------------------------------------------------------ c obtains forward path (in the permuted vector of nodes) c ------------------------------------------------------ call spermute (fnseg,fnode,is(fpnode),is(invp)) call segpath + (n,is(fpnode),fnseg,is(fstart),is(fpath), + is(path),is(found)) st(34)=st(34)+is(fstart+fnseg)-1 if(feqb.eq.equal) then needed= fpath + is(fstart+fnseg)-1 call checkst + (option,needed,maxneed,available,iout, + flag) if (flag.ne.0) return c ----------------------------- c initialize working rhs vector c ----------------------------- call zerout (is(work1),is(fpath),is(fstart),fnseg) c ---------------- c input rhs values c ---------------- call expand (fnseg,is(fpnode),fx,is(work1)) c ----- c solve c ----- call fastfb + (is(fpath),is(fstart),fnseg,is(fpath), + is(fstart),fnseg,is(iu),is(ju),is(iju), + is(u),is(d),is(work1)) c --------------------- c stores solutions in b c --------------------- call compress (fnseg,is(fpnode),b,is(work1)) st(31)= max0(int(st(31)),is(fstart+fnseg)-1) st(35)=st(35)+is(fstart+fnseg)-1 else bstart = fpath + is(fstart+fnseg)-1 bpnode = bstart+ bnseg+1 bpath = bpnode+ bnseg c ---------------------------------- c obtains backward path if necessary c ---------------------------------- call spermute (bnseg,bnode,is(bpnode),is(invp)) call segpath + (n,is(bpnode),bnseg,is(bstart), + is(bpath),is(path),is(found)) needed= bpath + is(bstart+bnseg)-1 call checkst + (option,needed,maxneed,available,iout, + flag) if (flag.ne.0) return call zerout (is(work1),is(bpath),is(bstart),bnseg) call zerout (is(work1),is(fpath),is(fstart),fnseg) call expand (fnseg,is(fpnode),fx,is(work1)) call fastfb + (is(fpath),is(fstart),fnseg,is(bpath), + is(bstart),bnseg,is(iu),is(ju),is(iju), + is(u),is(d),is(work1)) call compress (bnseg,is(bpnode),b,is(work1)) st(31)=max(int(st(31)),is(bstart+bnseg)-1) st(35)=st(35)+is(bstart+bnseg)-1 endif st(13)=st(13)+1 st(27)=st(27)+second()-t0 st(30)=max(int(st(30)),fnseg) st(31)=max(int(st(31)),bnseg) st(32)=st(32)+fnseg st(33)=st(33)+bnseg endif c --------------------------- elseif (option.eq.det) then c --------------------------- call mkdet (n,is(d),b) st(14)=st(14)+1 c ---------------------------- elseif (option.eq.ldet) then c ---------------------------- call mkldet (n,is(d),b) st(14)=st(14)+1 c ---------------------------- elseif (option.eq.check) then c ---------------------------- work1 = u + maxu*ratio call checkb (n,ia,ja,a,b,fx,tol,is(work1),iout,flag) c ---------------------------- elseif (option.eq.spin) then c ---------------------------- zja = u + maxu*ratio za = zja + maxu+n zd = za + (maxu+n)*ratio utia = zd + n*ratio utja = utia + n+1 work1 = utja + maxu+n work2 = work1 + n work3 = work2 + n needed = work3 + n*ratio call checkst (option,needed,maxneed,available,iout,flag) if (flag.ne.0) return t0=second() call fsinverse + (n,is(perm),is(invp),is(d),is(iju),is(ju), + is(iu),is(u),is(work3),is(za),is(zd), + is(zja),is(utia),is(utja), + is(work1),is(work2),ia,ja,a) st(26)=st(26)+second()-t0 st(15)=st(15)+1 c ---------------------------- elseif (option.eq.spin1) then c ---------------------------- work1 = u + maxu*ratio work2 = work1 + n*ratio work5 = work2 + n*ratio needed = work5 + maxu call checkst (option,needed,maxneed,available,iout,flag) if (flag.ne.0) return t0=second() call fsinverse1 + (n,is(perm),is(invp),is(d),is(iju),is(ju), + is(iu),is(u),is(work1),is(work2),is(work1), + is(work2),is(work5),ia,ja,a) st(26)=st(26)+second()-t0 st(15)=st(15)+1 c ---------------------------- elseif (option.eq.stat) then c ---------------------------- call figureout (n,available,iout,st) c ---- else c ---- write(iout,*) 'option no. ',option,' not implemented' flag=-99 endif endif st(20)=st(20)+second()-t_t return end c last modified 12/04/93 c c list of subroutines needed for FSPAK c Miguel PEREZ-ENCISO c Ignacy MISZTAL c c second c checkb c checkst c fastfb c fullfb c fgsfct c figureout c fratio c fsinverse c ival c loadap c loadu c lzero c makepath c mkdet c mkldet c permute c rowtocol c rowtocolu c rzero c segpath c spermute c unfold c zero c zerout c ----------------- subroutine checkb (n,ia,ja,a,b,x,tol,x1,iout,flag) c ----------------- c mpe, 18 Dec 1992 implicit none integer n,ia(1),ja(1),iout,flag,i,j,k real*8 dev,tol,a(1),b(1),x(1),x1(1) dev=0. do i=1,n x1(i)=0. enddo do i=1,n do k=ia(i),ia(i+1)-1 j=ja(k) x1(i)=x1(i)+a(k)*b(j) if(i.lt.j) x1(j)=x1(j)+a(k)*b(i) enddo enddo do i=1,n dev=dev+abs(x1(i)-x(i))/n enddo if(dev.gt.tol) then flag=1 write(iout,*) 'ERROR in CHECKB: not accurate solution' write(iout,*) 'observed deviation = ',dev write(iout,*) 'maximum allowed = ',tol endif return end c ------------------ subroutine checkst (option,needed,maxneed,available,iout,flag) c ------------------ integer option,needed,maxneed,available,iout,flag maxneed=max(maxneed,needed) if(maxneed.gt.available) then write(iout,*) ' insufficient storage in option ',option write(iout,*) ' avaliable ',available,'; needed ',needed flag=option*100+1 endif return end c ------------------- subroutine compress (nc,ipos,bc,b) c ------------------- c compreses nc positions of b into bc integer ipos(1) real*8 bc(1),b(1) do i=1,nc bc(i)=b(ipos(i)) enddo return end c ----------------- subroutine expand (nc,ipos,bc,b) c ----------------- c expand a compressed vector bc into b c MPE integer ipos(1) real*8 bc(1),b(1) do i=1,nc b(ipos(i))=bc(i) enddo return end c ----------------- subroutine fastfb c ----------------- + (fpath,fstart,fnseg,bpath,bstart,bnseg, + iu,ju,iju,u,d,b) c performs fast forward/backward on a triangular system iuju,iju,u c M_P-E Thu Aug 13, 1992 14:34:41 implicit none integer fpath(1),fstart(1),fnseg,bpath(1),bstart(1),bnseg, + iu(1),ju(1),iju(1),i,j,k,ik,iseg real*8 u(1),d(1),b(1),bi c fast forward in order within segments with reverse order btw segments do iseg=fnseg,1,-1 do ik=fstart(iseg),fstart(iseg+1)-1 i=fpath(ik) bi=b(i) do k=iu(i),iu(i+1)-1 j=ju(iju(i)+k-iu(i)) b(j)=b(j)+u(k)*bi enddo enddo c fast diagonal do ik=fstart(iseg),fstart(iseg+1)-1 i=fpath(ik) b(i)=b(i)*d(i) enddo enddo c fast backward in reverse order within segments with reverse order btw seg c segm do iseg=1,bnseg do ik=bstart(iseg+1)-1,bstart(iseg),-1 i=bpath(ik) bi=b(i) do k=iu(i),iu(i+1)-1 j=ju(iju(i)+k-iu(i)) bi=bi+u(k)*b(j) enddo b(i)=bi enddo enddo return end c ----------------- subroutine fullfb (n,iu,ju,iju,u,d,b) c ----------------- c performs full forward/backward on a triangular system iuju,iju,u c M_P-E Sat Jul 11, 1992 17:54:59 implicit none integer n,iu(1),ju(1),iju(1),i,j,k real*8 u(1),d(1),b(1),bi c full forward do i=1,n bi=b(i) do k=iu(i),iu(i+1)-1 j=ju(iju(i)+k-iu(i)) b(j)=b(j)+u(k)*bi enddo enddo c diagonal do i=1,n b(i)=b(i)*d(i) enddo c full backward do i=n,1,-1 bi=b(i) do k=iu(i),iu(i+1)-1 j=ju(iju(i)+k-iu(i)) bi=bi+u(k)*b(j) enddo b(i)=bi enddo return end c --------------- function fratio() c --------------- c returns the ratio of storage occupied by reals and integers; c the only two values returned are 1 and 2. Orginally written to c facilitate easy porting between Crays (ratio=1) to other c computers (ratio=2, if reals are real*8) c Ignacy Misztal integer fratio,ix(2),const parameter(const=12345) real*8 y equivalence (y,ix(1)) ix(2)=const y=0 if (ix(2).eq.const) then fratio=1 else fratio=2 endif end C Slightly modified and restructured C M_P-E Sat Jul 11, 1992 16:50:38 C From George & Liu (1981) C----- SUBROUTINE FGSFCT C*************************************************************** C*************************************************************** C****** GSFCT ..... GENERAL SPARSE SYMMETRIC FACT ****** C*************************************************************** C*************************************************************** C C PURPOSE - THIS SUBROUTINE PERFORMS THE SYMMETRIC C FACTORIZATION FOR A GENERAL SPARSE SYSTEM, STORED IN C THE COMPRESSED SUBSCRIPT DATA FORMAT. C 1 C INPUT PARAMETERS - 1 C NEQNS - NUMBER OF EQUATIONS. 1 C XLNZ - INDEX VECTOR FOR LNZ. XLNZ(I) POINTS TO THE 1 C START OF NONZEROS IN COLUMN I OF FACTOR L. 1 C (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT DATA 1 C STRUCTURE FOR FACTOR L. 1 C 1 C UPDATED PARAMETERS - 1 C LNZ - ON INPUT, CONTAINS NONZEROS OF A, AND ON 1 C RETURN, THE NONZEROS OF L. 2 C DIAG - THE DIAGONAL OF L OVERWRITES THAT OF A. 2 C IFLAG - THE ERROR FLAG. IT IS SET TO 1 IF A ZERO OR 2 C NEGATIVE SQUARE ROOT OCCURS DURING THE 2 C FACTORIZATION. 2 C OPS - A DOUBLE PRECISION COMMON PARAMETER THAT IS 2 C INCREMENTED BY THE NUMBER OF OPERATIONS 2 C PERFORMED BY THE SUBROUTINE. 2 C 2 C WORKING PARAMETERS - 2 C LINK - AT STEP J, THE LIST IN 3 C LINK(J), LINK(LINK(J)), ........... 3 C CONSISTS OF THOSE COLUMNS THAT WILL MODIFY 3 C THE COLUMN L(*,J). 3 C FIRST - TEMPORARY VECTOR TO POINT TO THE FIRST 3 C NONZERO IN EACH COLUMN THAT WILL BE USED 3 C NEXT FOR MODIFICATION. 3 C TEMP - A TEMPORARY VECTOR TO ACCUMULATE MODIFICATIONS. 3 C 3 C*************************************************************** 3 C 4 SUBROUTINE FGSFCT ( NEQNS, XLNZ, LNZ, XNZSUB, NZSUB, DIAG, 1 LINK, FIRST, TEMP, IFLAG ) 4 C 4 C*************************************************************** 4 C 4 REAL*8 COUNT, OPS COMMON /SPKOPS/ OPS 4 REAL*8 DIAG(1), LNZ(1), TEMP(1), DIAGJ, LJK INTEGER LINK(1), NZSUB(1) 4 INTEGER FIRST(1), XLNZ(1), XNZSUB(1), 5 1 I, IFLAG, II, ISTOP, ISTRT, ISUB, J, 5 1 K, KFIRST, NEQNS, NEWK 5 C 5 C*************************************************************** 5 C 5 C ------------------------------ 5 C INITIALIZE WORKING VECTORS ... 5 C ------------------------------ 5 DO I = 1, NEQNS LINK(I) = 0 6 TEMP(I) = 0.0 6 ENDDO C -------------------------------------------- 6 C COMPUTE COLUMN L(*,J) FOR J = 1,...., NEQNS. 6 C -------------------------------------------- 6 DO J = 1, NEQNS C ------------------------------------------- 6 C FOR EACH COLUMN L(*,K) THAT AFFECTS L(*,J). 6 C ------------------------------------------- 6 DIAGJ = 0.0 7 NEWK = LINK(J) 7 K = NEWK DO WHILE (K.NE.0) NEWK = LINK(K) 7 C --------------------------------------- 7 C OUTER PRODUCT MODIFICATION OF L(*,J) BY 7 C L(*,K) STARTING AT FIRST(K) OF L(*,K). 7 C --------------------------------------- 7 KFIRST = FIRST(K) 7 LJK = LNZ(KFIRST) 8 DIAGJ = DIAGJ + LJK*LJK 8 OPS = OPS + 1.00 8 ISTRT = KFIRST + 1 8 ISTOP = XLNZ(K+1) - 1 8 IF ( ISTOP .GE. ISTRT ) THEN C ------------------------------------------ 8 C BEFORE MODIFICATION, UPDATE VECTORS FIRST, 8 C AND LINK FOR FUTURE MODIFICATION STEPS. 8 C ------------------------------------------ 8 FIRST(K) = ISTRT 9 I = XNZSUB(K) + (KFIRST-XLNZ(K)) + 1 9 ISUB = NZSUB(I) 9 LINK(K) = LINK(ISUB) 9 LINK(ISUB) = K 9 C --------------------------------------- 9 C THE ACTUAL MOD IS SAVED IN VECTOR TEMP. 9 C --------------------------------------- 9 cdir$ ivdep DO II = ISTRT, ISTOP 9 ISUB = NZSUB(I) 9 TEMP(ISUB) = TEMP(ISUB) + LNZ(II)*LJK 10 I = I + 1 10 ENDDO 10 COUNT = ISTOP - ISTRT + 1 10 OPS = OPS + COUNT ENDIF K = NEWK ENDDO C ---------------------------------------------- 10 C APPLY THE MODIFICATIONS ACCUMULATED IN TEMP TO 10 C COLUMN L(*,J). 10 C ---------------------------------------------- 10 DIAGJ = DIAG(J) - DIAGJ 11 IF ( DIAGJ .LE. 0.0E0 ) THEN C ------------------------------------------------------ C ERROR - ZERO OR NEGATIVE SQUARE ROOT IN FACTORIZATION. C ------------------------------------------------------ IFLAG = J RETURN ENDIF DIAGJ = SQRT(DIAGJ) DIAG(J) = DIAGJ ISTRT = XLNZ(J) 11 ISTOP = XLNZ(J+1) - 1 11 IF ( ISTOP .GE. ISTRT ) THEN FIRST(J) = ISTRT 11 I = XNZSUB(J) 11 ISUB = NZSUB(I) 11 LINK(J) = LINK(ISUB) 12 LINK(ISUB) = J 12 DO II = ISTRT, ISTOP ISUB = NZSUB(I) 12 LNZ(II) = ( LNZ(II)-TEMP(ISUB) ) / DIAGJ 12 TEMP(ISUB) = 0.0E0 12 I = I + 1 12 ENDDO COUNT = ISTOP - ISTRT + 1 12 OPS = OPS + COUNT ENDIF ENDDO RETURN END 13 c -------------------- subroutine figureout (n,avail,iout,st) c -------------------- c mpe integer n,avail real*8 st(*) write(iout,'(a,f25.6)')' *************' write(iout,'(a,f25.6)')' *** FSPAK ***' write(iout,'(a,f25.6)')' *************' write(iout,'(a,f25.6)')' MPE / IM ' write(iout,'(a,f25.6)')' August 1992' write(iout,'(a,f25.6)') write(iout,'(a,f25.6)')' SPARSE STATISTICS' write(iout,'(a,I25)') ' RANK OF COEFF MATRIX =',n write(iout,'(a,I25)') ' STORAGE AVAILABLE =',avail write(iout,'(a,I25)') ' MAXIMUM NEEDED =' + ,IVAL(ST(3)) write(iout,'(a,I25)') ' NZE IN UPPER TRIANGULAR =' + ,IVAL(ST(4))+N write(iout,'(a,I25)') ' NZE IN FACTOR =' + ,IVAL(ST(2)) write(iout,'(a,i25)') ' NO. OF CALLS NUM FACT =' + ,INT(ST(11)) write(iout,'(a,i25)') ' NO. OF CALLS SOLVE =' + ,INT(ST(12)) write(iout,'(a,i25)') ' NO. OF CALLS SPARS SOLV =' + ,INT(ST(13)) write(iout,'(a,i25)') ' NO. OF CALLS DET / LDET =' + ,INT(ST(14)) write(iout,'(a,i25)') ' NO. OF CALLS SPARS INV =' + ,INT(ST(15)) if(st(13).ne.0) then write(iout,'(a,i25)') ' MAX NO. NODES =' + ,INT(ST(30)) write(iout,'(a,i25)') ' MAX PATH LENGTH =' + ,INT(ST(31)) write(iout,'(a,f25.3)')' AVG NO. NODES FPATH =' + ,ST(32)/ST(13) write(iout,'(a,f25.3)')' AVG NO. NODES BPATH =' + ,ST(33)/ST(13) write(iout,'(a,f25.3)')' AVG FPATH LENGTH =' + ,ST(34)/ST(13) write(iout,'(a,f25.3)')' AVG BPATH LENGTH =' + ,ST(35)/ST(13) endif write(iout,'(a,f25.6)')' TOTAL CPU TIME IN FSPAK =',st(20) write(iout,'(a,f25.6)')' TIME FOR FINDING ORDER =',st(21) write(iout,'(a,f25.6)')' TIME FOR SYMBOLIC FAC =',st(22) write(iout,'(a,f25.6)')' TIME FOR NUMERICAL FAC =',st(24) write(iout,'(a,f25.6)')' TIME FOR SOLVE =',st(25) write(iout,'(a,f25.6)')' TIME FOR SPARSE SOLVE =',st(27) write(iout,'(a,f25.6)')' TIME FOR SPARSE INVERSE =',st(26) return end c ------------- function ival (st) c ------------- integer st ival=st return end c ----------------- subroutine loadap c ----------------- + (n,ia,ja,a,iap,jap,ap,d,ip,tmp) c this subroutine loads half sparse stored elements of ia,ja,a into half sp c sparse stored c structure iap,jap,ap for non-diagonal elements and into d for diagonal el C elements; c elements are reordered according to permutation vectors p,ip c M_P-E Tue Jun 30, 1992 10:23:01 implicit none integer n,ia(1),ja(1),iap(1),jap(1),ip(1),tmp(1),i,j,k,kk real*8 a(1),ap(1),d(1) do i=1,n tmp(i)=0 enddo c compute no. of entries per upper stored row do k=1,n i=ip(k) do kk=ia(k),ia(k+1)-1 j=ip(ja(kk)) if(i.eq.j .and. i.le.n) then d(i)=a(kk) elseif(i.gt.j .and. i.le.n) then tmp(j)=tmp(j)+1 elseif(i.lt.j .and. j.le.n) then tmp(i)=tmp(i)+1 endif enddo enddo c fill iap iap(1)=1 do i=1,n iap(i+1)=iap(i)+tmp(i) tmp(i)=0 enddo c fill jap,ap do k=1,n i=ip(k) do kk=ia(k),ia(k+1)-1 j=ip(ja(kk)) if(i.gt.j .and. i.le.n) then jap(iap(j)+tmp(j))=i ap(iap(j)+tmp(j))=a(kk) tmp(j)=tmp(j)+1 elseif(i.lt.j .and. i.le.n) then jap(iap(i)+tmp(i))=j ap(iap(i)+tmp(i))=a(kk) tmp(i)=tmp(i)+1 endif enddo enddo return end c ---------------- subroutine loadu (n,iap,jap,ap,iu,ju,iju,u,tmp) c ---------------- c sparse stored matrix iap,jap,ap is copied into sparse compressed storage C GE MATRIX IU,JU,IJU,U c where u allows storage for fill-ins c M_P-E Tue Jun 30, 1992 11:01:23 implicit none integer n,iap(1),jap(1),iu(1),ju(1),iju(1),i,j,k real*8 ap(1),u(1),tmp(1) do i=1,n tmp(i)=0. enddo do i=1,n c - - scatters do k=iap(i),iap(i+1)-1 j=jap(k) tmp(j)=ap(k) enddo c - - loads u and zero out tmp do k=iu(i),iu(i+1)-1 j=ju(iju(i)+k-iu(i)) u(k)=tmp(j) tmp(j)=0. enddo enddo return end c ---------------- subroutine lzero (n,l) c ---------------- c initialize a logical vector c MPE logical l(1) do i=1,n l(i)=.false. enddo return end c ------------------- subroutine makepath c ------------------- + (n,path,iu,ju,iju) c obtains factorization path stored as a linked list where path(i) c is next node of the path of element i c M. Perez-Enciso, Madison,Thu May 21,1992 integer n,path(n),iu(n+1),ju(1),iju(n) do i=1,n path(i)=0 if(iu(i+1).gt.iu(i)) path(i)=ju(iju(i)) enddo return end c ----------------- subroutine merlin (n,iu,u,d) c ----------------- c MPE integer n,iu(1) real*8 u(1),d(1) do i=1,n do j=iu(i),iu(i+1)-1 u(j)=-u(j)/d(i) enddo d(i)=1./(d(i)*d(i)) enddo return end c ---------------- subroutine mkdet (n,d,b) c ---------------- c computes determinant, M_P-E integer i,n real*8 b,d(1) b=1. do i=1,n b=b/d(i) enddo return end c ----------------- subroutine mkldet (n,d,b) c ----------------- c computes log determinant c M_P-E Mon Aug 10, 1992 15:43:45 integer i,n real*8 b,d(1) b=0. do i=1,n b=b+log(1./d(i)) enddo return end c ------------------ subroutine permute (n,vec1,vec2,p) c ------------------ c permutes vector vec1 according to permutation vector p into vec2 c M_P-E Wed Jul 1, 1992 11:12:01 integer n,i,p(1) real*8 vec1(1),vec2(1) do i=1,n vec2(p(i))=vec1(i) enddo return end c ------------------- subroutine rowtocol(aia,aja,aa,bia,bja,ba,tmp,n) c ------------------- c Matrix a, in sparse i-j-a form and stored row-wise, is converted to b, c stored column-wise; after the move, the rows of b are sorted in c increasing order. c a and b are have n rows and columns, and tmp is c is a temporary integer vector of size n. c I. Misztal integer n,aia(1),aja(1),bia(1),bja(1),tmp(1),i,j,k real*8 aa(1),ba(1) c count the number of entries in each row of a call zero(n,tmp) do i=1,n do j=aia(i),aia(i+1)-1 tmp(aja(j))=tmp(aja(j))+1 enddo enddo c create the row count for b bia(1)=1 do i=1,n bia(i+1)=bia(i)+tmp(i) enddo c load a into b call zero(n,tmp) do i=1,n do j=aia(i),aia(i+1)-1 k=bia(aja(j))+tmp(aja(j)) bja(k)=i ba(k)=aa(j) tmp(aja(j))=tmp(aja(j))+1 enddo enddo end c -------------------- subroutine rowtocolu (aia,aju,aiju,aa,bia,bja,ba,tmp,n) c -------------------- c Matrix a, in sparse i-ju-iju -a form and stored row-wise, is converted to C TO B, c stored column-wise; after the move, the rows of b are sorted in c increasing order. c a and b are have n rows and columns, and tmp is c is a temporary integer vector of size n. c I. Misztal integer n,aia(1),bia(1),bja(1),tmp(1),i,j,k, * aju(1),aiju(1) real*8 aa(1),ba(1) c count the number of entries in each row of a call zero(n,tmp) do i=1,n do j=aia(i),aia(i+1)-1 tmp(aju(aiju(i)-aia(i)+j))=tmp(aju(aiju(i)-aia(i)+j))+1 enddo enddo c create the row count for b bia(1)=1 do i=1,n bia(i+1)=bia(i)+tmp(i) enddo c load a into b call zero(n,tmp) do i=1,n do j=aia(i),aia(i+1)-1 k=bia(aju(aiju(i)-aia(i)+j))+tmp(aju(aiju(i)-aia(i)+j)) bja(k)=i ba(k)=aa(j) tmp(aju(aiju(i)-aia(i)+j))=tmp(aju(aiju(i)-aia(i)+j))+1 enddo enddo end c ---------------- subroutine rzero (n,b) c ---------------- c initialize a real vector c MPE real*8 b(1) do i=1,n b(i)=0. enddo return end c ------------------ subroutine segpath c ------------------ + (n,node,nseg,sstart,spath,path,found) c computes a segmented path for the nseg nodes stored in node c M. Perez-Enciso, Madison, Thu May 21, 1992 implicit none integer n,nseg, node(nseg),sstart(nseg+1),spath(1),path(n) + ,ik,k,iseg logical found(n) ik=0 do iseg=1,nseg k=node(iseg) sstart(iseg)=ik+1 do while(k.ne.0.and..not.found(k)) ik=ik+1 spath(ik)=k found(k)=.true. k=path(k) enddo enddo sstart(nseg+1)=ik+1 c zero out boolean working vector do ik=1,sstart(nseg+1)-1 found(spath(ik))=.false. enddo return end c ------------------- subroutine spermute (n,vect1,vect2,perm) c ------------------- c permutes an integer compressed vector c MPE integer vect1(1),vect2(1),perm(1) do i=1,n vect2(i)=perm(vect1(i)) enddo return end c ----------------- subroutine unfold (n,uia,uja,ia,ja,tmp) c ----------------- c this subroutine copies a half sparse stored matrix into a sparse c stored matrix, excluding diagonals c option=0: no reordering performed, only ia and ja created c =1: reordering and sorting performed, a,d also created c I. Misztal, M. Perez-Enciso implicit none integer n,uia(1),uja(1),ia(1),ja(1) + ,i,j,k,tmp(1) c count the number of entries above and below diagonal do i=1,n tmp(i)=0 enddo do i=1,n do k=uia(i),uia(i+1)-1 j=uja(k) if(i.ne.j .and. j.le.n .and. j.ne.0) then tmp(i)=tmp(i)+1 tmp(j)=tmp(j)+1 endif enddo enddo c fill ia ( <--> xadj) ia(1)=1 do i=1,n ia(i+1)=ia(i)+tmp(i) tmp(i)=0 enddo c fill ja ( <--> adjncy ) do i=1,n do k=uia(i),uia(i+1)-1 j=uja(k) if(j.ne.i .and. j.le.n .and. j.ne.0) then ja(ia(i)+tmp(i))=j ja(ia(j)+tmp(j))=i tmp(i)=tmp(i)+1 tmp(j)=tmp(j)+1 endif enddo enddo return end c --------------- subroutine zero (n,x) c --------------- integer n,x(n),i do i=1,n x(i)=0 enddo end c ----------------- subroutine zerout c ----------------- + (b,path,sstart,nseg) integer path(1),sstart(1),nseg real*8 b(1) do i=1,sstart(nseg+1)-1 b(path(i))=0. enddo return end C IM 10/24/93 C*********************************************************************** C Sinverse1 -- inverse OF SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEM OF C LINEAR EQUATIONS Z=inv(M) GIVEN UT-D-U FACTORIZATION OF M; C C Replaces a in input matrix ia-ja-a with the corresponding c inverse elements C C Compared to Sinverse, this version has smaller memory C requirements and does not expand ia-ja-a to full storage C*********************************************************************** SUBROUTINE fsinverse1 * (N, P, IP, D, IJU,JU,IU,U, TMPu,tmpz,tmpi,tmpi1, * jja,ia,ja,a) implicit none INTEGER P(1), IP(1), IJU(1), JU(1), IU(1), TMPi1(1), * tmpi(1),jja(1),ia(1),ja(1),n,j1,k1 real*8 D(1), U(1), a(1), TMPu(1),tmpz(1) integer i,j,k,iiu C C ADDITIONAL PARAMETERS C C TMPu,tmpz - real ONE-DIMENSIONAL WORK ARRAYs; DIMENSION = N C jja - integer on-dimensional work array of order iu(n+1)-1 C tmpi,tmpi1 - integer one-dimensional work array; dimension=N, c tmpi and tmpi1 can occupy the same storage as c tmpu and tmpz C C----------------------------------------------------------------------- iiu(i,j)=ju(iju(i)-iu(i)+j) c c CAUTION: U is really -U C----invert in site, z overwrites u C tmpu keeps unfolded current column of u; C tmpz keeps unfolded current column of z; do i=1,n tmpz(i)=0 tmpu(i)=0 enddo do i=n,1,-1 cdir$ ivdep do j=iu(i), iu(i+1)-1 j1=iiu(i,j) tmpu(j1)=u(j) tmpz(j1)=u(j)*d(j1) enddo c-------off-diagonal elements do j=iu(i),iu(i+1)-1 j1=iiu(i,j) c tmpz(j1)=tmpz(j1)+tmpu(j1)*d(j1) cdir$ ivdep do k=iu(j1),iu(j1+1)-1 k1=iiu(j1,k) tmpz(j1)=tmpz(j1)+tmpu(k1)*u(k) tmpz(k1)=tmpz(k1)+tmpu(j1)*u(k) enddo enddo c------diagonal element last cdir$ ivdep do j=iu(i),iu(i+1)-1 d(i)=d(i)+u(j)*tmpz(iiu(i,j)) enddo c store inverse and zero all nonzeroes in tmp and tmp1 cdir$ ivdep do j=iu(i), iu(i+1)-1 tmpu(iiu(i,j))=0 u(j)=tmpz(iiu(i,j)) tmpz(iiu(i,j))=0 enddo enddo c replace a in ia-ja-a with inverse elements do i=1,n tmpi1(i)=i enddo c first partly permute ia-ja-a call pperm(n,ip,ia,ja,a,tmpi,jja) do i=1,n do j=iu(ip(i)),iu(ip(i)+1)-1 tmpu(p(iiu(ip(i),j)))=u(j) enddo tmpu(i)=d(ip(i)) do j=ia(i),ia(i+1)-1 a(j)=tmpu(ja(j)) enddo enddo c reverse the permutation call pperm(n,tmpi1,ia,ja,a,tmpi,jja) end C IM - 10-19-1993 C********************************************************************* C Partly permute upper-triangular matrix ia-ja-a using permutation vector C p. C The permutation is such that it does not change the indices of elements, C but insures, that after permutation all elements would be in the C upper diagonal. C********************************************************************* subroutine pperm(n,p,ia,ja,a,tmpia,tmprow) implicit none integer n,p(1),ia(1),ja(1),tmpia(1),tmprow(1),i,j,k real*8 a(1),x c extra parameters c tmpia - integer vector of size n c tmprow - integer vector of size ia(n+1)-n c zero temporary row vector for reorder matrix do i=1,n tmpia(i)=0 enddo c set ja to store column entries and tmprow of row entries do i=1,n do j=ia(i),ia(i+1)-1 if (p(i).gt.p(ja(j))) then k=ja(j) ja(j)=i tmprow(j)=k else tmprow(j)=i k=i endif tmpia(k)=tmpia(k)+1 enddo enddo c set ia to permuted ia do i=1,n ia(i+1)=ia(i)+tmpia(i) tmpia(i)=ia(i+1) enddo c assign rows to appropriate addresses do i=ia(n+1)-1,ia(1),-1 j=tmprow(i) tmpia(j)=tmpia(j)-1 tmprow(i)=tmpia(j) enddo c final permutation do i=1,ia(n+1)-1 10 j=tmprow(i) if (i.ne.j) then c swap entries i and j k=ja(i) ja(i)=ja(j) ja(j)=k k=tmprow(i) tmprow(i)=tmprow(j) tmprow(j)=k x=a(i) a(i)=a(j) a(j)=x goto 10 endif enddo end C IM - 1/10/92 - 10/24/93 C*********************************************************************** C Sinverse -- inverse OF SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEM OF C LINEAR EQUATIONS Z=inv(M) GIVEN UT-D-U FACTORIZATION OF M; C C Expands input matrix ia-ja-a to full storage and repalces c cvalues of a with those of the inverse C*********************************************************************** SUBROUTINE fsinverse * (N, P, IP, D, IJU,JU,IU,U, TMP,za,zd,zja, * utia,utja,tmpi,tmpi1,ia,ja,a) INTEGER P(1), IP(1), IJU(1), JU(1), IU(1), TMPi1(1), * zja(1), tmpi(1), * utia(1),utja(1),ia(1),ja(1) real*8 D(1), U(1), za(1), a(1), TMP(1),zd(1) integer i,j,k,l,freeuz,nnzero C C ADDITIONAL PARAMETERS C C TMP - real ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C iu,zja,za - inverse, the lower triangular form, in the i-j-a C form; diagonals in ZD C utia,utja - U in the i-j-a format, indices only; sorted within rows C tmpi,tmpi1 - integer one-dimensional work array; dimension=N C C----------------------------------------------------------------------- C c CAUTION: U is really -U c create indices for U transposed, sorted, za used as temporary call rowtocolu(iu,ju,iju,u,utia,utja,za,tmpi,n) c create I index for Z the same as in U; J index is c added during computations C----invert, tmpi contains number of filled elements in each row of z, C tmp keeps unfolded current column of z; c remember: U is -U do i=1,n tmpi(i)=0 tmp(i)=0 enddo do i=n,1,-1 c------scatter existing elements of column i of z into tmp do j=iu(i), iu(i)+tmpi(i)-1 tmp(zja(j))=za(j) enddo c------diagonal element first freeuz=iu(i)+tmpi(i) tmp(i)=d(i) do j=iu(i),iu(i+1)-1 tmp(i)=tmp(i)+u(j)*tmp(ju(iju(i)-iu(i)+j)) enddo zd(i)=tmp(i) nnzero=1 tmpi1(nnzero)=i c-------off-diagonal elements do j=utia(i+1)-1,utia(i),-1 freeuz=iu(utja(j))+tmpi(utja(j)) l=utja(j) cdir$ ivdep do k=iu(l),iu(l+1)-1 tmp(l)=tmp(l)+u(k)*tmp(ju(iju(l)-iu(l)+k)) enddo za(freeuz)=tmp(l) zja(freeuz)=i tmpi(utja(j))=tmpi(utja(j))+1 nnzero=nnzero+1 tmpi1(nnzero)=l enddo c zero all nonzeroes in tmp do j=1,nnzero tmp(tmpi1(j))=0 enddo enddo c in Z, zero all entries not in input matrix A c ----- since A is upper diagonal, create index structure for A transpose call rowtocol(ia,ja,a,utia,utja,a,tmpi,n) call zero(n,tmpi) call zero(n,tmpi1) do i=1,n c --------scatter row i of A (upper + lower triangle do j=ia(i),ia(i+1)-1 tmpi(ja(j))=1 enddo do j=utia(i),utia(i+1)-1 tmpi(utja(j))=1 enddo do j=iu(ip(i)),iu(ip(i)+1)-1 if (tmpi(p(zja(j))).eq.0) then za(j)=0 endif enddo c--------zero tmpi and count number of entries in upper+lower diagonal a do j=ia(i),ia(i+1)-1 tmpi(ja(j))=0 k=ja(j) if (i.ne.k) tmpi1(i)=tmpi1(i)+1 tmpi1(k)=tmpi1(k)+1 enddo do j=utia(i),utia(i+1)-1 tmpi(utja(j))=0 enddo enddo c create row structure for upper+lower diagonal a ia(1)=1 do i=1,n ia(i+1)=ia(i)+tmpi1(i) enddo c copy z to a, from half to full stored call zero(n,tmpi) do i=1,n k=p(i) a(ia(k)+tmpi(k))=zd(i) ja(ia(k)+tmpi(k))=k tmpi(k)=tmpi(k)+1 do j=iu(i),iu(i+1)-1 if (za(j).ne.0) then k=p(zja(j)) l=p(i) if (k.eq.l) then a(ia(k)+tmpi(k))=za(j) ja(ia(k)+tmpi(k))=k tmpi(k)=tmpi(k)+1 elseif (l.ne.k .and. k.ne.0) then a(ia(l)+tmpi(l))=za(j) a(ia(k)+tmpi(k))=za(j) ja(ia(l)+tmpi(l))=k ja(ia(k)+tmpi(k))=l tmpi(l)=tmpi(l)+1 tmpi(k)=tmpi(k)+1 endif endif enddo enddo return end C----- SUBROUTINE GENQMD C**************************************************************** 1. C**************************************************************** 2. C********** GENQMD ..... QUOT MIN DEGREE ORDERING ********* 3. C**************************************************************** 4. C**************************************************************** 5. C 6. C PURPOSE - THIS ROUTINE IMPLEMENTS THE MINIMUM DEGREE 7. C ALGORITHM. IT MAKES USE OF THE IMPLICIT REPRESENT- 8. C ATION OF THE ELIMINATION GRAPHS BY QUOTIENT GRAPHS, 9. C AND THE NOTION OF INDISTINGUISHABLE NODES. 10. C CAUTION - THE ADJACENCY VECTOR ADJNCY WILL BE 11. C DESTROYED. 12. C 13. C INPUT PARAMETERS - 14. C NEQNS - NUMBER OF EQUATIONS. 15. C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. 16. C 17. C OUTPUT PARAMETERS - 18. C PERM - THE MINIMUM DEGREE ORDERING. 19. C INVP - THE INVERSE OF PERM. 20. C 21. C WORKING PARAMETERS - 22. C DEG - THE DEGREE VECTOR. DEG(I) IS NEGATIVE MEANS 23. C NODE I HAS BEEN NUMBERED. 24. C MARKER - A MARKER VECTOR, WHERE MARKER(I) IS 25. C NEGATIVE MEANS NODE I HAS BEEN MERGED WITH 26. C ANOTHER NODE AND THUS CAN BE IGNORED. 27. C RCHSET - VECTOR USED FOR THE REACHABLE SET. 28. C NBRHD - VECTOR USED FOR THE NEIGHBORHOOD SET. 29. C QSIZE - VECTOR USED TO STORE THE SIZE OF 30. C INDISTINGUISHABLE SUPERNODES. 31. C QLINK - VECTOR TO STORE INDISTINGUISHABLE NODES, 32. C I, QLINK(I), QLINK(QLINK(I)) ... ARE THE 33. C MEMBERS OF THE SUPERNODE REPRESENTED BY I. 34. C 35. C PROGRAM SUBROUTINES - 36. C QMDRCH, QMDQT, QMDUPD. 37. C 38. C**************************************************************** 39. C 40. C 41. SUBROUTINE GENQMD ( NEQNS, XADJ, ADJNCY, PERM, INVP, DEG, 42. 1 MARKER, RCHSET, NBRHD, QSIZE, QLINK, 43. 1 NOFSUB ) 44. C 45. C**************************************************************** 46. C 47. INTEGER ADJNCY(1), PERM(1), INVP(1), DEG(1), MARKER(1), 48. 1 RCHSET(1), NBRHD(1), QSIZE(1), QLINK(1) 49. INTEGER XADJ(1), INODE, IP, IRCH, J, MINDEG, NDEG, 50. 1 NEQNS, NHDSZE, NODE, NOFSUB, NP, NUM, NUMP1, 51. 1 NXNODE, RCHSZE, SEARCH, THRESH 52. C 53. C**************************************************************** 54. C 55. C ----------------------------------------------------- 56. C INITIALIZE DEGREE VECTOR AND OTHER WORKING VARIABLES. 57. C ----------------------------------------------------- 58. MINDEG = NEQNS 59. NOFSUB = 0 60. DO 100 NODE = 1, NEQNS 61. PERM(NODE) = NODE 62. INVP(NODE) = NODE 63. MARKER(NODE) = 0 64. QSIZE(NODE) = 1 65. QLINK(NODE) = 0 66. NDEG = XADJ(NODE+1) - XADJ(NODE) 67. DEG(NODE) = NDEG 68. IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 69. 100 CONTINUE 70. NUM = 0 71. C ----------------------------------------------------- 72. C PERFORM THRESHOLD SEARCH TO GET A NODE OF MIN DEGREE. 73. C VARIABLE SEARCH POINTS TO WHERE SEARCH SHOULD START. 74. C ----------------------------------------------------- 75. 200 SEARCH = 1 76. THRESH = MINDEG 77. MINDEG = NEQNS 78. 300 NUMP1 = NUM + 1 79. IF ( NUMP1 .GT. SEARCH ) SEARCH = NUMP1 80. DO 400 J = SEARCH, NEQNS 81. NODE = PERM(J) 82. IF ( MARKER(NODE) .LT. 0 ) GOTO 400 83. NDEG = DEG(NODE) 84. IF ( NDEG .LE. THRESH ) GO TO 500 85. IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 86. 400 CONTINUE 87. GO TO 200 88. C --------------------------------------------------- 89. C NODE HAS MINIMUM DEGREE. FIND ITS REACHABLE SETS BY 90. C CALLING QMDRCH. 91. C --------------------------------------------------- 92. 500 SEARCH = J 93. NOFSUB = NOFSUB + DEG(NODE) 94. MARKER(NODE) = 1 95. CALL QMDRCH (NODE, XADJ, ADJNCY, DEG, MARKER, 96. 1 RCHSZE, RCHSET, NHDSZE, NBRHD ) 97. C ------------------------------------------------ 98. C ELIMINATE ALL NODES INDISTINGUISHABLE FROM NODE. 99. C THEY ARE GIVEN BY NODE, QLINK(NODE), .... 100. C ------------------------------------------------ 101. NXNODE = NODE 102. 600 NUM = NUM + 1 103. NP = INVP(NXNODE) 104. IP = PERM(NUM) 105. PERM(NP) = IP 106. INVP(IP) = NP 107. PERM(NUM) = NXNODE 108. INVP(NXNODE) = NUM 109. DEG(NXNODE) = - 1 110. NXNODE = QLINK(NXNODE) 111. IF (NXNODE .GT. 0) GOTO 600 112. C 113. IF ( RCHSZE .LE. 0 ) GO TO 800 114. C ------------------------------------------------ 115. C UPDATE THE DEGREES OF THE NODES IN THE REACHABLE 116. C SET AND IDENTIFY INDISTINGUISHABLE NODES. 117. C ------------------------------------------------ 118. CALL QMDUPD ( XADJ, ADJNCY, RCHSZE, RCHSET, DEG, 119. 1 QSIZE, QLINK, MARKER, RCHSET(RCHSZE+1), 120. 1 NBRHD(NHDSZE+1) ) 121. C ------------------------------------------- 122. C RESET MARKER VALUE OF NODES IN REACH SET. 123. C UPDATE THRESHOLD VALUE FOR CYCLIC SEARCH. 124. C ALSO CALL QMDQT TO FORM NEW QUOTIENT GRAPH. 125. C ------------------------------------------- 126. MARKER(NODE) = 0 127. DO 700 IRCH = 1, RCHSZE 128. INODE = RCHSET(IRCH) 129. IF ( MARKER(INODE) .LT. 0 ) GOTO 700 130. MARKER(INODE) = 0 131. NDEG = DEG(INODE) 132. IF ( NDEG .LT. MINDEG ) MINDEG = NDEG 133. IF ( NDEG .GT. THRESH ) GOTO 700 134. MINDEG = THRESH 135. THRESH = NDEG 136. SEARCH = INVP(INODE) 137. 700 CONTINUE 138. IF ( NHDSZE .GT. 0 ) CALL QMDQT ( NODE, XADJ, 139. 1 ADJNCY, MARKER, RCHSZE, RCHSET, NBRHD ) 140. 800 IF ( NUM .LT. NEQNS ) GO TO 300 141. RETURN 142. END 143. C----- SUBROUTINE QMDUPD C**************************************************************** 1. C**************************************************************** 2. C********** QMDUPD ..... QUOT MIN DEG UPDATE *********** 3. C**************************************************************** 4. C**************************************************************** 5. C 6. C PURPOSE - THIS ROUTINE PERFORMS DEGREE UPDATE FOR A SET 7. C OF NODES IN THE MINIMUM DEGREE ALGORITHM. 8. C 9. C INPUT PARAMETERS - 10. C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. 11. C (NLIST, LIST) - THE LIST OF NODES WHOSE DEGREE HAS TO 12. C BE UPDATED. 13. C 14. C UPDATED PARAMETERS - 15. C DEG - THE DEGREE VECTOR. 16. C QSIZE - SIZE OF INDISTINGUISHABLE SUPERNODES. 17. C QLINK - LINKED LIST FOR INDISTINGUISHABLE NODES. 18. C MARKER - USED TO MARK THOSE NODES IN REACH/NBRHD SETS. 19. C 20. C WORKING PARAMETERS - 21. C RCHSET - THE REACHABLE SET. 22. C NBRHD - THE NEIGHBORHOOD SET. 23. C 24. C PROGRAM SUBROUTINES - 25. C QMDMRG. 26. C 27. C**************************************************************** 28. C 29. SUBROUTINE QMDUPD ( XADJ, ADJNCY, NLIST, LIST, DEG, 30. 1 QSIZE, QLINK, MARKER, RCHSET, NBRHD ) 31. C 32. C**************************************************************** 33. C 34. INTEGER ADJNCY(1), LIST(1), DEG(1), MARKER(1), 35. 1 RCHSET(1), NBRHD(1), QSIZE(1), QLINK(1) 36. INTEGER XADJ(1), DEG0, DEG1, IL, INHD, INODE, IRCH, 37. 1 J, JSTRT, JSTOP, MARK, NABOR, NHDSZE, NLIST, 38. 1 NODE, RCHSZE, ROOT 39. C 40. C**************************************************************** 41. C 42. C ------------------------------------------------ 43. C FIND ALL ELIMINATED SUPERNODES THAT ARE ADJACENT 44. C TO SOME NODES IN THE GIVEN LIST. PUT THEM INTO 45. C (NHDSZE, NBRHD). DEG0 CONTAINS THE NUMBER OF 46. C NODES IN THE LIST. 47. C ------------------------------------------------ 48. IF ( NLIST .LE. 0 ) RETURN 49. DEG0 = 0 50. NHDSZE = 0 51. DO 200 IL = 1, NLIST 52. NODE = LIST(IL) 53. DEG0 = DEG0 + QSIZE(NODE) 54. JSTRT = XADJ(NODE) 55. JSTOP = XADJ(NODE+1) - 1 56. DO 100 J = JSTRT, JSTOP 57. NABOR = ADJNCY(J) 58. IF ( MARKER(NABOR) .NE. 0 .OR. 59. 1 DEG(NABOR) .GE. 0 ) GO TO 100 60. MARKER(NABOR) = - 1 61. NHDSZE = NHDSZE + 1 62. NBRHD(NHDSZE) = NABOR 63. 100 CONTINUE 64. 200 CONTINUE 65. C -------------------------------------------- 66. C MERGE INDISTINGUISHABLE NODES IN THE LIST BY 67. C CALLING THE SUBROUTINE QMDMRG. 68. C -------------------------------------------- 69. IF ( NHDSZE .GT. 0 ) 70. 1 CALL QMDMRG ( XADJ, ADJNCY, DEG, QSIZE, QLINK, 71. 1 MARKER, DEG0, NHDSZE, NBRHD, RCHSET, 72. 1 NBRHD(NHDSZE+1) ) 73. C ---------------------------------------------------- 74. C FIND THE NEW DEGREES OF THE NODES THAT HAVE NOT BEEN 75. C MERGED. 76. C ---------------------------------------------------- 77. DO 600 IL = 1, NLIST 78. NODE = LIST(IL) 79. MARK = MARKER(NODE) 80. IF ( MARK .GT. 1 .OR. MARK .LT. 0 ) GO TO 600 81. MARKER(NODE) = 2 82. CALL QMDRCH ( NODE, XADJ, ADJNCY, DEG, MARKER, 83. 1 RCHSZE, RCHSET, NHDSZE, NBRHD ) 84. DEG1 = DEG0 85. IF ( RCHSZE .LE. 0 ) GO TO 400 86. DO 300 IRCH = 1, RCHSZE 87. INODE = RCHSET(IRCH) 88. DEG1 = DEG1 + QSIZE(INODE) 89. MARKER(INODE) = 0 90. 300 CONTINUE 91. 400 DEG(NODE) = DEG1 - 1 92. IF ( NHDSZE .LE. 0 ) GO TO 600 93. DO 500 INHD = 1, NHDSZE 94. INODE = NBRHD(INHD) 95. MARKER(INODE) = 0 96. 500 CONTINUE 97. 600 CONTINUE 98. RETURN 99. END 100. C----- SUBROUTINE QMDQT C************************************************************* 1. C************************************************************* 2. C******* QMDQT ..... QUOT MIN DEG QUOT TRANSFORM ******* 3. C************************************************************* 4. C************************************************************* 5. C 6. C PURPOSE - THIS SUBROUTINE PERFORMS THE QUOTIENT GRAPH 7. C TRANSFORMATION AFTER A NODE HAS BEEN ELIMINATED. 8. C 9. C INPUT PARAMETERS - 10. C ROOT - THE NODE JUST ELIMINATED. IT BECOMES THE 11. C REPRESENTATIVE OF THE NEW SUPERNODE. 12. C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. 13. C (RCHSZE, RCHSET) - THE REACHABLE SET OF ROOT IN THE 14. C OLD QUOTIENT GRAPH. 15. C NBRHD - THE NEIGHBORHOOD SET WHICH WILL BE MERGED 16. C WITH ROOT TO FORM THE NEW SUPERNODE. 17. C MARKER - THE MARKER VECTOR. 18. C 19. C UPDATED PARAMETER - 20. C ADJNCY - BECOMES THE ADJNCY OF THE QUOTIENT GRAPH. 21. C 22. C************************************************************* 23. C 24. SUBROUTINE QMDQT ( ROOT, XADJ, ADJNCY, MARKER, 25. 1 RCHSZE, RCHSET, NBRHD ) 26. C 27. C************************************************************* 28. C 29. INTEGER ADJNCY(1), MARKER(1), RCHSET(1), NBRHD(1) 30. INTEGER XADJ(1), INHD, IRCH, J, JSTRT, JSTOP, LINK, 31. 1 NABOR, NODE, RCHSZE, ROOT 32. C 33. C************************************************************* 34. C 35. IRCH = 0 36. INHD = 0 37. NODE = ROOT 38. 100 JSTRT = XADJ(NODE) 39. JSTOP = XADJ(NODE+1) - 2 40. IF ( JSTOP .LT. JSTRT ) GO TO 300 41. C ------------------------------------------------ 42. C PLACE REACH NODES INTO THE ADJACENT LIST OF NODE 43. C ------------------------------------------------ 44. DO 200 J = JSTRT, JSTOP 45. IRCH = IRCH + 1 46. ADJNCY(J) = RCHSET(IRCH) 47. IF ( IRCH .GE. RCHSZE ) GOTO 400 48. 200 CONTINUE 49. C ---------------------------------------------- 50. C LINK TO OTHER SPACE PROVIDED BY THE NBRHD SET. 51. C ---------------------------------------------- 52. 300 LINK = ADJNCY(JSTOP+1) 53. NODE = - LINK 54. IF ( LINK .LT. 0 ) GOTO 100 55. INHD = INHD + 1 56. NODE = NBRHD(INHD) 57. ADJNCY(JSTOP+1) = - NODE 58. GO TO 100 59. C ------------------------------------------------------- 60. C ALL REACHABLE NODES HAVE BEEN SAVED. END THE ADJ LIST. 61. C ADD ROOT TO THE NBR LIST OF EACH NODE IN THE REACH SET. 62. C ------------------------------------------------------- 63. 400 ADJNCY(J+1) = 0 64. DO 600 IRCH = 1, RCHSZE 65. NODE = RCHSET(IRCH) 66. IF ( MARKER(NODE) .LT. 0 ) GOTO 600 67. JSTRT = XADJ(NODE) 68. JSTOP = XADJ(NODE+1) - 1 69. DO 500 J = JSTRT, JSTOP 70. NABOR = ADJNCY(J) 71. IF ( MARKER(NABOR) .GE. 0 ) GO TO 500 72. ADJNCY(J) = ROOT 73. GOTO 600 74. 500 CONTINUE 75. 600 CONTINUE 76. RETURN 77. END 78. C----- SUBROUTINE QMDMRG C**************************************************************** 1. C**************************************************************** 2. C********** QMDMRG ..... QUOT MIN DEG MERGE *********** 3. C**************************************************************** 4. C**************************************************************** 5. C 6. C PURPOSE - THIS ROUTINE MERGES INDISTINGUISHABLE NODES IN 7. C THE MINIMUM DEGREE ORDERING ALGORITHM. 8. C IT ALSO COMPUTES THE NEW DEGREES OF THESE 9. C NEW SUPERNODES. 10. C 11. C INPUT PARAMETERS - 12. C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. 13. C DEG0 - THE NUMBER OF NODES IN THE GIVEN SET. 14. C (NHDSZE, NBRHD) - THE SET OF ELIMINATED SUPERNODES 15. C ADJACENT TO SOME NODES IN THE SET. 16. C 17. C UPDATED PARAMETERS - 18. C DEG - THE DEGREE VECTOR. 19. C QSIZE - SIZE OF INDISTINGUISHABLE NODES. 20. C QLINK - LINKED LIST FOR INDISTINGUISHABLE NODES. 21. C MARKER - THE GIVEN SET IS GIVEN BY THOSE NODES WITH 22. C MARKER VALUE SET TO 1. THOSE NODES WITH DEGREE 23. C UPDATED WILL HAVE MARKER VALUE SET TO 2. 24. C 25. C WORKING PARAMETERS - 26. C RCHSET - THE REACHABLE SET. 27. C OVRLP - TEMP VECTOR TO STORE THE INTERSECTION OF TWO 28. C REACHABLE SETS. 29. C 30. C**************************************************************** 31. C 32. SUBROUTINE QMDMRG ( XADJ, ADJNCY, DEG, QSIZE, QLINK, 33. 1 MARKER, DEG0, NHDSZE, NBRHD, RCHSET, 34. 1 OVRLP ) 35. C 36. C**************************************************************** 37. C 38. INTEGER ADJNCY(1), DEG(1), QSIZE(1), QLINK(1), 39. 1 MARKER(1), RCHSET(1), NBRHD(1), OVRLP(1) 40. INTEGER XADJ(1), DEG0, DEG1, HEAD, INHD, IOV, IRCH, 41. 1 J, JSTRT, JSTOP, LINK, LNODE, MARK, MRGSZE, 42. 1 NABOR, NHDSZE, NODE, NOVRLP, RCHSZE, ROOT 43. C 44. C**************************************************************** 45. C 46. C ------------------ 47. C INITIALIZATION ... 48. C ------------------ 49. IF ( NHDSZE .LE. 0 ) RETURN 50. DO 100 INHD = 1, NHDSZE 51. ROOT = NBRHD(INHD) 52. MARKER(ROOT) = 0 53. 100 CONTINUE 54. C ------------------------------------------------- 55. C LOOP THROUGH EACH ELIMINATED SUPERNODE IN THE SET 56. C (NHDSZE, NBRHD). 57. C ------------------------------------------------- 58. DO 1400 INHD = 1, NHDSZE 59. ROOT = NBRHD(INHD) 60. MARKER(ROOT) = - 1 61. RCHSZE = 0 62. NOVRLP = 0 63. DEG1 = 0 64. 200 JSTRT = XADJ(ROOT) 65. JSTOP = XADJ(ROOT+1) - 1 66. C ---------------------------------------------- 67. C DETERMINE THE REACHABLE SET AND ITS INTERSECT- 68. C ION WITH THE INPUT REACHABLE SET. 69. C ---------------------------------------------- 70. DO 600 J = JSTRT, JSTOP 71. NABOR = ADJNCY(J) 72. ROOT = - NABOR 73. IF (NABOR) 200, 700, 300 74. C 75. 300 MARK = MARKER(NABOR) 76. IF ( MARK ) 600, 400, 500 77. 400 RCHSZE = RCHSZE + 1 78. RCHSET(RCHSZE) = NABOR 79. DEG1 = DEG1 + QSIZE(NABOR) 80. MARKER(NABOR) = 1 81. GOTO 600 82. 500 IF ( MARK .GT. 1 ) GOTO 600 83. NOVRLP = NOVRLP + 1 84. OVRLP(NOVRLP) = NABOR 85. MARKER(NABOR) = 2 86. 600 CONTINUE 87. C -------------------------------------------- 88. C FROM THE OVERLAPPED SET, DETERMINE THE NODES 89. C THAT CAN BE MERGED TOGETHER. 90. C -------------------------------------------- 91. 700 HEAD = 0 92. MRGSZE = 0 93. DO 1100 IOV = 1, NOVRLP 94. NODE = OVRLP(IOV) 95. JSTRT = XADJ(NODE) 96. JSTOP = XADJ(NODE+1) - 1 97. DO 800 J = JSTRT, JSTOP 98. NABOR = ADJNCY(J) 99. IF ( MARKER(NABOR) .NE. 0 ) GOTO 800 100. MARKER(NODE) = 1 101. GOTO 1100 102. 800 CONTINUE 103. C ----------------------------------------- 104. C NODE BELONGS TO THE NEW MERGED SUPERNODE. 105. C UPDATE THE VECTORS QLINK AND QSIZE. 106. C ----------------------------------------- 107. MRGSZE = MRGSZE + QSIZE(NODE) 108. MARKER(NODE) = - 1 109. LNODE = NODE 110. 900 LINK = QLINK(LNODE) 111. IF ( LINK .LE. 0 ) GOTO 1000 112. LNODE = LINK 113. GOTO 900 114. 1000 QLINK(LNODE) = HEAD 115. HEAD = NODE 116. 1100 CONTINUE 117. IF ( HEAD .LE. 0 ) GOTO 1200 118. QSIZE(HEAD) = MRGSZE 119. DEG(HEAD) = DEG0 + DEG1 - 1 120. MARKER(HEAD) = 2 121. C -------------------- 122. C RESET MARKER VALUES. 123. C -------------------- 124. 1200 ROOT = NBRHD(INHD) 125. MARKER(ROOT) = 0 126. IF ( RCHSZE .LE. 0 ) GOTO 1400 127. DO 1300 IRCH = 1, RCHSZE 128. NODE = RCHSET(IRCH) 129. MARKER(NODE) = 0 130. 1300 CONTINUE 131. 1400 CONTINUE 132. RETURN 133. END 134. C----- SUBROUTINE QMDRCH C*************************************************************** 1. C*************************************************************** 2. C********* QMDRCH ..... QUOT MIN DEG REACH SET ********** 3. C*************************************************************** 4. C*************************************************************** 5. C 6. C PURPOSE - THIS SUBROUTINE DETERMINES THE REACHABLE SET OF 7. C A NODE THROUGH A GIVEN SUBSET. THE ADJACENCY STRUCTURE 8. C IS ASSUMED TO BE STORED IN A QUOTIENT GRAPH FORMAT. 9. C 10. C INPUT PARAMETERS - 11. C ROOT - THE GIVEN NODE NOT IN THE SUBSET. 12. C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE PAIR. 13. C DEG - THE DEGREE VECTOR. DEG(I) LT 0 MEANS THE NODE 14. C BELONGS TO THE GIVEN SUBSET. 15. C 16. C OUTPUT PARAMETERS - 17. C (RCHSZE, RCHSET) - THE REACHABLE SET. 18. C (NHDSZE, NBRHD) - THE NEIGHBORHOOD SET. 19. C 20. C UPDATED PARAMETERS - 21. C MARKER - THE MARKER VECTOR FOR REACH AND NBRHD SETS. 22. C GT 0 MEANS THE NODE IS IN REACH SET. 23. C LT 0 MEANS THE NODE HAS BEEN MERGED WITH 24. C OTHERS IN THE QUOTIENT OR IT IS IN NBRHD SET. 25. C 26. C*************************************************************** 27. C 28. SUBROUTINE QMDRCH ( ROOT, XADJ, ADJNCY, DEG, MARKER, 29. 1 RCHSZE, RCHSET, NHDSZE, NBRHD ) 30. C 31. C*************************************************************** 32. C 33. INTEGER ADJNCY(1), DEG(1), MARKER(1), 34. 1 RCHSET(1), NBRHD(1) 35. INTEGER XADJ(1), I, ISTRT, ISTOP, J, JSTRT, JSTOP, 36. 1 NABOR, NHDSZE, NODE, RCHSZE, ROOT 37. C 38. C*************************************************************** 39. C 40. C ----------------------------------------- 41. C LOOP THROUGH THE NEIGHBORS OF ROOT IN THE 42. C QUOTIENT GRAPH. 43. C ----------------------------------------- 44. NHDSZE = 0 45. RCHSZE = 0 46. ISTRT = XADJ(ROOT) 47. ISTOP = XADJ(ROOT+1) - 1 48. IF ( ISTOP .LT. ISTRT ) RETURN 49. DO 600 I = ISTRT, ISTOP 50. NABOR = ADJNCY(I) 51. IF ( NABOR .EQ. 0 ) RETURN 52. IF ( MARKER(NABOR) .NE. 0 ) GO TO 600 53. IF ( DEG(NABOR) .LT. 0 ) GO TO 200 54. C ------------------------------------- 55. C INCLUDE NABOR INTO THE REACHABLE SET. 56. C ------------------------------------- 57. RCHSZE = RCHSZE + 1 58. RCHSET(RCHSZE) = NABOR 59. MARKER(NABOR) = 1 60. GO TO 600 61. C ------------------------------------- 62. C NABOR HAS BEEN ELIMINATED. FIND NODES 63. C REACHABLE FROM IT. 64. C ------------------------------------- 65. 200 MARKER(NABOR) = -1 66. NHDSZE = NHDSZE + 1 67. NBRHD(NHDSZE) = NABOR 68. 300 JSTRT = XADJ(NABOR) 69. JSTOP = XADJ(NABOR+1) - 1 70. DO 500 J = JSTRT, JSTOP 71. NODE = ADJNCY(J) 72. NABOR = - NODE 73. IF (NODE) 300, 600, 400 74. 400 IF ( MARKER(NODE) .NE. 0 ) GO TO 500 75. RCHSZE = RCHSZE + 1 76. RCHSET(RCHSZE) = NODE 77. MARKER(NODE) = 1 78. 500 CONTINUE 79. 600 CONTINUE 80. RETURN 81. END 82. C----- SUBROUTINE SMBFCT C**************************************************************** 1. C**************************************************************** 2. C********* SMBFCT ..... SYMBOLIC FACTORIZATION ******** 3. C**************************************************************** 4. C**************************************************************** 5. C 6. C PURPOSE - THIS ROUTINE PERFORMS SYMBOLIC FACTORIZATION 7. C ON A PERMUTED LINEAR SYSTEM AND IT ALSO SETS UP THE 8. C COMPRESSED DATA STRUCTURE FOR THE SYSTEM. 9. C 10. C INPUT PARAMETERS - 11. C NEQNS - NUMBER OF EQUATIONS. 12. C (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. 13. C (PERM, INVP) - THE PERMUTATION VECTOR AND ITS INVERSE. 14. C 15. C UPDATED PARAMETERS - 16. C MAXSUB - SIZE OF THE SUBSCRIPT ARRAY NZSUB. ON RETURN, 17. C IT CONTAINS THE NUMBER OF SUBSCRIPTS USED 18. C 19. C OUTPUT PARAMETERS - 20. C XLNZ - INDEX INTO THE NONZERO STORAGE VECTOR LNZ. 21. C (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT VECTORS. 22. C MAXLNZ - THE NUMBER OF NONZEROS FOUND. 23. C FLAG - ERROR FLAG. POSITIVE VALUE INDICATES THAT. 24. C NZSUB ARRAY IS TOO SMALL. 25. C 26. C WORKING PARAMETERS - 27. C MRGLNK - A VECTOR OF SIZE NEQNS. AT THE KTH STEP, 28. C MRGLNK(K), MRGLNK(MRGLNK(K)) , ......... 29. C IS A LIST CONTAINING ALL THOSE COLUMNS L(*,J) 30. C WITH J LESS THAN K, SUCH THAT ITS FIRST OFF- 31. C DIAGONAL NONZERO IS L(K,J). THUS, THE 32. C NONZERO STRUCTURE OF COLUMN L(*,K) CAN BE FOUND 33. C BY MERGING THAT OF SUCH COLUMNS L(*,J) WITH 34. C THE STRUCTURE OF A(*,K). 35. C RCHLNK - A VECTOR OF SIZE NEQNS. IT IS USED TO ACCUMULATE 36. C THE STRUCTURE OF EACH COLUMN L(*,K). AT THE 37. C END OF THE KTH STEP, 38. C RCHLNK(K), RCHLNK(RCHLNK(K)), ........ 39. C IS THE LIST OF POSITIONS OF NONZEROS IN COLUMN K 40. C OF THE FACTOR L. 41. C MARKER - AN INTEGER VECTOR OF LENGTH NEQNS. IT IS USED 42. C TO TEST IF MASS SYMBOLIC ELIMINATION CAN BE 43. C PERFORMED. THAT IS, IT IS USED TO CHECK WHETHER 44. C THE STRUCTURE OF THE CURRENT COLUMN K BEING 45. C PROCESSED IS COMPLETELY DETERMINED BY THE SINGLE 46. C COLUMN MRGLNK(K). 47. C 48. C**************************************************************** 49. C 50. SUBROUTINE SMBFCT ( NEQNS, XADJ, ADJNCY, PERM, INVP, 51. 1 XLNZ, MAXLNZ, XNZSUB, NZSUB, MAXSUB, 52. 1 RCHLNK, MRGLNK, MARKER, FLAG ) 53. C 54. C**************************************************************** 55. C 56. INTEGER ADJNCY(1), INVP(1), MRGLNK(1), NZSUB(1), 57. 1 PERM(1), RCHLNK(1), MARKER(1) 58. INTEGER XADJ(1), XLNZ(1), XNZSUB(1), 59. 1 FLAG, I, INZ, J, JSTOP, JSTRT, K, KNZ, 60. 1 KXSUB, MRGK, LMAX, M, MAXLNZ, MAXSUB, 61. 1 NABOR, NEQNS, NODE, NP1, NZBEG, NZEND, 62. 1 RCHM, MRKFLG 63. C 64. C**************************************************************** 65. C 66. C ------------------ 67. C INITIALIZATION ... 68. C ------------------ 69. NZBEG = 1 70. NZEND = 0 71. XLNZ(1) = 1 72. DO 100 K = 1, NEQNS 73. MRGLNK(K) = 0 74. MARKER(K) = 0 75. 100 CONTINUE 76. C -------------------------------------------------- 77. C FOR EACH COLUMN ......... . KNZ COUNTS THE NUMBER 78. C OF NONZEROS IN COLUMN K ACCUMULATED IN RCHLNK. 79. C -------------------------------------------------- 80. NP1 = NEQNS + 1 81. DO 1500 K = 1, NEQNS 82. KNZ = 0 83. MRGK = MRGLNK(K) 84. MRKFLG = 0 85. MARKER(K) = K 86. IF (MRGK .NE. 0 ) MARKER(K) = MARKER(MRGK) 87. XNZSUB(K) = NZEND 88. NODE = PERM(K) 89. JSTRT = XADJ(NODE) 90. JSTOP = XADJ(NODE+1) - 1 91. IF (JSTRT.GT.JSTOP) GO TO 1500 92. C ------------------------------------------- 93. C USE RCHLNK TO LINK THROUGH THE STRUCTURE OF 94. C A(*,K) BELOW DIAGONAL 95. C ------------------------------------------- 96. RCHLNK(K) = NP1 97. DO 300 J = JSTRT, JSTOP 98. NABOR = ADJNCY(J) 99. NABOR = INVP(NABOR) 100. IF ( NABOR .LE. K ) GO TO 300 101. RCHM = K 102. 200 M = RCHM 103. RCHM = RCHLNK(M) 104. IF ( RCHM .LE. NABOR ) GO TO 200 105. KNZ = KNZ+1 106. RCHLNK(M) = NABOR 107. RCHLNK(NABOR) = RCHM 108. IF ( MARKER(NABOR) .NE. MARKER(K) ) MRKFLG = 1 109. 300 CONTINUE 110. C -------------------------------------- 111. C TEST FOR MASS SYMBOLIC ELIMINATION ... 112. C -------------------------------------- 113. LMAX = 0 114. IF ( MRKFLG .NE. 0 .OR. MRGK .EQ. 0 ) GO TO 350 115. IF ( MRGLNK(MRGK) .NE. 0 ) GO TO 350 116. XNZSUB(K) = XNZSUB(MRGK) + 1 117. KNZ = XLNZ(MRGK+1) - (XLNZ(MRGK) + 1) 118. GO TO 1400 119. C ----------------------------------------------- 120. C LINK THROUGH EACH COLUMN I THAT AFFECTS L(*,K). 121. C ----------------------------------------------- 122. 350 I = K 123. 400 I = MRGLNK(I) 124. IF (I.EQ.0) GO TO 800 125. INZ = XLNZ(I+1) - (XLNZ(I)+1) 126. JSTRT = XNZSUB(I) + 1 127. JSTOP = XNZSUB(I) + INZ 128. IF (INZ.LE.LMAX) GO TO 500 129. LMAX = INZ 130. XNZSUB(K) = JSTRT 131. C ----------------------------------------------- 132. C MERGE STRUCTURE OF L(*,I) IN NZSUB INTO RCHLNK. 133. C ----------------------------------------------- 134. 500 RCHM = K 135. DO 700 J = JSTRT, JSTOP 136. NABOR = NZSUB(J) 137. 600 M = RCHM 138. RCHM = RCHLNK(M) 139. IF (RCHM.LT.NABOR) GO TO 600 140. IF (RCHM.EQ.NABOR) GO TO 700 141. KNZ = KNZ+1 142. RCHLNK(M) = NABOR 143. RCHLNK(NABOR) = RCHM 144. RCHM = NABOR 145. 700 CONTINUE 146. GO TO 400 147. C ------------------------------------------------------ 148. C CHECK IF SUBSCRIPTS DUPLICATE THOSE OF ANOTHER COLUMN. 149. C ------------------------------------------------------ 150. 800 IF (KNZ.EQ.LMAX) GO TO 1400 151. C ----------------------------------------------- 152. C OR IF TAIL OF K-1ST COLUMN MATCHES HEAD OF KTH. 153. C ----------------------------------------------- 154. IF (NZBEG.GT.NZEND) GO TO 1200 155. I = RCHLNK(K) 156. DO 900 JSTRT=NZBEG,NZEND 157. IF (NZSUB(JSTRT)-I) 900, 1000, 1200 158. 900 CONTINUE 159. GO TO 1200 160. 1000 XNZSUB(K) = JSTRT 161. DO 1100 J=JSTRT,NZEND 162. IF (NZSUB(J).NE.I) GO TO 1200 163. I = RCHLNK(I) 164. IF (I.GT.NEQNS) GO TO 1400 165. 1100 CONTINUE 166. NZEND = JSTRT - 1 167. C ---------------------------------------- 168. C COPY THE STRUCTURE OF L(*,K) FROM RCHLNK 169. C TO THE DATA STRUCTURE (XNZSUB, NZSUB). 170. C ---------------------------------------- 171. 1200 NZBEG = NZEND + 1 172. NZEND = NZEND + KNZ 173. IF (NZEND.GT.MAXSUB) GO TO 1600 174. I = K 175. DO 1300 J=NZBEG,NZEND 176. I = RCHLNK(I) 177. NZSUB(J) = I 178. MARKER(I) = K 179. 1300 CONTINUE 180. XNZSUB(K) = NZBEG 181. MARKER(K) = K 182. C -------------------------------------------------------- 183. C UPDATE THE VECTOR MRGLNK. NOTE COLUMN L(*,K) JUST FOUND 184. C IS REQUIRED TO DETERMINE COLUMN L(*,J), WHERE 185. C L(J,K) IS THE FIRST NONZERO IN L(*,K) BELOW DIAGONAL. 186. C -------------------------------------------------------- 187. 1400 IF (KNZ.LE.1) GO TO 1500 188. KXSUB = XNZSUB(K) 189. I = NZSUB(KXSUB) 190. MRGLNK(K) = MRGLNK(I) 191. MRGLNK(I) = K 192. 1500 XLNZ(K+1) = XLNZ(K) + KNZ 193. MAXLNZ = XLNZ(NEQNS) - 1 194. MAXSUB = XNZSUB(NEQNS) 195. XNZSUB(NEQNS+1) = XNZSUB(NEQNS) 196. FLAG = 0 197. RETURN 198. C ---------------------------------------------------- 199. C ERROR - INSUFFICIENT STORAGE FOR NONZERO SUBSCRIPTS. 200. C ---------------------------------------------------- 201. 1600 FLAG = 1 202. RETURN 203. END 204. C----- SUBROUTINE GSFCT C*************************************************************** 1. C*************************************************************** 2. C****** GSFCT ..... GENERAL SPARSE SYMMETRIC FACT ****** 3. C*************************************************************** 4. C*************************************************************** 5. C 6. C PURPOSE - THIS SUBROUTINE PERFORMS THE SYMMETRIC 7. C FACTORIZATION FOR A GENERAL SPARSE SYSTEM, STORED IN 8. C THE COMPRESSED SUBSCRIPT DATA FORMAT. 9. C 10. C INPUT PARAMETERS - 11. C NEQNS - NUMBER OF EQUATIONS. 12. C XLNZ - INDEX VECTOR FOR LNZ. XLNZ(I) POINTS TO THE 13. C START OF NONZEROS IN COLUMN I OF FACTOR L. 14. C (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT DATA 15. C STRUCTURE FOR FACTOR L. 16. C 17. C UPDATED PARAMETERS - 18. C LNZ - ON INPUT, CONTAINS NONZEROS OF A, AND ON 19. C RETURN, THE NONZEROS OF L. 20. C DIAG - THE DIAGONAL OF L OVERWRITES THAT OF A. 21. C IFLAG - THE ERROR FLAG. IT IS SET TO 1 IF A ZERO OR 22. C NEGATIVE SQUARE ROOT OCCURS DURING THE 23. C FACTORIZATION. 24. C OPS - A DOUBLE PRECISION COMMON PARAMETER THAT IS 25. C INCREMENTED BY THE NUMBER OF OPERATIONS 26. C PERFORMED BY THE SUBROUTINE. 27. C 28. C WORKING PARAMETERS - 29. C LINK - AT STEP J, THE LIST IN 30. C LINK(J), LINK(LINK(J)), ........... 31. C CONSISTS OF THOSE COLUMNS THAT WILL MODIFY 32. C THE COLUMN L(*,J). 33. C FIRST - TEMPORARY VECTOR TO POINT TO THE FIRST 34. C NONZERO IN EACH COLUMN THAT WILL BE USED 35. C NEXT FOR MODIFICATION. 36. C TEMP - A TEMPORARY VECTOR TO ACCUMULATE MODIFICATIONS. 37. C 38. C*************************************************************** 39. C 40. SUBROUTINE GSFCT ( NEQNS, XLNZ, LNZ, XNZSUB, NZSUB, DIAG, 41. 1 LINK, FIRST, TEMP, IFLAG ) 42. C 43. C*************************************************************** 44. C 45. DOUBLE PRECISION COUNT, OPS 46. COMMON /SPKOPS/ OPS 47. REAL DIAG(1), LNZ(1), TEMP(1), DIAGJ, LJK 48. INTEGER LINK(1), NZSUB(1) 49. INTEGER FIRST(1), XLNZ(1), XNZSUB(1), 50. 1 I, IFLAG, II, ISTOP, ISTRT, ISUB, J, 51. 1 K, KFIRST, NEQNS, NEWK 52. C 53. C*************************************************************** 54. C 55. C ------------------------------ 56. C INITIALIZE WORKING VECTORS ... 57. C ------------------------------ 58. DO 100 I = 1, NEQNS 59. LINK(I) = 0 60. TEMP(I) = 0.0E0 61. 100 CONTINUE 62. C -------------------------------------------- 63. C COMPUTE COLUMN L(*,J) FOR J = 1,...., NEQNS. 64. C -------------------------------------------- 65. DO 600 J = 1, NEQNS 66. C ------------------------------------------- 67. C FOR EACH COLUMN L(*,K) THAT AFFECTS L(*,J). 68. C ------------------------------------------- 69. DIAGJ = 0.0E0 70. NEWK = LINK(J) 71. 200 K = NEWK 72. IF ( K .EQ. 0 ) GO TO 400 73. NEWK = LINK(K) 74. C --------------------------------------- 75. C OUTER PRODUCT MODIFICATION OF L(*,J) BY 76. C L(*,K) STARTING AT FIRST(K) OF L(*,K). 77. C --------------------------------------- 78. KFIRST = FIRST(K) 79. LJK = LNZ(KFIRST) 80. DIAGJ = DIAGJ + LJK*LJK 81. OPS = OPS + 1.0D0 82. ISTRT = KFIRST + 1 83. ISTOP = XLNZ(K+1) - 1 84. IF ( ISTOP .LT. ISTRT ) GO TO 200 85. C ------------------------------------------ 86. C BEFORE MODIFICATION, UPDATE VECTORS FIRST, 87. C AND LINK FOR FUTURE MODIFICATION STEPS. 88. C ------------------------------------------ 89. FIRST(K) = ISTRT 90. I = XNZSUB(K) + (KFIRST-XLNZ(K)) + 1 91. ISUB = NZSUB(I) 92. LINK(K) = LINK(ISUB) 93. LINK(ISUB) = K 94. C --------------------------------------- 95. C THE ACTUAL MOD IS SAVED IN VECTOR TEMP. 96. C --------------------------------------- 97. DO 300 II = ISTRT, ISTOP 98. ISUB = NZSUB(I) 99. TEMP(ISUB) = TEMP(ISUB) + LNZ(II)*LJK 100. I = I + 1 101. 300 CONTINUE 102. COUNT = ISTOP - ISTRT + 1 103. OPS = OPS + COUNT 104. GO TO 200 105. C ---------------------------------------------- 106. C APPLY THE MODIFICATIONS ACCUMULATED IN TEMP TO 107. C COLUMN L(*,J). 108. C ---------------------------------------------- 109. 400 DIAGJ = DIAG(J) - DIAGJ 110. IF ( DIAGJ .LE. 0.0E0 ) GO TO 700 111. DIAGJ = SQRT(DIAGJ) 112. DIAG(J) = DIAGJ 113. ISTRT = XLNZ(J) 114. ISTOP = XLNZ(J+1) - 1 115. IF ( ISTOP .LT. ISTRT ) GO TO 600 116. FIRST(J) = ISTRT 117. I = XNZSUB(J) 118. ISUB = NZSUB(I) 119. LINK(J) = LINK(ISUB) 120. LINK(ISUB) = J 121. DO 500 II = ISTRT, ISTOP 122. ISUB = NZSUB(I) 123. LNZ(II) = ( LNZ(II)-TEMP(ISUB) ) / DIAGJ 124. TEMP(ISUB) = 0.0E0 125. I = I + 1 126. 500 CONTINUE 127. COUNT = ISTOP - ISTRT + 1 128. OPS = OPS + COUNT 129. 600 CONTINUE 130. RETURN 131. C ------------------------------------------------------ 132. C ERROR - ZERO OR NEGATIVE SQUARE ROOT IN FACTORIZATION. 133. C ------------------------------------------------------ 134. 700 IFLAG = 1 135. RETURN 136. END 137. C*************************************************************** C*************************************************************** C**** GENMMD ..... MULTIPLE MINIMUM EXTERNAL DEGREE **** C*************************************************************** C*************************************************************** C C AUTHOR - JOSEPH W.H. LIU C DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. C C PURPOSE - THIS ROUTINE IMPLEMENTS THE MINIMUM DEGREE C ALGORITHM. IT MAKES USE OF THE IMPLICIT REPRESENTATION C OF ELIMINATION GRAPHS BY QUOTIENT GRAPHS, AND THE C NOTION OF INDISTINGUISHABLE NODES. IT ALSO IMPLEMENTS C THE MODIFICATIONS BY MULTIPLE ELIMINATION AND MINIMUM C EXTERNAL DEGREE. C --------------------------------------------- C CAUTION - THE ADJACENCY VECTOR ADJNCY WILL BE C DESTROYED. C --------------------------------------------- C C INPUT PARAMETERS - C NEQNS - NUMBER OF EQUATIONS. C (XADJ,ADJNCY) - THE ADJACENCY STRUCTURE. C DELTA - TOLERANCE VALUE FOR MULTIPLE ELIMINATION. C MAXINT - MAXIMUM MACHINE REPRESENTABLE (SHORT) INTEGER C (ANY SMALLER ESTIMATE WILL DO) FOR MARKING C NODES. C C OUTPUT PARAMETERS - C PERM - THE MINIMUM DEGREE ORDERING. C INVP - THE INVERSE OF PERM. C NOFSUB - AN UPPER BOUND ON THE NUMBER OF NONZERO C SUBSCRIPTS FOR THE COMPRESSED STORAGE SCHEME. C C WORKING PARAMETERS - C DHEAD - VECTOR FOR HEAD OF DEGREE LISTS. C INVP - USED TEMPORARILY FOR DEGREE FORWARD LINK. C PERM - USED TEMPORARILY FOR DEGREE BACKWARD LINK. C QSIZE - VECTOR FOR SIZE OF SUPERNODES. C LLIST - VECTOR FOR TEMPORARY LINKED LISTS. C MARKER - A TEMPORARY MARKER VECTOR. C C PROGRAM SUBROUTINES - C MMDELM, MMDINT, MMDNUM, MMDUPD. C C*************************************************************** C SUBROUTINE GENMMD ( NEQNS, XADJ, ADJNCY, INVP, PERM, 1 DELTA, DHEAD, QSIZE, LLIST, MARKER, 1 MAXINT, NOFSUB ) C C*************************************************************** C C INTEGER*2 ADJNCY(1), DHEAD(1) , INVP(1) , LLIST(1) , INTEGER*4 ADJNCY(1), DHEAD(1) , INVP(1) , LLIST(1) , 1 MARKER(1), PERM(1) , QSIZE(1) INTEGER*4 XADJ(1) INTEGER*4 DELTA , EHEAD , I , MAXINT, MDEG , 1 MDLMT , MDNODE, NEQNS , NEXTMD, NOFSUB, 1 NUM, TAG C C*************************************************************** C IF ( NEQNS .LE. 0 ) RETURN C C ------------------------------------------------ C INITIALIZATION FOR THE MINIMUM DEGREE ALGORITHM. C ------------------------------------------------ NOFSUB = 0 CALL MMDINT ( NEQNS, XADJ, ADJNCY, DHEAD, INVP, PERM, 1 QSIZE, LLIST, MARKER ) C C ---------------------------------------------- C NUM COUNTS THE NUMBER OF ORDERED NODES PLUS 1. C ---------------------------------------------- NUM = 1 C C ----------------------------- C ELIMINATE ALL ISOLATED NODES. C ----------------------------- NEXTMD = DHEAD(1) 100 CONTINUE IF ( NEXTMD .LE. 0 ) GO TO 200 MDNODE = NEXTMD NEXTMD = INVP(MDNODE) MARKER(MDNODE) = MAXINT INVP(MDNODE) = - NUM NUM = NUM + 1 GO TO 100 C 200 CONTINUE C ---------------------------------------- C SEARCH FOR NODE OF THE MINIMUM DEGREE. C MDEG IS THE CURRENT MINIMUM DEGREE; C TAG IS USED TO FACILITATE MARKING NODES. C ---------------------------------------- IF ( NUM .GT. NEQNS ) GO TO 1000 TAG = 1 DHEAD(1) = 0 MDEG = 2 300 CONTINUE IF ( DHEAD(MDEG) .GT. 0 ) GO TO 400 MDEG = MDEG + 1 GO TO 300 400 CONTINUE C ------------------------------------------------- C USE VALUE OF DELTA TO SET UP MDLMT, WHICH GOVERNS C WHEN A DEGREE UPDATE IS TO BE PERFORMED. C ------------------------------------------------- MDLMT = MDEG + DELTA EHEAD = 0 C 500 CONTINUE MDNODE = DHEAD(MDEG) IF ( MDNODE .GT. 0 ) GO TO 600 MDEG = MDEG + 1 IF ( MDEG .GT. MDLMT ) GO TO 900 GO TO 500 600 CONTINUE C ---------------------------------------- C REMOVE MDNODE FROM THE DEGREE STRUCTURE. C ---------------------------------------- NEXTMD = INVP(MDNODE) DHEAD(MDEG) = NEXTMD IF ( NEXTMD .GT. 0 ) PERM(NEXTMD) = - MDEG INVP(MDNODE) = - NUM NOFSUB = NOFSUB + MDEG + QSIZE(MDNODE) - 2 IF ( NUM+QSIZE(MDNODE) .GT. NEQNS ) GO TO 1000 C ---------------------------------------------- C ELIMINATE MDNODE AND PERFORM QUOTIENT GRAPH C TRANSFORMATION. RESET TAG VALUE IF NECESSARY. C ---------------------------------------------- TAG = TAG + 1 IF ( TAG .LT. MAXINT ) GO TO 800 TAG = 1 DO 700 I = 1, NEQNS IF ( MARKER(I) .LT. MAXINT ) MARKER(I) = 0 700 CONTINUE 800 CONTINUE CALL MMDELM ( MDNODE, XADJ, ADJNCY, DHEAD, INVP, 1 PERM, QSIZE, LLIST, MARKER, MAXINT, 1 TAG ) NUM = NUM + QSIZE(MDNODE) LLIST(MDNODE) = EHEAD EHEAD = MDNODE IF ( DELTA .GE. 0 ) GO TO 500 900 CONTINUE C ------------------------------------------- C UPDATE DEGREES OF THE NODES INVOLVED IN THE C MINIMUM DEGREE NODES ELIMINATION. C ------------------------------------------- IF ( NUM .GT. NEQNS ) GO TO 1000 CALL MMDUPD ( EHEAD, NEQNS, XADJ, ADJNCY, DELTA, MDEG, 1 DHEAD, INVP, PERM, QSIZE, LLIST, MARKER, 1 MAXINT, TAG ) GO TO 300 C 1000 CONTINUE CALL MMDNUM ( NEQNS, PERM, INVP, QSIZE ) RETURN C END C*************************************************************** C*************************************************************** C*** MMDINT ..... MULT MINIMUM DEGREE INITIALIZATION *** C*************************************************************** C*************************************************************** C C AUTHOR - JOSEPH W.H. LIU C DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. C C PURPOSE - THIS ROUTINE PERFORMS INITIALIZATION FOR THE C MULTIPLE ELIMINATION VERSION OF THE MINIMUM DEGREE C ALGORITHM. C C INPUT PARAMETERS - C NEQNS - NUMBER OF EQUATIONS. C (XADJ,ADJNCY) - ADJACENCY STRUCTURE. C C OUTPUT PARAMETERS - C (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. C QSIZE - SIZE OF SUPERNODE (INITIALIZED TO ONE). C LLIST - LINKED LIST. C MARKER - MARKER VECTOR. C C*************************************************************** C SUBROUTINE MMDINT ( NEQNS, XADJ, ADJNCY, DHEAD, DFORW, 1 DBAKW, QSIZE, LLIST, MARKER ) C C*************************************************************** C C INTEGER*2 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , INTEGER*4 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , 1 LLIST(1) , MARKER(1), QSIZE(1) INTEGER*4 XADJ(1) INTEGER*4 FNODE , NDEG , NEQNS , NODE C C*************************************************************** C DO 100 NODE = 1, NEQNS DHEAD(NODE) = 0 QSIZE(NODE) = 1 MARKER(NODE) = 0 LLIST(NODE) = 0 100 CONTINUE C ------------------------------------------ C INITIALIZE THE DEGREE DOUBLY LINKED LISTS. C ------------------------------------------ DO 200 NODE = 1, NEQNS NDEG = XADJ(NODE+1) - XADJ(NODE) + 1 FNODE = DHEAD(NDEG) DFORW(NODE) = FNODE DHEAD(NDEG) = NODE IF ( FNODE .GT. 0 ) DBAKW(FNODE) = NODE DBAKW(NODE) = - NDEG 200 CONTINUE RETURN C END C*************************************************************** C*************************************************************** C** MMDELM ..... MULTIPLE MINIMUM DEGREE ELIMINATION *** C*************************************************************** C*************************************************************** C C AUTHOR - JOSEPH W.H. LIU C DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. C C PURPOSE - THIS ROUTINE ELIMINATES THE NODE MDNODE OF C MINIMUM DEGREE FROM THE ADJACENCY STRUCTURE, WHICH C IS STORED IN THE QUOTIENT GRAPH FORMAT. IT ALSO C TRANSFORMS THE QUOTIENT GRAPH REPRESENTATION OF THE C ELIMINATION GRAPH. C C INPUT PARAMETERS - C MDNODE - NODE OF MINIMUM DEGREE. C MAXINT - ESTIMATE OF MAXIMUM REPRESENTABLE (SHORT) C INTEGER. C TAG - TAG VALUE. C C UPDATED PARAMETERS - C (XADJ,ADJNCY) - UPDATED ADJACENCY STRUCTURE. C (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. C QSIZE - SIZE OF SUPERNODE. C MARKER - MARKER VECTOR. C LLIST - TEMPORARY LINKED LIST OF ELIMINATED NABORS. C C*************************************************************** C SUBROUTINE MMDELM ( MDNODE, XADJ, ADJNCY, DHEAD, DFORW, 1 DBAKW, QSIZE, LLIST, MARKER, MAXINT, 1 TAG ) C C*************************************************************** C C INTEGER*2 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , INTEGER*4 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , 1 LLIST(1) , MARKER(1), QSIZE(1) INTEGER*4 XADJ(1) INTEGER*4 ELMNT , I , ISTOP , ISTRT , J , 1 JSTOP , JSTRT , LINK , MAXINT, MDNODE, 1 NABOR , NODE , NPV , NQNBRS, NXNODE, 1 PVNODE, RLMT , RLOC , RNODE , TAG , 1 XQNBR C C*************************************************************** C C ----------------------------------------------- C FIND REACHABLE SET AND PLACE IN DATA STRUCTURE. C ----------------------------------------------- MARKER(MDNODE) = TAG ISTRT = XADJ(MDNODE) ISTOP = XADJ(MDNODE+1) - 1 C ------------------------------------------------------- C ELMNT POINTS TO THE BEGINNING OF THE LIST OF ELIMINATED C NABORS OF MDNODE, AND RLOC GIVES THE STORAGE LOCATION C FOR THE NEXT REACHABLE NODE. C ------------------------------------------------------- ELMNT = 0 RLOC = ISTRT RLMT = ISTOP DO 200 I = ISTRT, ISTOP NABOR = ADJNCY(I) IF ( NABOR .EQ. 0 ) GO TO 300 IF ( MARKER(NABOR) .GE. TAG ) GO TO 200 MARKER(NABOR) = TAG IF ( DFORW(NABOR) .LT. 0 ) GO TO 100 ADJNCY(RLOC) = NABOR RLOC = RLOC + 1 GO TO 200 100 CONTINUE LLIST(NABOR) = ELMNT ELMNT = NABOR 200 CONTINUE 300 CONTINUE C ----------------------------------------------------- C MERGE WITH REACHABLE NODES FROM GENERALIZED ELEMENTS. C ----------------------------------------------------- IF ( ELMNT .LE. 0 ) GO TO 1000 ADJNCY(RLMT) = - ELMNT LINK = ELMNT 400 CONTINUE JSTRT = XADJ(LINK) JSTOP = XADJ(LINK+1) - 1 DO 800 J = JSTRT, JSTOP NODE = ADJNCY(J) LINK = - NODE IF ( NODE ) 400, 900, 500 500 CONTINUE IF ( MARKER(NODE) .GE. TAG .OR. 1 DFORW(NODE) .LT. 0 ) GO TO 800 MARKER(NODE) = TAG C --------------------------------- C USE STORAGE FROM ELIMINATED NODES C IF NECESSARY. C --------------------------------- 600 CONTINUE IF ( RLOC .LT. RLMT ) GO TO 700 LINK = - ADJNCY(RLMT) RLOC = XADJ(LINK) RLMT = XADJ(LINK+1) - 1 GO TO 600 700 CONTINUE ADJNCY(RLOC) = NODE RLOC = RLOC + 1 800 CONTINUE 900 CONTINUE ELMNT = LLIST(ELMNT) GO TO 300 1000 CONTINUE IF ( RLOC .LE. RLMT ) ADJNCY(RLOC) = 0 C -------------------------------------------------------- C FOR EACH NODE IN THE REACHABLE SET, DO THE FOLLOWING ... C -------------------------------------------------------- LINK = MDNODE 1100 CONTINUE ISTRT = XADJ(LINK) ISTOP = XADJ(LINK+1) - 1 DO 1700 I = ISTRT, ISTOP RNODE = ADJNCY(I) LINK = - RNODE IF ( RNODE ) 1100, 1800, 1200 1200 CONTINUE C -------------------------------------------- C IF RNODE IS IN THE DEGREE LIST STRUCTURE ... C -------------------------------------------- PVNODE = DBAKW(RNODE) IF ( PVNODE .EQ. 0 .OR. 1 PVNODE .EQ. (-MAXINT) ) GO TO 1300 C ------------------------------------- C THEN REMOVE RNODE FROM THE STRUCTURE. C ------------------------------------- NXNODE = DFORW(RNODE) IF ( NXNODE .GT. 0 ) DBAKW(NXNODE) = PVNODE IF ( PVNODE .GT. 0 ) DFORW(PVNODE) = NXNODE NPV = - PVNODE IF ( PVNODE .LT. 0 ) DHEAD(NPV) = NXNODE 1300 CONTINUE C ---------------------------------------- C PURGE INACTIVE QUOTIENT NABORS OF RNODE. C ---------------------------------------- JSTRT = XADJ(RNODE) JSTOP = XADJ(RNODE+1) - 1 XQNBR = JSTRT DO 1400 J = JSTRT, JSTOP NABOR = ADJNCY(J) IF ( NABOR .EQ. 0 ) GO TO 1500 IF ( MARKER(NABOR) .GE. TAG ) GO TO 1400 ADJNCY(XQNBR) = NABOR XQNBR = XQNBR + 1 1400 CONTINUE 1500 CONTINUE C ---------------------------------------- C IF NO ACTIVE NABOR AFTER THE PURGING ... C ---------------------------------------- NQNBRS = XQNBR - JSTRT IF ( NQNBRS .GT. 0 ) GO TO 1600 C ----------------------------- C THEN MERGE RNODE WITH MDNODE. C ----------------------------- QSIZE(MDNODE) = QSIZE(MDNODE) + QSIZE(RNODE) QSIZE(RNODE) = 0 MARKER(RNODE) = MAXINT DFORW(RNODE) = - MDNODE DBAKW(RNODE) = - MAXINT GO TO 1700 1600 CONTINUE C -------------------------------------- C ELSE FLAG RNODE FOR DEGREE UPDATE, AND C ADD MDNODE AS A NABOR OF RNODE. C -------------------------------------- DFORW(RNODE) = NQNBRS + 1 DBAKW(RNODE) = 0 ADJNCY(XQNBR) = MDNODE XQNBR = XQNBR + 1 IF ( XQNBR .LE. JSTOP ) ADJNCY(XQNBR) = 0 C 1700 CONTINUE 1800 CONTINUE RETURN C END C*************************************************************** C*************************************************************** C***** MMDUPD ..... MULTIPLE MINIMUM DEGREE UPDATE ***** C*************************************************************** C*************************************************************** C C AUTHOR - JOSEPH W.H. LIU C DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. C C PURPOSE - THIS ROUTINE UPDATES THE DEGREES OF NODES C AFTER A MULTIPLE ELIMINATION STEP. C C INPUT PARAMETERS - C EHEAD - THE BEGINNING OF THE LIST OF ELIMINATED C NODES (I.E., NEWLY FORMED ELEMENTS). C NEQNS - NUMBER OF EQUATIONS. C (XADJ,ADJNCY) - ADJACENCY STRUCTURE. C DELTA - TOLERANCE VALUE FOR MULTIPLE ELIMINATION. C MAXINT - MAXIMUM MACHINE REPRESENTABLE (SHORT) C INTEGER. C C UPDATED PARAMETERS - C MDEG - NEW MINIMUM DEGREE AFTER DEGREE UPDATE. C (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. C QSIZE - SIZE OF SUPERNODE. C LLIST - WORKING LINKED LIST. C MARKER - MARKER VECTOR FOR DEGREE UPDATE. C TAG - TAG VALUE. C C*************************************************************** C SUBROUTINE MMDUPD ( EHEAD, NEQNS, XADJ, ADJNCY, DELTA, 1 MDEG, DHEAD, DFORW, DBAKW, QSIZE, 1 LLIST, MARKER, MAXINT, TAG ) C C*************************************************************** C C INTEGER*2 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , INTEGER*4 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , 1 LLIST(1) , MARKER(1), QSIZE(1) INTEGER*4 XADJ(1) INTEGER*4 DEG , DEG0 , DELTA , EHEAD , ELMNT , 1 ENODE , FNODE , I , IQ2 , ISTOP , 1 ISTRT , J , JSTOP , JSTRT , LINK , 1 MAXINT, MDEG , MDEG0 , MTAG , NABOR , 1 NEQNS , NODE , Q2HEAD, QXHEAD, TAG C C*************************************************************** C MDEG0 = MDEG + DELTA ELMNT = EHEAD 100 CONTINUE C ------------------------------------------------------- C FOR EACH OF THE NEWLY FORMED ELEMENT, DO THE FOLLOWING. C (RESET TAG VALUE IF NECESSARY.) C ------------------------------------------------------- IF ( ELMNT .LE. 0 ) RETURN MTAG = TAG + MDEG0 IF ( MTAG .LT. MAXINT ) GO TO 300 TAG = 1 DO 200 I = 1, NEQNS IF ( MARKER(I) .LT. MAXINT ) MARKER(I) = 0 200 CONTINUE MTAG = TAG + MDEG0 300 CONTINUE C --------------------------------------------- C CREATE TWO LINKED LISTS FROM NODES ASSOCIATED C WITH ELMNT: ONE WITH TWO NABORS (Q2HEAD) IN C ADJACENCY STRUCTURE, AND THE OTHER WITH MORE C THAN TWO NABORS (QXHEAD). ALSO COMPUTE DEG0, C NUMBER OF NODES IN THIS ELEMENT. C --------------------------------------------- Q2HEAD = 0 QXHEAD = 0 DEG0 = 0 LINK = ELMNT 400 CONTINUE ISTRT = XADJ(LINK) ISTOP = XADJ(LINK+1) - 1 DO 700 I = ISTRT, ISTOP ENODE = ADJNCY(I) LINK = - ENODE IF ( ENODE ) 400, 800, 500 C 500 CONTINUE IF ( QSIZE(ENODE) .EQ. 0 ) GO TO 700 DEG0 = DEG0 + QSIZE(ENODE) MARKER(ENODE) = MTAG C ---------------------------------- C IF ENODE REQUIRES A DEGREE UPDATE, C THEN DO THE FOLLOWING. C ---------------------------------- IF ( DBAKW(ENODE) .NE. 0 ) GO TO 700 C --------------------------------------- C PLACE EITHER IN QXHEAD OR Q2HEAD LISTS. C --------------------------------------- IF ( DFORW(ENODE) .EQ. 2 ) GO TO 600 LLIST(ENODE) = QXHEAD QXHEAD = ENODE GO TO 700 600 CONTINUE LLIST(ENODE) = Q2HEAD Q2HEAD = ENODE 700 CONTINUE 800 CONTINUE C -------------------------------------------- C FOR EACH ENODE IN Q2 LIST, DO THE FOLLOWING. C -------------------------------------------- ENODE = Q2HEAD IQ2 = 1 900 CONTINUE IF ( ENODE .LE. 0 ) GO TO 1500 IF ( DBAKW(ENODE) .NE. 0 ) GO TO 2200 TAG = TAG + 1 DEG = DEG0 C ------------------------------------------ C IDENTIFY THE OTHER ADJACENT ELEMENT NABOR. C ------------------------------------------ ISTRT = XADJ(ENODE) NABOR = ADJNCY(ISTRT) IF ( NABOR .EQ. ELMNT ) NABOR = ADJNCY(ISTRT+1) C ------------------------------------------------ C IF NABOR IS UNELIMINATED, INCREASE DEGREE COUNT. C ------------------------------------------------ LINK = NABOR IF ( DFORW(NABOR) .LT. 0 ) GO TO 1000 DEG = DEG + QSIZE(NABOR) GO TO 2100 1000 CONTINUE C -------------------------------------------- C OTHERWISE, FOR EACH NODE IN THE 2ND ELEMENT, C DO THE FOLLOWING. C -------------------------------------------- ISTRT = XADJ(LINK) ISTOP = XADJ(LINK+1) - 1 DO 1400 I = ISTRT, ISTOP NODE = ADJNCY(I) LINK = - NODE IF ( NODE .EQ. ENODE ) GO TO 1400 IF ( NODE ) 1000, 2100, 1100 C 1100 CONTINUE IF ( QSIZE(NODE) .EQ. 0 ) GO TO 1400 IF ( MARKER(NODE) .GE. TAG ) GO TO 1200 C ------------------------------------- C CASE WHEN NODE IS NOT YET CONSIDERED. C ------------------------------------- MARKER(NODE) = TAG DEG = DEG + QSIZE(NODE) GO TO 1400 1200 CONTINUE C ---------------------------------------- C CASE WHEN NODE IS INDISTINGUISHABLE FROM C ENODE. MERGE THEM INTO A NEW SUPERNODE. C ---------------------------------------- IF ( DBAKW(NODE) .NE. 0 ) GO TO 1400 IF ( DFORW(NODE) .NE. 2 ) GO TO 1300 QSIZE(ENODE) = QSIZE(ENODE) + 1 QSIZE(NODE) QSIZE(NODE) = 0 MARKER(NODE) = MAXINT DFORW(NODE) = - ENODE DBAKW(NODE) = - MAXINT GO TO 1400 1300 CONTINUE C -------------------------------------- C CASE WHEN NODE IS OUTMATCHED BY ENODE. C -------------------------------------- IF ( DBAKW(NODE) .EQ.0 ) 1 DBAKW(NODE) = - MAXINT 1400 CONTINUE GO TO 2100 1500 CONTINUE C ------------------------------------------------ C FOR EACH ENODE IN THE QX LIST, DO THE FOLLOWING. C ------------------------------------------------ ENODE = QXHEAD IQ2 = 0 1600 CONTINUE IF ( ENODE .LE. 0 ) GO TO 2300 IF ( DBAKW(ENODE) .NE. 0 ) GO TO 2200 TAG = TAG + 1 DEG = DEG0 C --------------------------------- C FOR EACH UNMARKED NABOR OF ENODE, C DO THE FOLLOWING. C --------------------------------- ISTRT = XADJ(ENODE) ISTOP = XADJ(ENODE+1) - 1 DO 2000 I = ISTRT, ISTOP NABOR = ADJNCY(I) IF ( NABOR .EQ. 0 ) GO TO 2100 IF ( MARKER(NABOR) .GE. TAG ) GO TO 2000 MARKER(NABOR) = TAG LINK = NABOR C ------------------------------ C IF UNELIMINATED, INCLUDE IT IN C DEG COUNT. C ------------------------------ IF ( DFORW(NABOR) .LT. 0 ) GO TO 1700 DEG = DEG + QSIZE(NABOR) GO TO 2000 1700 CONTINUE C ------------------------------- C IF ELIMINATED, INCLUDE UNMARKED C NODES IN THIS ELEMENT INTO THE C DEGREE COUNT. C ------------------------------- JSTRT = XADJ(LINK) JSTOP = XADJ(LINK+1) - 1 DO 1900 J = JSTRT, JSTOP NODE = ADJNCY(J) LINK = - NODE IF ( NODE ) 1700, 2000, 1800 C 1800 CONTINUE IF ( MARKER(NODE) .GE. TAG ) 1 GO TO 1900 MARKER(NODE) = TAG DEG = DEG + QSIZE(NODE) 1900 CONTINUE 2000 CONTINUE 2100 CONTINUE C ------------------------------------------- C UPDATE EXTERNAL DEGREE OF ENODE IN DEGREE C STRUCTURE, AND MDEG (MIN DEG) IF NECESSARY. C ------------------------------------------- DEG = DEG - QSIZE(ENODE) + 1 FNODE = DHEAD(DEG) DFORW(ENODE) = FNODE DBAKW(ENODE) = - DEG IF ( FNODE .GT. 0 ) DBAKW(FNODE) = ENODE DHEAD(DEG) = ENODE IF ( DEG .LT. MDEG ) MDEG = DEG 2200 CONTINUE C ---------------------------------- C GET NEXT ENODE IN CURRENT ELEMENT. C ---------------------------------- ENODE = LLIST(ENODE) IF ( IQ2 .EQ. 1 ) GO TO 900 GO TO 1600 2300 CONTINUE C ----------------------------- C GET NEXT ELEMENT IN THE LIST. C ----------------------------- TAG = MTAG ELMNT = LLIST(ELMNT) GO TO 100 C END C*************************************************************** C*************************************************************** C***** MMDNUM ..... MULTI MINIMUM DEGREE NUMBERING ***** C*************************************************************** C*************************************************************** C C AUTHOR - JOSEPH W.H. LIU C DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. C C PURPOSE - THIS ROUTINE PERFORMS THE FINAL STEP IN C PRODUCING THE PERMUTATION AND INVERSE PERMUTATION C VECTORS IN THE MULTIPLE ELIMINATION VERSION OF THE C MINIMUM DEGREE ORDERING ALGORITHM. C C INPUT PARAMETERS - C NEQNS - NUMBER OF EQUATIONS. C QSIZE - SIZE OF SUPERNODES AT ELIMINATION. C C UPDATED PARAMETERS - C INVP - INVERSE PERMUTATION VECTOR. ON INPUT, C IF QSIZE(NODE)=0, THEN NODE HAS BEEN MERGED C INTO THE NODE -INVP(NODE); OTHERWISE, C -INVP(NODE) IS ITS INVERSE LABELLING. C C OUTPUT PARAMETERS - C PERM - THE PERMUTATION VECTOR. C C*************************************************************** C SUBROUTINE MMDNUM ( NEQNS, PERM, INVP, QSIZE ) C C*************************************************************** C C INTEGER*2 INVP(1) , PERM(1) , QSIZE(1) INTEGER*4 INVP(1) , PERM(1) , QSIZE(1) INTEGER*4 FATHER, NEQNS , NEXTF , NODE , NQSIZE, 1 NUM , ROOT C C*************************************************************** C DO 100 NODE = 1, NEQNS NQSIZE = QSIZE(NODE) IF ( NQSIZE .LE. 0 ) PERM(NODE) = INVP(NODE) IF ( NQSIZE .GT. 0 ) PERM(NODE) = - INVP(NODE) 100 CONTINUE C ------------------------------------------------------ C FOR EACH NODE WHICH HAS BEEN MERGED, DO THE FOLLOWING. C ------------------------------------------------------ DO 500 NODE = 1, NEQNS IF ( PERM(NODE) .GT. 0 ) GO TO 500 C ----------------------------------------- C TRACE THE MERGED TREE UNTIL ONE WHICH HAS C NOT BEEN MERGED, CALL IT ROOT. C ----------------------------------------- FATHER = NODE 200 CONTINUE IF ( PERM(FATHER) .GT. 0 ) GO TO 300 FATHER = - PERM(FATHER) GO TO 200 300 CONTINUE C ----------------------- C NUMBER NODE AFTER ROOT. C ----------------------- ROOT = FATHER NUM = PERM(ROOT) + 1 INVP(NODE) = - NUM PERM(ROOT) = NUM C ------------------------ C SHORTEN THE MERGED TREE. C ------------------------ FATHER = NODE 400 CONTINUE NEXTF = - PERM(FATHER) IF ( NEXTF .LE. 0 ) GO TO 500 PERM(FATHER) = - ROOT FATHER = NEXTF GO TO 400 500 CONTINUE C ---------------------- C READY TO COMPUTE PERM. C ---------------------- DO 600 NODE = 1, NEQNS NUM = - INVP(NODE) INVP(NODE) = NUM PERM(NUM) = NODE 600 CONTINUE RETURN C END