! ! W.Quapp, Feb 2010 ! ! energy + gradient + Hessian from function mod nfk ! Programm bestimmt GE ! ! nn = Dimension ! le = Chainlength ! pstl = pred-Steplength ! eps = Convergence in RGF corrector ! rsel in corrector = tang INTEGER nn,le,iterNTmax DOUBLE PRECISION pstl PARAMETER (nn=2,le=24,iterNTmax=50,pstl=0.15d0) DOUBLE PRECISION kette(le+1,nn),eps,tgt,tt,energy,x,y, 1 rsel(nn),grad(nn),H(nn,nn),Proj(nn,nn), 2 xstart(nn),sw(nn),norm,dstep(nn),detk,t2(nn),t3(nn), 4 projNm1(nn-1,nn),dk(nn,nn+1),rgeq(nn-1,nn), 5 rgrad(nn-1),tang(nn),told(nn),unit(nn,nn),vec_ddot,det, 7 t1old(nn),t1(nn),rgnorm,point(nn),eival(nn),eivec(nn,nn) INTEGER k,j,i,jnr,istatus,icalc,iterNT logical omat_inv ccccccccccccccccccccccccccccccccccccccccc OPEN(8,FILE='dim2kette.dat') OPEN(7,FILE='chain.dat') OPEN(11,FILE='dim2dir.dat') OPEN(53,FILE='dim2Eall.dat') OPEN(25,FILE='dim2told.dat') OPEN(26,FILE='dim2t1old.dat') c protocol file: all interesting output of the vri-search OPEN(44,FILE='dim2protocol.txt') cccccccccccccccccccccccccccccccccccccccc 2 FORMAT(2F20.14) 50 FORMAT(1X,3F20.10) N=nn jnr=1 iterNT=0 istatus=2 Write(44,*) write(44,*)' Search TASC solution' write(*,*) ' Search TASC solution' eps=1.0d-10 c eps is criterion for rgf convergence to given rsel c jnr is the currend node number write(*,*) 'eps = ',eps write(44,*)'eps = ',eps xstart(1)= 2.85d0 xstart(2)= -0.1125d0 do i=1,N point(i)=xstart(i) kette(1,i)=xstart(i) enddo c other examples: (0.10, -0.10), (1.0, -1.30), (1.80, -2.20) WRITE(44,50)(xstart(I),I=1,N) WRITE(*,50) (xstart(I),I=1,N) WRITE(44,*) 'ketten-point start is ' WRITE(44,*) '{', xstart(1), ',', xstart(2) ,' },' WRITE(7,*) '{' WRITE(7,*) '{', xstart(1), ',', xstart(2) ,' },' call mat_diag(unit,nn,1.0d0) rsel(1)=-1.0d0 rsel(2)=-1.25d0 tgt=norm(rsel,nn) call vec_dmult_constant(rsel,nn,1.0d0/tgt,rsel) write(44,*) 'r_selected ' write(44,2)(rsel(i),i=1,nn) call vec_dmult_constant(rsel,nn,1.0d0,told) call vec_dmult_constant(rsel,nn,1.0d0,t1old) call vec_dmult_constant(rsel,nn,1.0d0,t2) call vec_dmult_constant(rsel,nn,1.0d0,t3) goto 11 c cccccccccccccccccccccccccccccccccccccccccccccccc 1 continue c cccccccccccccccccccccccccccccccccccccccccccccccc do i=1,nn point(i)=kette(jnr,i) enddo write(44,*) ' point ' WRITE(44,2) (point(I),I=1,nn) 11 x=point(1) y=point(2) c 0.03- modified NFK test surface, see CPL 2004 energy=-9.0d0*Dexp(-(-3.0d0+x)**2 - y**2) 1 -9.0d0*Dexp(-( 3.0d0+x)**2 - y**2) 2 + x*y +0.03d0*(x**2 + y**2)**2 write(44,*)'Energy ',energy write(53,*) energy c rewind 14 c Read(14,*) (grad(i),i=1,nn) grad(1)=18.0d0*(Dexp(-(-3.0d0+x)**2-y**2)*(-3.+x) + 1 Dexp(-( 3.0d0+x)**2-y**2)*( 3.+x))+ y+ 2 0.12d0*x*(x**2 + y**2) grad(2)= 1 18.0d0*y*(Dexp(-(-3.0d0+x)**2-y**2)+Dexp(-(3.0d0+x)**2-y**2)) 1 + x + 0.12d0*(x**2 + y**2)*y c rewind 15 c Do 15 k=1,nn c15 Read(15,*) (H(k,i),i=1,nn) h(1,1)= 1 18.0d0*(Dexp(-(-3.0d0+x)**2-y**2)+Dexp(-(3.0d0+x)**2-y**2))- 2 36.0d0*(Dexp(-(-3.0d0+x)**2-y**2)*(-3.0d0+ x)+ 3 Dexp(-( 3.0d0+x)**2 -y**2)*(3.0d0+ x)) + 4 0.12d0*(3.0d0*x**2 + y**2) h(1,2)= 1 - 36.0d0*Dexp(-(-3.0d0+x)**2-y**2)*(-3.0d0+x)*y+ 2 1.d0 + 0.24d0*x*y 3 - 36.0d0*Dexp(-( 3.0d0+x)**2-y**2)*( 3.0d0+x)*y h(2,2)= 1 18.0d0*(Dexp(-(-3.0d0+x)**2 -y**2)+Dexp(-(3.0d0+x)**2-y**2))- 2 36.d0*y**2*(Dexp(-(-3.d0+x)**2-y**2)+Dexp(-(3.d0+x)**2-y**2))+ 3 0.12d0*( x**2 + 3.0d0*y**2) h(2,1)= h(1,2) write(44,*) ' Hessian ' write(44,2) ((H(I,K),I=1,nn),K=1,nn) detHe=det(H,nn) write(44,*) 'det of Hessian', detHe write(*,*) 'det of Hessian', detHe c use Proj for a help matrix call eigen(H,nn,nn,eival,eivec,Proj,0) write(44,*)'eigenvalues of H' write(44,*)(eival(i),i=1,nn) write(44,*)'eigenvectors of H' do k=1,nn write(44,*) 'vec(',k,') = ',(eivec(k,i),i=1,nn) enddo write(44,*) ' grad ' WRITE(44,2) (grad(I),I=1,nn) rgnorm=norm(grad,nn) WRITE(44,*) ' full gradient norm ',rgnorm cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc new: projNm1: n-1 linearly independent lines from Proj call ortcompl(rsel,projNm1,unit,nn) write(44,*) ' projNm1 ' Do 13, k=1,nn-1 13 WRITE(44,2) (projNm1(k,I),I=1,nn) cccccccccccccccccccccccccccccccccccccccccccccccccccc c test print : is o.k. if(jnr.eq.1) then write(44,*) ' rsel ' WRITE(44,2) (rsel(I),I=1,nn) call vec_dinit(sw,nn,0.0d0) call matmult(projNm1,nn-1,rsel,1,sw,nn) write(44,*) ' Proj * rsel ' WRITE(44,*) (sw(I),I=1,nn) call vec_dinit(sw,nn,0.0d0) call matmult(projNm1,nn-1,grad,1,sw,nn) write(44,*) ' Proj * grad ' WRITE(44,*) (sw(I),I=1,nn) endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccc call matmult(projNm1,nn-1,H,nn,rgeq,nn) write(44,*) 'print rgeq ' Do 14, k=1,nn-1 14 WRITE(44,2) (rgeq(k,I),I=1,nn) call vec_dinit(rgrad,nn-1,0.0d0) call matmult(projNm1,nn-1,grad,1,rgrad,nn) write(44,*) ' rgrad ' WRITE(44,2) (rgrad(I),I=1,nn-1) rgnorm=norm(rgrad,nn-1) WRITE(44,*) ' PES reduced grad norm ',rgnorm c Tangent ccccccccccccccccccccccccccccccccccccccccccccc call gettang(rgeq,tang,nn) write(44,*) 'tang from rgeq is ' WRITE(44,2) (tang(I),I=1,nn) c ccccccccccccccccccccccccccccccccccccccccccccc call matmult(rgeq,nn-1,tang,1,sw,nn) write(44,*) ' test: rgeq*tang gives ' write(44,2) (sw(I),I=1,nn) tgt=vec_ddot(tang,told,nn) write(44,*)' tang * told ',tgt if (tgt.lt.0.0d0) 1 call vec_dmult_constant(tang,nn,-1.0d0,tang) cccccccccccccccccccccccccccccccccccccccccccccc call buildk(dk,rgeq,tang,rgrad,nn) cccccccccccccccccccccccccccccccccccccccccccccc write(44,*) ' K mat is ' Do 71, k=1,nn 71 WRITE(44,*) (dk(k,I),I=1,nn+1) detk=det(dk,nn) write(*,*) ' detK is ', detk write(44,*) ' detK is ', detk if (detk.lt.0.0d0) then call vec_dmult_constant(tang,nn,-1.0d0,t1) else call vec_dcopy(tang,t1,nn) endif tg=vec_ddot(t1,t1old,nn) if (tg.lt.0.0d0 .and. jnr.gt. le/2 ) then write(*,*) ' bif of branches' write(44,*) ' bif of branches' c to jump over turning points of NTs call vec_dmult_constant(tang,nn,-1.0d0,tang) call buildk(dk,rgeq,tang,rgrad,nn) endif call linsolve(dk,dstep,nn) WRITE(44,*) 'hypo Corr - dstep by linsolve ' WRITE(44,2)(dstep(i),i=1,nn) ccccccc 1.order goto: 333 continue ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccc Abbruch IF( 1 rgnorm.lt.eps .or. ! Corrector o.k. 3 iterNT.ge.iterNTmax ) THEN write(44,*) write(*,*) write(44,*)' Predictor step ---------------- pstl=',pstl write(*,*) ' Predictor step ---------------- pstl=',pstl write(44,*) 'Nr j in vri2 is =',jnr write(*,*) ' |rgnorm| ', rgnorm write(44,*) ' |rgnorm| ', rgnorm iterNT=1 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(jnr.lt.le+1) THEN cc c do predictorstep c this line is TASC !!! ccccccccccccccccccc call vec_dcopy(tang,rsel,nn) c ccccccccccccccccccccccccccccccccccccccccc do i=1,nn point(i)=kette(jnr,i) enddo cc predictor step by rgf3 c dk(nn,nn+1)=pstl c call linsolve(dk,dstep,nn) c write(44,*)'pred-step by rgf3 ' c or: call vec_dmult_constant(tang,nn,pstl,dstep) write(44,*)'pred-step by tangent ' tgt=norm(dstep,nn) write(*,*) write(44,*)'pred-step length would be',tgt write(44,2)(dstep(i),i=1,nn) do i=1,nn point(i)=point(i)+dstep(i)*pstl/tgt enddo kette(jnr+1,1)=point(1) kette(jnr+1,2)=point(2) WRITE(44,*) 'ketten-point nach predictor is ' WRITE(44,22) kette(jnr+1,1), kette(jnr+1,2) WRITE(7,22) kette(jnr+1,1), kette(jnr+1,2) 22 Format(4H { ,F20.14,3H , ,F20.14,5H },) jnr=jnr+1 call vec_dcopy(t2,t3,nn) call vec_dcopy(tang,t2,nn) GOTO 66 c noch Abbruch else istatus=10 goto 66 ENDIF ccc AbbruchEnde ENDIF ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCCCCCCCCCCCCCCCCCC tt=norm(dstep,nn) write(44,*) ' dstep norm whould be ',tt write(44,*) ' dstep is ' ccccc Notbremse -- emergency break : c if(tt.gt.pstl) c 1 call vec_dmult_constant(dstep,nn,pstl/tt,dstep) write(44,2) (dstep(I),I=1,nn) DO 25 i=1,nn point(i)=point(i)+dstep(i) 25 CONTINUE do i=1,nn kette(jnr,i)=point(i) enddo WRITE(44,*) 'ketten-point nach dstep addition' WRITE(44,22) kette(jnr,1),kette(jnr,2) WRITE(7,22) kette(jnr,1),kette(jnr,2) write(44,*) 'Nr j in vri2 is =',jnr write(*,*) 'j vri2 = ',jnr,' iterNT = ',iterNT write(*,*) ' hypo |dstep| ', tt,' rg ',rgnorm write(44,*) 'hypo |dstep| ', tt,' rg ',rgnorm iterNT=iterNT+1 cccccccccccccccccccccccccccccccccccccccccccccccccccccc 66 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccc WRITE(44,*) jnr WRITE(44,2) (kette(jnr,i),i=1,nn) rewind 8 DO 29 k=1,le+1 29 WRITE(8,2) (kette(k,i),i=1,nn) c if(istatus.lt.10) goto 1 c WRITE(7,23) kette(jnr,1),kette(jnr,2) 23 Format(4H { ,F20.14,3H , ,F20.14,5H }}) END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccc cccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine ginvers(in,NAT,out,omat_inv) c calculates the inverse matrix of IN real*8 in,out logical omat_inv integer NAT,i,j,k dimension in(NAT,NAT),out(NAT,NAT) do 1,k=1,NAT do 1,j=1,NAT 1 out(j,k)=in(j,k) c do 10,k=1,NAT if (out(k,k).eq.0.0d0) then omat_inv=.false. return endif out(k,k)=1.0d0/out(k,k) do 20,i=1,NAT if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k)) 20 continue do 30,j=1,NAT if (j.ne.k) then do 40,i=1,NAT if (i.ne.k) out(i,j)=out(i,j)+out(i,k)*out(k,j) 40 continue out(k,j)=out(k,j)*out(k,k) endif 30 continue 10 continue omat_inv=.true. return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine mat_diag(d,n,v) c create a NxN diagonal matrix D with entries V real*8 d,v integer n,i dimension d(n,n) if(n .le. 0) return call vec_dinit(d,n*n,0.0d0) do i=1,n d(i,i)=v end do return end cccccccccccccccccccccccccccccccc ! ! Multiplikation Ax=b ! SUBROUTINE multiply(b,A,x,nn) DOUBLE PRECISION b(nn),A(nn,nn),x(nn) DO 20 i=1,nn b(i)=0 DO 10 j=1,nn b(i)=b(i)+A(i,j)*x(j) 10 CONTINUE 20 CONTINUE END ! ! Gibt die Norm von v aus ! DOUBLE PRECISION FUNCTION norm(v,nn) DOUBLE PRECISION v(nn),s s=0.d0 DO 10 j=1,nn s=s+v(j)**2 10 CONTINUE IF (s.LT.1.D-23) THEN norm=0.d0 ELSE norm=DSQRT(s) ENDIF END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine ortcompl(vec,compl,dmet,ndim) c calculate an Orthonormalcomplement COMPL to vector VEC real*8 basis,vec,compl,swap,det,dmgs,dmet integer i,j,ndim dimension vec(ndim),compl(ndim-1,ndim),basis(ndim,ndim), * swap(ndim),dmet(ndim,ndim) call vec_dinit(basis,ndim*ndim,0.0d0) do 10,i=1,ndim basis(i,i)=1.0d0 10 continue do 20,i=1,ndim if (vec(i).ne.0.0d0) then do 30,j=1,ndim swap(j)=basis(j,1) basis(j,1)=vec(j) if (i.ne.1) basis(j,i)=swap(j) 30 continue goto 40 endif 20 continue 40 continue c note: dmgs is the Gram-Schmidt method det=dmgs(basis,dmet,ndim) call vec_dinit(compl,(ndim-1)*ndim,0.0d0) call mat_trans(basis,ndim,ndim,basis) do 60,i=2,ndim call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim) 60 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function dmgs(a,dmet,dim) c modified Gram-Schmidt algorithm return det(a) real*8 det,rr,a,q,r,vec_mmdot,dmet,dzero integer dim,i,k dimension a(dim,dim),q(dim,dim),r(dim,dim),dmet(dim,dim) data dzero/0.0d0/ call vec_dinit(q,dim*dim,0.0d0) call vec_dinit(r,dim*dim,0.0d0) det=1.0d0 call mat_dcopy(a,q,dim,dim) do 10,k=1,dim if (k.eq.1) goto 21 do 20,i=1,k-1 r(i,k)=vec_mmdot(i,q,dim,k,q,dim,dim,dmet) call vec_madd(k,q,dim,1.0d0,i,q,dim,-r(i,k),k,q,dim,dim) 20 continue 21 continue rr=vec_mmdot(k,q,dim,k,q,dim,dim,dmet) r(k,k)=dsqrt(rr) det=det*r(k,k) if(rr.lt.dzero) goto 33 call vec_madd(k,q,dim,1.0d0/r(k,k),k,q,dim,0.0d0,k,q,dim,dim) 10 continue call mat_dcopy(q,a,dim,dim) dmgs=det return 33 CONTINUE cc call cerr(-23) return end ccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function vec_mmdot &(nri,matin,din,nra,matout,dout,dim,dmet) c multiply the (nri) col of (matin) and the (nra) col of c (matout). (dim) is the number of the rows with metric dmet real*8 matin,matout,prod,v1,v2,dmet,vec_mdot integer nri,din,nra,dout,dim dimension matin(dim,din),matout(dim,dout),dmet(dim,dim), + v1(dim),v2(dim) ccc if ((nri.gt.din).or.(nra.gt.dout)) call cerr(-100) call vec_mat_centry(v1,matin,nri,din,dim,'o') call vec_mat_centry(v2,matout,nra,dout,dim,'o') prod=vec_mdot(v1,dmet,v2,dim) vec_mmdot=prod return end cccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function vec_mdot(vec1,metric,vec2,dim) c metric scalar product of VEC1 and VEC2 with respect to METRIC integer dim real*8 vec1,metric,vec2,cov2,vec_ddot dimension vec1(dim),vec2(dim),cov2(dim),metric(dim,dim) call matmult(metric,dim,vec2,1,cov2,dim) vec_mdot=vec_ddot(vec1,cov2,dim) return end ccccccccccccccccccccccccccccccccccccccccccccccc subroutine gettang(curveq,tangent,nn) c This routine calculate the tangent on a curve c whereas the tangent is unique determined by: c c H*t=0 c |t|=1 c det[H,t]>0 c real*8 tangent,curveq,q,ss1,det,zero,norm integer nn,i,j dimension tangent(nn),curveq(nn-1,nn),q(nn*nn),ss1(nn,nn) dimension gmat(nn,nn) data zero /0.0d0/ call qr(curveq,q,nn-1) call vec_dcopy(q(nn*(nn-1)+1),tangent(1),nn) c test do 1,i=1,nn c 1 tangent(i)=q(nn,i) c wq call vec_dmult_constant + (tangent,nn,1.0d0/norm(tangent,nn),tangent) cc call mat_dcopy(curveq,ss1,nn-1,nn) do 10,j = 1,nn-1 do 20,i = 1,nn ss1(j,i) = curveq(j,i) 20 continue 10 continue do 30,i = 1,nn 30 ss1(nn,i)=tangent(i) if (det(ss1,nn).lt.zero) then call vec_dmult_constant(tangent,nn,-1.0d0,tangent) endif return end ccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine qr(gex,q,n) C QR-Zerlegung von gex integer i,j,k,n real*8 s1,s2,s,gex,q,r,qq,rr dimension gex(n,n+1),q(n+1,n+1),r(n+1,n), * qq(n+1,n+1),rr(n+1,n) do 10,i=1,n do 20,j=1,n+1 r(j,i)=gex(i,j) 20 continue 10 continue do 30,i=1,n+1 do 40,j=1,n+1 q(i,j)=0.0d0 40 continue q(i,i)=1.0d0 30 continue do 50,i=1,n do 60,j=i+1,n+1 s1=r(i,i) s2=r(j,i) s=dsqrt(s1*s1+s2*s2) if (s.gt.0.0d0) then s1=s1/s s2=s2/s else s1=0.0d0 s2=0.0d0 endif call vec_dcopy(r,rr,(n+1)*n) do 70,k=1,n rr(i,k)=s1*r(i,k)+s2*r(j,k) rr(j,k)=s1*r(j,k)-s2*r(i,k) 70 continue call vec_dcopy(rr,r,(n+1)*n) call vec_dcopy(q,qq,(n+1)*(n+1)) do 80,k=1,n+1 qq(i,k)=s1*q(i,k)+s2*q(j,k) qq(j,k)=s1*q(j,k)-s2*q(i,k) 80 continue call vec_dcopy(qq,q,(n+1)*(n+1)) 60 continue 50 continue do 90,i=1,n+1 do 100,j=1,n+1 q(i,j)=qq(j,i) 100 continue 90 continue c write(44,*) 'in QR ' c do i=1,n+1 c write(44,*) (q(i,j),j=1,n+1) c enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine tri_mat(dm1,nr1,nc1,dm2,nc2,dm3,nc3,dres) c dm1[nr1,nc1]*dm2[nc1,nc2]*dm3[nc2,nc3]=dres[nr1,nc3] real*8 dm1,dm2,dm3,dres,r1 integer nr1,nc1,nc2,nc3 dimension dm1(nr1,nc1),dm2(nc1,nc2),dm3(nc2,nc3), + r1(nr1,nc2),dres(nr1,nc3) call matmult(dm1,nr1,dm2,nc2,r1,nc1) call matmult(r1,nr1,dm3,nc3,dres,nc2) return end ccccccccccccccc subroutine mat_dcopy(arra1,arra2,ndim,mdim) c Copies matrix ARRA1 to matrix ARRA2 real*8 arra1,arra2 integer i,j,ndim,mdim dimension arra1(ndim,mdim),arra2(ndim,mdim) do 10,i = 1,ndim do 20,j = 1,mdim arra2(i,j) = arra1(i,j) 20 continue 10 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine buildk(dk,req,tang,rgr,nn) c build the K-Matrix real*8 dk,req,tang,rgr integer nn,i,j dimension dk(nn,nn+1),req(nn-1,nn),tang(nn),rgr(nn-1) do 10,i=1,nn-1 do 11,j=1,nn dk(i,j)=req(i,j) 11 continue 10 continue do 20,j=1,nn dk(nn,j)=tang(j) 20 continue do 30,i=1,nn-1 dk(i,nn+1)=-rgr(i) 30 continue dk(nn,nn+1)=0.0d0 return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C subroutine linsolve(matein,x,n) c solves MATEIN * X = 0 real*8 matein,dre,x,w integer n dimension matein(n,n+1),dre(n,n+1),x(n),w(n) call householder(matein,dre,w,n,n+1) c write(44,*) 'in linsolve ' c do k=1,n c write(44,*) (dre(k,i),i=1,n+1) c enddo call loesung(dre,x,n) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dcopy(arra1,arra2,ndim) c Copies ARRA1 to ARRA2 real*8 arra1,arra2 integer i,ndim dimension arra1(ndim),arra2(ndim) if(ndim .le. 0) return do 10 i = 1,ndim arra2(i) = arra1(i) 10 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccc subroutine householder(matein,mataus,w,n,sp) c householder transformation integer k,l,m,n,sp real*8 alpha,rho,sk,s,w,matein,mataus dimension matein(n,sp),mataus(n,sp),w(n) call vec_dcopy(matein,mataus,n*sp) do 10,k=1,n-1 s=0.0d0 do 20,l=k,n s=s+mataus(l,k)*mataus(l,k) 20 continue if (mataus(k,k).lt.0) then alpha=dsqrt(s) else alpha=-dsqrt(s) endif rho=dsqrt(s-mataus(k,k)*mataus(k,k)+ * (mataus(k,k)-alpha)*(mataus(k,k)-alpha)) if(rho .lt. 1.0d-63) then rho= 1.0d-63 w(k)= 0.0d0 goto 22 endif w(k)=(mataus(k,k)-alpha)/rho 22 continue do 30,l=k+1,n w(l)=mataus(l,k)/rho 30 continue do 40,l=k,sp sk=0.0d0 do 55,m=k,n sk=sk+w(m)*mataus(m,l) 55 continue do 60,m=k,n mataus(m,l)=mataus(m,l)-2*sk*w(m) 60 continue 40 continue 10 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine loesung(matein,x,n) c intern linear equation routine real*8 matein,x integer i,j,n,k dimension matein(n,n+1),x(n) do 10,i=1,n 10 x(i)=matein(i,n+1) x(n)=x(n)/matein(n,n) do 20,i=1,n-1 k=n-i do 30,j=k+1,n 30 x(k)=x(k)-matein(k,j)*x(j) x(k)=x(k)/matein(k,k) 20 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_mat_centry(vec,mat,npos,ndim,ncol,zkind) c Copies VEC 'i'n or 'o'ut the NPOS column of matrix MAT real*8 vec,mat integer npos,ndim,ncol,i character zkind dimension vec(ndim),mat(ndim,ncol) if(npos.gt.ncol) return if (zkind.eq.'i') then do i=1,ndim mat(i,npos)=vec(i) end do endif if (zkind.eq.'o') then do i=1,ndim vec(i)=mat(i,npos) enddo endif return end cccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function vec_ddot(arra1,arra2,ndim) c cartesian scalar product of ARRA1 and ARRA2 real*8 arra1,arra2 integer ndim,i dimension arra1(ndim),arra2(ndim) vec_ddot = 0 if(ndim .le. 0)return do 70 i = 1,ndim vec_ddot = vec_ddot + arra1(i) * arra2(i) 70 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dmult_add(arra1,arra2,ndim,zahl,arra4) c For each Komponent: ARRA4 = ARRA1 + ZAHL*ARRA2 real*8 arra1,arra2,arra4,zahl integer ndim,i dimension arra1(ndim),arra2(ndim),arra4(ndim) if(ndim .le. 0)return do 145 i = 1,ndim arra4(i) = arra1(i)+arra2(i)*zahl 145 continue return end ccccccccccccccccccccccccccccccccccccccc subroutine matmult(m1,d1,m2,d2,mout,dim) c matrix multiplication: m1[d1,dim]*m2[dim,d2]=mout[d1,d2] real*8 m1,m2,mout,swap integer d1,d2,dim,i,j,k dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2) call vec_dinit(swap,d1*d2,0.0d0) do 10,i=1,d1 do 20,j=1,d2 do 30,k=1,dim swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j) 30 continue 20 continue 10 continue call vec_dcopy(swap,mout,d1*d2) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dinit(arra1,ndim,wert) c Set all elements of ARRA1 to the double value of WERT real*8 arra1,wert integer ndim,i dimension arra1(ndim) if(ndim .le. 0) return do i=1,ndim arra1(i)=wert end do return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine mat_trans(matin,nrow,ncol,matout) c transpose MATOUT of matrix MATIN real*8 matin,matout,swap integer nrow,ncol,i,j dimension matin(nrow,ncol),matout(ncol,nrow),swap(ncol,nrow) do 10,i=1,nrow do 20,j=1,ncol swap(j,i)=matin(i,j) 20 continue 10 continue call mat_dcopy(swap,matout,ncol,nrow) return end ccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_mrowcopy(nri,matin,din,nra,matout,dout,dim) c copies the (nri) row of (matin) with (din) row and (dim) col c to the (nra) row of (matout) with (dout) row and (dim) col real*8 matin,matout integer nri,din,nra,dout,dim,i dimension matin(din,dim),matout(dout,dim) ccc if ((nri.gt.din).or.(nra.gt.dout)) call cerr(17) do 10,i=1,dim matout(nra,i)=matin(nri,i) 10 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_madd * (nr1,mat1,d1,fact1,nr2,mat2,d2,fact2,nrx,matx,dx,dim) c matx(nrx)=fact1*mat1(nr1)+fact2*mat2(nr2) c d1,d2,dx - columns c dim - rows real*8 mat1,mat2,matx,fact1,fact2,vec integer nr1,nr2,nrx,d1,d2,dx,dim,i dimension mat1(dim,d1),mat2(dim,d2),matx(dim,dx),vec(dim,1) call vec_dinit(vec,dim,0.0d0) do 10,i=1,dim vec(i,1)=fact1*mat1(i,nr1)+fact2*mat2(i,nr2) 10 continue call vec_mat_centry(vec,matx,nrx,dim,dx,'i') return end cccccccccccccccccccccccccccccccccccccccccccccccccc subroutine vec_dmult_constant(arra1,ndim,zahl,arra2) c ARRA2 = ZAHL * ARRA1 real*8 arra1,arra2,zahl integer ndim,i dimension arra1(ndim),arra2(ndim) if(ndim .le. 0)return do i = 1,ndim arra2(i) = arra1(i)*zahl end do return end ccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 function det(matein,nn) c calculation of the determinant of MATEIN integer i,nn real*8 matein,mataus,w dimension matein(nn,nn),mataus(nn,nn),w(nn) call householder(matein,mataus,w,nn,nn) det=1.0d0 do 10,i=1,nn if (dabs(mataus(i,i)).lt.1.0d-63) then det= 1.0d-63 goto 20 endif if (dabs(det) .lt.1.0d-63) then det= 1.0d-63 goto 20 endif 10 det=det*mataus(i,i) if (iand(nn,1).eq.0) det=-det 20 continue return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine eigen (a,m,n,d,vec,e,iff) c eigenvalues and -vectors of A real*8 a,d,vec,e,eps,tol,h,g,s,f,b,p,r,c integer m,iff,n,i,j,ni,ii,l,k,j1 dimension a(m,m), d(m), vec(m,m) dimension e(m) eps = 0.278d-16 tol = 0.211758d-21 if (n.eq.1) then d(1) = a(1,1) vec(1,1) = 1.0d0 else do 30 i = 1 , n c householder's reduction c simulation of loop do 150 i=n,2,(-1) do 20 j = 1 , i vec(i,j) = a(i,j) 20 continue 30 continue do 120 ni = 2 , n ii = n + 2 - ni do 110 i = ii , ii l = i - 2 h = 0.0d0 g = vec(i,i-1) if (l.gt.0) then do 40 k = 1 , l h = h + vec(i,k)**2 40 continue s = h + g*g if (s.lt.tol) then h = 0.0d0 else if (h.gt.0) then l = l + 1 f = g g = dsqrt(s) if (f.gt.0) then g = -g end if h = s - f*g vec(i,i-1) = f - g f = 0.0d0 do 70 j = 1 , l vec(j,i) = vec(i,j)/h s = 0.0d0 do 50 k = 1 , j s = s + vec(j,k)*vec(i,k) 50 continue j1 = j + 1 if (j1.le.l) then do 60 k = j1 , l s = s + vec(k,j)*vec(i,k) 60 continue end if e(j) = s/h f = f + s*vec(j,i) 70 continue f = f/(h+h) do 80 j = 1 , l e(j) = e(j) - f*vec(i,j) 80 continue do 100 j = 1 , l f = vec(i,j) s = e(j) do 90 k = 1 , j vec(j,k) = vec(j,k) - f*e(k) - vec(i,k)*s 90 continue 100 continue end if end if c accumulation of transformation matrices d(i) = h e(i-1) = g 110 continue 120 continue d(1) = vec(1,1) vec(1,1) = 1.0d0 do 170 i = 2 , n l = i - 1 if (d(i).gt.0) then do 150 j = 1 , l s = 0.0d0 do 130 k = 1 , l s = s + vec(i,k)*vec(k,j) 130 continue do 140 k = 1 , l vec(k,j) = vec(k,j) - s*vec(k,i) 140 continue 150 continue end if d(i) = vec(i,i) vec(i,i) = 1.0d0 do 160 j = 1 , l c diagonalization of the tridiagonal matrix vec(i,j) = 0.0d0 vec(j,i) = 0.0d0 160 continue 170 continue b = 0.0d0 f = 0.0d0 e(n) = 0.0d0 do 250 l = 1 , n c test for splitting h = eps*(dabs(d(l))+dabs(e(l))) if (h.gt.b) b = h do 180 j = l , n c test for convergence if (dabs(e(j)).le.b) go to 190 180 continue 190 if (j.eq.l) then d(l) = d(l) + f else 200 p = (d(l+1)-d(l))*0.5d0/e(l) r = dsqrt(p*p+1.0d0) if (p.lt.0) then p = p - r else p = p + r end if h = d(l) - e(l)/p do 210 i = l , n c qr transformation d(i) = d(i) - h 210 continue f = f + h p = d(j) c simulation of loop do 330 i=j-1,l,(-1) c = 1.0d0 s = 0.0d0 j1 = j - 1 do 240 ni = l , j1 ii = l + j1 - ni do 230 i = ii , ii c protection against underflow of exponents g = c*e(i) h = c*p if (dabs(p).lt.dabs(e(i))) then c = p/e(i) r = dsqrt(c*c+1.0d0) e(i+1) = s*e(i)*r s = 1.0d0/r c = c/r else c = e(i)/p r = dsqrt(c*c+1.0d0) e(i+1) = s*p*r s = c/r c = 1.0d0/r end if p = c*d(i) - s*g d(i+1) = h + s*(c*g+s*d(i)) do 220 k = 1 , n h = vec(k,i+1) vec(k,i+1) = vec(k,i)*s + h*c vec(k,i) = vec(k,i)*c - h*s 220 continue 230 continue 240 continue e(l) = s*p c convergence d(l) = c*p if (dabs(e(l)).gt.b) go to 200 c ordering of eigenvalues d(l) = d(l) + f end if 250 continue if (iff.lt.1) then ni = n - 1 do 280 i = 1 , ni k = i p = d(i) j1 = i + 1 do 260 j = j1 , n if (d(j).lt.p) then k = j p = d(j) end if 260 continue if (k.ne.i) then d(k) = d(i) d(i) = p do 270 j = 1 , n p = vec(j,i) vec(j,i) = vec(j,k) vec(j,k) = p 270 continue end if c special treatment of case n = 1 280 continue end if end if return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! subroutine mat_adj(matin,ndim,matout) c adjoint matrix of matin real*8 matin,matout,matscr,det integer ndim,i,j,k,l,ii,jj dimension matin(ndim,ndim),matout(ndim,ndim), 1 matscr(ndim-1,ndim-1) call vec_dinit(matscr,(ndim-1)*(ndim-1),0.0d0) do 10, i=1 , ndim do 20, j=1 , ndim do 30, ii=1 , ndim-1 do 40, jj=1 , ndim-1 if (ii.lt.i) then k=ii else k=ii+1 endif if (jj.lt.j) then l=jj else l=jj+1 endif matscr(ii,jj)=matin(k,l) if(dabs(matscr(ii,jj)).lt.1.0d-63) then matscr(ii,jj)=1.0d-63 endif 40 continue 30 continue matout(i,j)=det(matscr,ndim-1) if(dabs(matout(i,j)).lt.1.0d-70) goto 20 matout(i,j)= (-1.0d0)**(i+j)*matout(i,j) 20 continue 10 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccc