C ************************** IMPORTANT ***************************** C * C All software in this file is (C) Copyright Marco Marletta, 1992. * C All rights reserved. * C * C ******************************************************************* SUBROUTINE sl12f(elam1,eps,a,b,k,n,tol,hamat,sysbc1,xmesh,mesh, & work,iwork,lwk,ilwk,nxtrap,ifail) C This subroutine solves eigenvalue problems for linear Hamiltonian C systems in which the sub-matrix S22 is positive definite. The C problems concerned are reducible to the form of a first order matrix C differential equation C Theta' = i Theta Omega(x,lambda,Theta) C where Theta is a unitary matrix which must satisfy prescribed C boundary conditions at two endpoints x=a and x=b. C The method use by SL12F is to approximate the coefficient matrices in C the Hamiltonian system by piecewise constants and transform the problem C to one in Sturm Liouville form using a transformation due to Reid (1973). C C Variables: C elam1: double precision. On entry elam1 should be set to an approximation C to the eigenvalue to be found. The value of elam1 is not important in C determining the accuracy of the final result but it may affect the C run times. On exit, elam1 specifies the approximation to the eigenvalue C computed by SL13F. C eps: double precision. On entry eps should be set to a positive real C number of an order of magnitude such that the true eigenvalue lies C within [elam1-C*eps,elam1+C*eps], where by elam1 we mean the initial C approximation provided by the user, and C is a constant of modest size. C Again the value of eps is not critical. On exit, eps specifies a C (generally pessimistic) assessment of the error in the final C approximation elam1. C a,b: double precision. On entry, a and b specify the finite regular C endpoints of the interval on which the problem is posed. Unchanged C on exit. C k: integer. The index of the eigenvalue sought. For many problems k=0 C specifies the lowest eigenvalue (e.g. for problems derived from C formally self-adjoint even-order differential operators). However C for other problems there exist eigenvalues for any integer value of C the index k, positive or negative. Unchanged on exit. C n: integer. The dimension of the system, i.e. the size of the matrix C Theta. Unchanged on exit. C tol: double precision. On entry tol should have a positive value. C By a process of Wynn extrapolation, SL13F will endeavour to ensure C that the error in the final approximation elam1 is less than tol. C Unchanged on exit. C xmesh: double precision array of dimension (0:mesh). Used as workspace. C mesh: integer. On entry, the first dimension of xmesh as declared in the C calling (sub)program. The user should assign mesh a large enough value C to ensure that SL12F is able to store a sufficiently fine mesh to meet C the user's accuracy requirements. For most purposes a mesh of 2000 C intervals is adequate (i.e. setting mesh = 2000 should be sufficient). C On exit, mesh specifies the number of mesh intervals actually used. C work: double precision array of dimension iwork. Used as workspace. C iwork: integer. The first dimension of work as declared in the calling C (sub)program, iwork >= n*(12+64*n). Unchanged on exit. C lwk: integer array of dimension ilwk. Used as workspace. C ilwk: integer. The first dimension of lwk as declared in the calling C (sub)program, ilwk >= n. Unchanged on exit. C nxtrap: integer. On entry, nxtrap need not be set. On exit, nxtrap C specifies the number of Wynn extrapolations which were required C to achieve the user's specified accuracy. C ifail: integer. On entry, ifail must have one of the values -1,0,1. C On successful exit, ifail will have the value 0. Any other value C indicates a failure. C C User-supplied subroutines HAMAT and SYSBC1: C C Hamat: C C SUBROUTINE hamat(x,elam,s11,s12,s21,s22,n) C INTEGER n C DOUBLE PRECISION x,elam,s11(n,n),s12(n,n),s21(n,n),s22(n,n) C C This routine must specify in s11, s12, s21 and s22 the values of the C four coefficient matrix functions of x and lambda (elam) occurring C in the Hamiltonian system. It must not change the values of x, elam, C or n. C C Sysbc1: C C SUBROUTINE sysbc1(iend,ising,x,u,v,ndim,n,elam) C INTEGER iend,ndim,n C LOGICAL ising C DOUBLE PRECISION x,elam,u(ndim,n),v(ndim,n) C C If iend=0 on entry then sysbc1 must specify the boundary condition at C x=a, while if iend=1 then sysbc1 must specify the boundary condition C at x=b. The parameter ising must be assigned the value .false. before C returning. The parameters x, elam, ndim and n must not be changed by C sysbc1. The specification of boundary conditions is achieved as follows: C a boundary condition C transpose(A1).U + transpose(A2).V = 0 C is specified by setting C U = A2, V = -A1 C and returning. C C END of brief documentation. C C .. Parameters .. INTEGER ismax,ismall,nsampl PARAMETER (ismax=10,ismall=300,nsampl=20) C .. C .. Scalar Arguments .. DOUBLE PRECISION a,b,elam1,eps,tol INTEGER ifail,ilwk,iwork,k,mesh,n,nxtrap C .. C .. Array Arguments .. DOUBLE PRECISION work(((24+2*nsampl)*n+12)*n),xmesh(0:mesh) INTEGER lwk(ilwk) C .. C .. Local Scalars .. DOUBLE PRECISION c,el,epso,err,eu,h,toloc INTEGER i,ic,ifo,imatch,istor,j,kloc,kntmax,kp,newmsh,nmesh,nrec, & nsmall CHARACTER srname*6 C .. C .. Local Arrays .. DOUBLE PRECISION argsmp(nsampl),dsampl(nsampl),elam(ismax), & xsmall(0:ismall) INTEGER nj(ismax) CHARACTER rec(3)*80 C .. C .. External Functions .. DOUBLE PRECISION x02ajf INTEGER p01abf EXTERNAL x02ajf,p01abf C .. C .. External Subroutines .. EXTERNAL chosc2,extrap,hamesh,solvsy C .. C .. Intrinsic Functions .. INTRINSIC abs,max,sqrt,float C .. C .. Subroutine Arguments .. EXTERNAL hamat,sysbc1 C .. srname = 'SL12F' ifo = ifail C First do the input checks. IF (a.GE.b.OR.n.LT.2.OR.tol.LE.0.0D0. & OR.iwork.LT.n*(12+64*n).OR. & ilwk.LT.n.OR.eps.LE.0.0D0) THEN C We have an input error ifail = 1 IF (a.GE.b) THEN WRITE (REC,FMT=9000) a,b 9900 FORMAT ('** A must be less than B, you have',/, & '** A = ',g18.8,' B = ',g18.8) nrec = 2 GOTO 70 END IF IF (n.LT.2) THEN WRITE (REC,FMT=9920) n 9920 FORMAT ('** N must be at least 2, you have',/, & '** N = ',i8) nrec = 2 GOTO 70 END IF IF (tol.LE.0.0D0) THEN WRITE (REC,FMT=9930) tol 9930 FORMAT ('** TOL must be strictly positive, you have',/, & '** TOL = ',g18.8) nrec = 2 GOTO 70 END IF IF (iwork.LT.n*(12+64*n)) THEN WRITE (REC,FMT=9940) & n*(12+64*n),iwork 9940 FORMAT ('** IWORK must be at least ',i8,' you have',/, & '** IWORK = ',i8) nrec = 2 GOTO 70 END IF IF (ilwk.LT.n) THEN WRITE (REC,FMT=9950) n,ilwk 9950 FORMAT ('** ILWK must be at least',i8,' you have',/, & '** ILWK = ',i8) nrec = 2 GOTO 70 END IF IF (eps.LE.0.0D0) THEN WRITE (REC,FMT=9955) eps 9955 FORMAT ('** EPS must be strictly positive, you have',/, & '** EPS = ',g18.8) nrec = 2 GOTO 70 END IF END IF nxtrap = 0 n2 = n*n iul = 1 ivl = iul+n2 iur = ivl+n2 ivr = iur+n2 iwk = ivr+n2 isampl = iwk+20*n2+12*n epso = eps err = epso c = 0.50d0* (a+b) toloc = min(1.0D-2,sqrt(tol)) ifail = 0 C Set up a mesh using a fairly ad-hoc method CALL hamesh(a,b,c,n,elam1,xmesh,mesh,nmesh,ic,toloc,work,hamat, & ifail) IF (ifail.NE.0) THEN C There are two cases: C IFAIL = 10: mesh refined as much as possible and still cannot meet C accuracy requirement. C IFAIL = 11: mesh routine ran out of space. IF (ifail.EQ.11) THEN WRITE (rec,FMT=9000) 9000 FORMAT ('** Meshing routine ran out of space,',/, & '** increase MESH or increase TOL.') ifail = 5 ELSE WRITE (rec,FMT=9010) 9010 FORMAT ('** Meshing routine making little progress,',/, & '** stepsize too small. Increase TOL.') ifail = 4 END IF GO TO 70 END IF C If there are too many mesh points then we simply thin the mesh out by C the crudest method possible. 10 IF (nmesh.GT.ismall) THEN newmsh = nmesh/2 ic = ic/2 DO 20 j = 1,newmsh xmesh(j) = xmesh(2*j) 20 CONTINUE xmesh(newmsh) = xmesh(nmesh) nmesh = newmsh GO TO 10 END IF imatch = ic nj(1) = 1 nj(2) = 2 nj(3) = 3 nj(4) = 4 nj(5) = 6 nj(6) = 8 nj(7) = 10 nj(8) = 14 nj(9) = 18 nj(10) = 24 C nj(1) = 1 C nj(2) = 2 C nj(3) = 4 C nj(4) = 8 C nj(5) = 16 C nj(6) = 32 C nj(7) = 64 C nj(8) = 128 C nj(9) = 256 C nj(10) = 512 C Store the mesh in the small array: DO 30 j = 0,nmesh xsmall(j) = xmesh(j) 30 CONTINUE nsmall = nmesh istor = 1 C Solve the problem 40 ifail = 0 toloc = sqrt(x02ajf(0.0D0)) kntmax = 100 eps = max(epso/100.D0,10.0D0*tol,ABS(err)) CALL solvsy(xmesh,nmesh,ic,k,elam1,eps,el,eu,hamat,sysbc1,toloc, & work,lwk,n,kntmax,ifail) IF (ifail.NE.0) THEN C Some useful messages for the user: nrec = 2 IF (ifail.EQ.6) THEN ifail = 2 WRITE (rec,FMT=9020) 9020 FORMAT ('** Error finding Theta matrix from boundary',/, & '** conditions supplied.') ELSE IF (ifail.EQ.7) THEN ifail = 3 WRITE (rec,FMT=9030) 9030 FORMAT ('** Too many miss-distance evaluations.') nrec = 1 ELSE IF (ifail.EQ.14) THEN ifail = 6 WRITE (rec,FMT=9040) 9040 FORMAT ('** The matrix V+iU has become rank deficient',/, & '** Try reducing TOL.') ELSE IF (ifail.EQ.16) THEN ifail = 8 WRITE (rec,FMT=9050) 9050 FORMAT ('** The eigenvalues of a Theta matrix cannot',/, & '** be found, possibly due to ill-conditioning of', & /,'** V+iU. Try reducing TOL.') nrec = 3 ELSE IF (ifail.EQ.17) THEN ifail = 9 WRITE (rec,FMT=9060) 9060 FORMAT ('** An H-matrix has occured whose eigenvalues',/, & '** cannot be computed accurately.') ELSE IF (ifail.EQ.18) THEN ifail = 10 WRITE (rec,FMT=9070) 9070 FORMAT ('** A non-Hermitian H-matrix has occurred',/, & '** possibly due to roundoff error.') ELSE IF (ifail.EQ.19) THEN ifail = 7 WRITE (rec,FMT=9080) 9080 FORMAT ('** The matrix V+iU has become ill-conditioned') nrec = 1 ELSE IF (ifail.EQ.20) THEN ifail = 11 WRITE (rec,FMT=9090) 9090 FORMAT ('** The coefficient S22 cannot be diagonalised',/, & '** with sufficient accuracy.') ELSE IF (ifail.EQ.21) THEN ifail = 12 WRITE (rec,FMT=9100) 9100 FORMAT ('** The coefficient S22 is not positive definite.' & ) nrec = 1 END IF GO TO 70 END IF IF (istor.EQ.1) THEN C Now get a better matchpoint by calling the new matchpoint routine. CALL chosc2(a,b,ic,work(iul),work(ivl),n,work(iur),work(ivr), & n,el,eu,xmesh,nmesh,hamat,sysbc1,work(iwk),lwk, & n,k,nsampl,work(isampl),dsampl,argsmp, & ifail) imatch = ic END IF C WRITE (6,*) ' nmesh,ic:',nmesh,ic C WRITE (6,FMT=*) 'Eigenvalue found was',elam1 elam(istor) = elam1 istor = istor + 1 err = epso IF (istor.GT.2) THEN C We have some eigenvalues available and we can do some Richardson C extrapolation: CALL extrap(elam,work,work(istor),nj,istor-1,elam1,err) C WRITE (6,FMT=*) 'Best approx so far:',elam1 C WRITE (6,FMT=*) 'Current error:',err END IF IF (abs(err).LT.tol*max(1.0D0,abs(elam1)) .OR. istor.GT.ismax.OR. & 2*nmesh.GT.mesh) THEN C WRITE (6,FMT=*) 'Number of meshpoints used:',nmesh C WRITE (6,FMT=*) 'Number of extrapolations:',istor - 2 nxtrap = istor - 2 eps = err mesh = nmesh RETURN END IF C Now add more meshpoints and see what happens nmesh = nsmall*nj(istor) ic = imatch*nj(istor) kp = 1 DO 60 i = 1,nsmall h = (xsmall(i)-xsmall(i-1))/float(nj(istor)) DO 50 kloc = 1,nj(istor) xmesh(kp) = xmesh(kp-1) + h kp = kp + 1 50 CONTINUE 60 CONTINUE GO TO 40 70 ifail = p01abf(ifo,ifail,srname,nrec,rec) RETURN END C ---------------------------------------------------------------------- SUBROUTINE extrap(a,b,c,nj,n,best,err) C .. Scalar Arguments .. DOUBLE PRECISION best,err INTEGER n C .. C .. Array Arguments .. DOUBLE PRECISION a(n),b(n),c(n,n) INTEGER nj(n) C .. C .. Local Scalars .. DOUBLE PRECISION fac1 INTEGER j,k,p C .. C .. Intrinsic Functions .. INTRINSIC float C .. DO 10 j = 1,n b(j) = a(j) 10 CONTINUE C DO 30 k = 1,n - 1 C fac1 = 2.0D0**(2*k) C DO 20 j = n,k + 1,-1 C b(j) = (b(j)-b(j-1)/fac1)/ (1.0D0-1.0D0/fac1) C20 CONTINUE C30 CONTINUE C DO 30 p = 1,n DO 20 j = 1,n c(j,p) = (1.0D0/float(nj(j)))** (2*p) 20 CONTINUE 30 CONTINUE DO 70 k = 1,n - 1 DO 40 j = n,k + 1,-1 best = (b(j)- (c(j,1)/c(j-1,1))*b(j-1))/ & (1.0D0-c(j,1)/c(j-1,1)) IF (j.EQ.n.AND.k.EQ.n-1) err = best-b(j) b(j) = best 40 CONTINUE DO 60 p = n - k,1,-1 DO 50 j = n,k + 1,-1 c(j,p) = (c(j,p+1)- (c(j,1)/c(j-1,1))*c(j-1,p+1))/ & (1.0D0-c(j,1)/c(j-1,1)) 50 CONTINUE 60 CONTINUE 70 CONTINUE best = b(n) RETURN END C ----------------------------------------------------------------- SUBROUTINE hamesh(a,b,c,n,elam,xmesh,maxmsh,nmesh,imatch,tol,work, & hamat,ifail) C .. Scalar Arguments .. DOUBLE PRECISION a,b,c,elam,tol INTEGER ifail,imatch,maxmsh,n,nmesh C .. C .. Array Arguments .. DOUBLE PRECISION work(12*n*n),xmesh(0:maxmsh) C .. C .. Subroutine Arguments .. EXTERNAL hamat C .. C .. Local Scalars .. INTEGER ia1,ia2,ia3,ib1,ib2,ib3,ic1,ic2,ic3,id1,id2,id3 C .. C .. External Subroutines .. EXTERNAL hamish C .. ia1 = 1 ia2 = ia1 + n*n ia3 = ia2 + n*n ib1 = ia3 + n*n ib2 = ib1 + n*n ib3 = ib2 + n*n ic1 = ib3 + n*n ic2 = ic1 + n*n ic3 = ic2 + n*n id1 = ic3 + n*n id2 = id1 + n*n id3 = id2 + n*n CALL hamish(a,b,c,n,elam,xmesh,maxmsh,nmesh,imatch,tol,work(ia1), & work(ia2),work(ia3),work(ib1),work(ib2),work(ib3), & work(ic1),work(ic2),work(ic3),work(id1),work(id2), & work(id3),hamat,ifail) RETURN END C ---------------------------------------------------------------------- SUBROUTINE hamish(a,b,c,n,elam,xmesh,maxmsh,nmesh,imatch,tol,a1, & a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3,hamat,ifail) C This subroutine forms an initial mesh for a Hamiltonian system using C the Pruess equidistribution heuristic. C C Equidistribution heuristic: (h**3)*max(coefficient'') = const. C C .. Parameters .. DOUBLE PRECISION half,two,one,safe,trd,five,ten,safe3,p008,sixth, & zero,tenth,twelf PARAMETER (half=5.d-1,two=2.d0,one=1.d0,safe=0.8D0,trd=one/3.d0, & five=5.d0,ten=10.d0,safe3=7.29d-1,p008=8.d-3, & sixth=one/6.d0,zero=0.d0,tenth=one/ten,twelf=half*sixth) C .. C .. Scalar Arguments .. DOUBLE PRECISION a,b,c,elam,tol INTEGER ifail,imatch,maxmsh,n,nmesh C .. C .. Array Arguments .. DOUBLE PRECISION a1(n,n),a2(n,n),a3(n,n),b1(n,n),b2(n,n),b3(n,n), & c1(n,n),c2(n,n),c3(n,n),d1(n,n),d2(n,n),d3(n,n), & xmesh(0:maxmsh) C .. C .. Local Scalars .. DOUBLE PRECISION di,err,h,hmax,qmax,ratio,tolow,x1,x2,x3,xend,xo INTEGER i,ic,icnt,idi,imesh,irf,j,kntr,ii,jj LOGICAL phase1,repeat C .. C .. External Functions .. DOUBLE PRECISION x02ajf EXTERNAL x02ajf C .. C .. External Subroutines .. C .. C .. Intrinsic Functions .. INTRINSIC abs,log,max,min,sign,int C .. C .. Subroutine Arguments .. EXTERNAL hamat C .. tolow = safe*tol irf = int(log(x02ajf(1.0D0))*0.50D0/log(safe)) C First do the left hand half of the interval: xo = c/ (10.0D0+abs(c)) xend = a/ (10.0D0+abs(a)) imesh = maxmsh x1 = xo xmesh(imesh) = x1 imesh = imesh - 1 h = abs(xend-xo)/20 di = -1.0D0 idi = -1 hmax = h kntr = 0 icnt = 0 phase1 = .true. repeat = .false. 10 x2 = x1 + h*di IF (x2*di.GT.xend*di) x2 = xend C Get the error measure for the interval [x1,x2]: C CALL hamat(10.D0*x1/ (1.0D0-abs(x1)),elam,a1,b1,c1,d1,n) CALL hamat(10.0D0*x1/ (1.0D0-abs(x1)),elam,c1,d1,a1,b1,n) C CALL hamat(10.D0*x2/ (1.0D0-abs(x2)),elam,a2,b2,c2,d2,n) CALL hamat(10.0D0*x2/ (1.0D0-abs(x2)),elam,c2,d2,a2,b2,n) x3 = 0.50D0* (x1+x2) C CALL hamat(10.D0*x3/(1.0D0-abs(x3)),elam,a3,b3,c3,d3,n) CALL hamat(10.0D0*x3/ (1.0D0-abs(x3)),elam,c3,d3,a3,b3,n) DO 7000 jj = 1,n DO 7010 ii = 1,n c1(ii,jj)=-c1(ii,jj) c2(ii,jj)=-c2(ii,jj) c3(ii,jj)=-c3(ii,jj) 7010 CONTINUE 7000 CONTINUE C Form the three-point difference in the usual way: err = 0.0D0 qmax = 1.0D0 DO 30 j = 1,n DO 20 i = 1,n qmax = max(qmax,abs(a3(i,j)),abs(b3(i,j)),abs(c3(i,j)), & abs(d3(i,j))) a3(i,j) = a1(i,j) + a2(i,j) - 2.0D0*a3(i,j) b3(i,j) = b1(i,j) + b2(i,j) - 2.0D0*b3(i,j) c3(i,j) = c1(i,j) + c2(i,j) - 2.0D0*c3(i,j) d3(i,j) = d1(i,j) + d2(i,j) - 2.0D0*d3(i,j) err = max(err,abs(a3(i,j)),abs(b3(i,j)),abs(c3(i,j)), & abs(d3(i,j))) 20 CONTINUE 30 CONTINUE err = ((1.0D0-abs(x3))**2)*h*err/sqrt(max(one,qmax)) DO 35 j = 1,n DO 25 i = 1,n c3(i,j) = 0.50D0*(c2(i,j)-c1(i,j)) d3(i,j) = 0.50D0*(d2(i,j)-d1(i,j)) 25 CONTINUE 35 CONTINUE CALL f06yaf('N','N',n,n,n,1.0D0,a3,n,c3,n,0.0D0,b1,n) CALL f06yaf('T','N',n,n,n,1.0D0,a3,n,c3,n,0.0D0,b2,n) CALL f06yaf('N','N',n,n,n,1.0D0,a3,n,d3,n,0.0D0,d1,n) CALL f06yaf('T','N',n,n,n,1.0D0,a3,n,d3,n,0.0D0,d2,n) CALL f06yaf('N','N',n,n,n,1.0D0,c3,n,d3,n,0.0D0,b3,n) qmax = 0.0D0 DO 37 j = 1,n DO 27 i = 1,n qmax = max(qmax,abs(b1(i,j)),abs(b2(i,j)),abs(b3(i,j)), & abs(d1(i,j)),abs(d2(i,j))) 27 CONTINUE 37 CONTINUE err = max(err,h**2*qmax/sqrt(max(one,qmax))) C Now decide what to do at the next step: IF (err.GT.tol) THEN repeat = .true. ratio = safe IF (tol.LT.err*safe3) ratio = (tol/err)**trd h = h*ratio icnt = icnt + 1 IF (icnt.LT.irf) GO TO 10 ifail = 10 RETURN END IF IF (err.LT.tolow .AND. icnt.EQ.0) THEN ratio = five IF (err.GT.tolow*p008) ratio = (tolow/err)**trd h = h*ratio h = min(h,hmax) C The next piece of code gets the meshing routine 'on-scale' at the first step. IF (phase1) THEN repeat = .true. kntr = kntr + 1 IF (kntr.LT.5) GO TO 10 END IF END IF C General step has been completed successfully, prepare for next step. phase1 = .false. xmesh(imesh) = x2 IF (x2*di.LT.xend*di) THEN imesh = imesh + idi IF (imesh.LT.0.0D0 .OR. imesh.GT.maxmsh) THEN ifail = 11 RETURN END IF x1 = x2 GO TO 10 C Look to see if we have done the whole of the interval yet: ELSE IF (di*sign(1.0D0,a-c).GT.0.0D0) THEN C Shift what we have already got of the mesh into the first half of the array C xmesh: ic = maxmsh - imesh imatch = ic IF (ic.GE.maxmsh) THEN ifail = 11 RETURN END IF DO 40 j = 0,ic xmesh(j) = xmesh(imesh+j) 40 CONTINUE xo = c/ (10.0D0+abs(c)) xend = b/ (10.0D0+abs(b)) x1 = xo imesh = ic + 1 di = 1.0D0 idi = 1 kntr = 0 icnt = 0 phase1 = .true. repeat = .false. GO TO 10 END IF nmesh = imesh C Convert mesh back to original coordinates DO 50 i = 0,nmesh xmesh(i) = 10.D0*xmesh(i)/ (1.0D0-abs(xmesh(i))) 50 CONTINUE ifail = 0 RETURN END C --------------------------------------------------------------------- SUBROUTINE chosc2(a,b,ic,ul,vl,ndiml,ur,vr,ndimr,el,eu,xmesh, & nmesh,cofmat,sysbc,work,lwk,n,k,nsampl,wsampl,dlsamp, & argsmp,ifail) C Subroutine to determine the best matching point for a shooting algorithm C to solve a Hamiltonian System eigenvalue problem. Uses a lot of storage, C unfortunately. C .. Parameters .. INTEGER iparam PARAMETER (iparam=1) C .. C .. Scalar Arguments .. DOUBLE PRECISION a,b,el,eu INTEGER ic,ifail,k,n,ndiml,ndimr,nmesh,nsampl C .. C .. Array Arguments .. DOUBLE PRECISION argsmp(nsampl),dlsamp(nsampl),ul(ndiml,n), & ur(ndimr,n),vl(ndiml,n),vr(ndimr,n), & work(n*(18*n+12+ndiml+ndimr)), & wsampl(2*nsampl*n**2),xmesh(0:nmesh) INTEGER lwk(n) C .. C .. Subroutine Arguments .. EXTERNAL cofmat,sysbc C .. C .. Local Scalars .. DOUBLE PRECISION alfaj,alfmax,alfmin,alfsum,argc,argl,argr,c, & diff,dla11f,elam,shift,twopi INTEGER i,imissi,imissr,iri,irr,isampi,isampr,ix,iy,j,jmatch, & jsampl,jstep,n2 LOGICAL ising,lnor C .. C .. Local Arrays .. DOUBLE PRECISION params(iparam) C .. C .. External Functions .. DOUBLE PRECISION x01aaf EXTERNAL x01aaf C .. C .. External Subroutines .. EXTERNAL f02ajf,f06qff,f06yaf,spthta,sumvec,thstp4 C .. C .. Intrinsic Functions .. INTRINSIC abs,atan2,dble C .. ifail = 0 shift = dble(k) + 5.0D-1 twopi = 2.0D0*x01aaf(0.0D0) n2 = n*n ic = nmesh/2 jstep = max(1,nmesh/nsampl) nsmplo = min(nsampl,nmesh) C WRITE (6,*) 'nmesh,nsmplo,jstep:',nmesh,nsmplo,jstep elam = el diff = 100.0D0 C Set up initial condition at x=a: 10 CALL sysbc(0,ising,a,ul,vl,ndiml,n,elam) C Here we have used ul and vl as storage for y and pdy respectively. C Get the phase-angles of the initial Theta-matrix: CALL f06qff('G',n,n,ul,ndiml,work(n+1),n) CALL f06qff('G',n,n,vl,ndiml,work(n2+n+1),n) lnor = .false. ifail = 0 CALL spthta(work(n+1),work(n2+n+1),n,n,work(1),work(2*n2+n+1), & work(3*n2+n+1),lwk,lnor,ifail) IF (ifail.NE.0) THEN ifail = 5 RETURN END IF CALL sumvec(work,n,argl) argl = -argl alfsum = argl C Do case of lower eigenvalue approximationm first: params(1) = elam C Integrate from x=a to x=b storing the matrix Sigma and the scalar argl C at each of the sampling points, which are equally spaced between C x=a and x=b: isampr = 1 isampi = n*n + 1 DO 30 jsampl = 1,nsmplo ifail = 0 DO 20 i = (jsampl-1)*jstep,jsampl*jstep - 1 C WRITE (7,*) 'Stepping from ',i,' to ',i+1 CALL thstp4(xmesh(i),xmesh(i+1),ul,vl,ndiml,n,elam,argl, & alfsum,work,lwk,cofmat,ifail) 20 CONTINUE IF (ifail.NE.0) RETURN C Store the outcome of the integration: CALL f06qff('G',n,n,ul,ndiml,wsampl(isampr),n) CALL f06qff('G',n,n,vl,ndiml,wsampl(isampi),n) argsmp(jsampl) = argl isampr = isampi + n*n isampi = isampr + n*n 30 CONTINUE C Set up initial condition at x=b: CALL sysbc(1,ising,b,ur,vr,ndimr,n,elam) C Here we have used ur and vr as storage for y and pdy respectively. C Get the phase-angles of the initial Theta-matrix: CALL f06qff('G',n,n,ur,ndimr,work(n+1),n) CALL f06qff('G',n,n,vr,ndimr,work(n2+n+1),n) lnor = .true. ifail = 0 CALL spthta(work(n+1),work(n2+n+1),n,n,work(1),work(2*n2+n+1), & work(3*n2+n+1),lwk,lnor,ifail) IF (ifail.NE.0) THEN ifail = 5 RETURN END IF CALL sumvec(work,n,argr) argr = -argr alfsum = 0.0D0 DO 25 i = 1,n IF (work(i).NE.twopi) alfsum = alfsum + work(i) 25 CONTINUE alfsum = -alfsum C Integrate from x=b to x=c: jmatch = nsmplo DO 40 i = nmesh,nsmplo*jstep + 1,-1 C WRITE (7,*) 'Stepping from ',i,' to ',i-1 CALL thstp4(xmesh(i),xmesh(i-1),ur,vr,ndimr,n,elam,argr, & alfsum,work,lwk,cofmat,ifail) 40 CONTINUE DO 70 jsampl = nsmplo - 1,1,-1 ifail = 0 DO 50 i = (jsampl+1)*jstep,jsampl*jstep + 1,-1 C WRITE (7,*) 'Stepping from ',i,' to ',i-1 CALL thstp4(xmesh(i),xmesh(i-1),ur,vr,ndimr,n,elam,argr, & alfsum,work,lwk,cofmat,ifail) 50 CONTINUE IF (ifail.NE.0) RETURN C C Form the product of the sigma matrices at the centre. What we actually want C is Conj(ThetaR)*ThetaL = SigmaR*Conj(SigmaL) C = (ur + i*vr)*(transp(ul)-i*transp(vl)). C isampi = 1+n2+2*(jsampl-1)*n2 isampr = isampi - n*n imissr = 1 CALL f06yaf('N','T',n,n,n,1.0D0,ur,ndimr,wsampl(isampr),n, & 0.0D0,work(imissr),n) CALL f06yaf('N','T',n,n,n,1.0D0,vr,ndimr,wsampl(isampi),n, & 1.0D0,work(imissr),n) imissi = imissr + n2 CALL f06yaf('N','T',n,n,n,1.0D0,vr,ndimr,wsampl(isampr),n, & 0.0D0,work(imissi),n) CALL f06yaf('N','T',n,n,n,-1.0D0,ur,ndimr,wsampl(isampi),n, & 1.0D0,work(imissi),n) C Get eigenvalues using F02AJF: irr = imissi + n2 iri = irr + n ifail = -1 CALL f02ajf(work(imissr),n,work(imissi),n,n,work(irr), & work(iri),lwk,ifail) IF (ifail.NE.0) THEN ifail = 4 RETURN END IF C Now get the phase angles which are required to form the miss distance C function. alfaj = 0.0D0 argc = 0.0D0 alfmax = -100.0D0 alfmin = 100.0D0 DO 60 j = 1,n ix = irr + j - 1 iy = iri + j - 1 alfaj = atan2(work(iy),work(ix)) IF (alfaj.LT.0.0D0) alfaj = alfaj + twopi IF (alfaj.GT.alfmax) alfmax = alfaj IF (alfaj.LT.alfmin) alfmin = alfaj C WRITE (6,FMT=*) 'Eigenvalue from F02AJF',work(ix),work(iy) argc = argc + alfaj 60 CONTINUE argl = argsmp(jsampl) dla11f = (-argl+argr-argc)/twopi + n - shift IF (dla11f.GT.0.0D0 .AND. dla11f.LT.0.6D0) dla11f = alfmin/ & twopi IF (dla11f.LT.0.0D0 .AND. dla11f.GT. & (-0.6D0)) dla11f = alfmax/twopi - 1.0D0 IF (elam.EQ.el) THEN dlsamp(jsampl) = dla11f ELSE C Compare with the miss-distance at the lower lambda and see if the difference C is smaller than at other mesh-points. C WRITE (6,FMT=*) 'x,du-dl:',xmesh(i-1), C & dla11f - dlsamp(jsampl) IF (abs(dla11f-dlsamp(jsampl)).LE.diff) THEN jmatch = jsampl diff = abs(dla11f-dlsamp(jsampl)) ic = i-1 END IF END IF 70 CONTINUE IF (elam.EQ.el) THEN elam = eu GO TO 10 END IF C WRITE (6,FMT=*) 'jmatch,nsmplo,a,b:',jmatch,nsmplo,a,b c = xmesh(ic) C ic = jmatch*jstep C WRITE (6,FMT=*) 'Match-point chosen (c,ic): ',c,ic RETURN END C --------------------------------------------------------------------- SUBROUTINE solvsy(xmesh,nmesh,ic,k,elam,eps,el,eu,hamat,sysbc,tol, & work,lwk,n,kntmax,ifail) C This subroutine solves a regular matrix-vector Sturm-Liouville system, C provided there are no bugs in it. C First find a bracket for the eigenvalue: C .. Scalar Arguments .. DOUBLE PRECISION el,elam,eps,eu,tol INTEGER ic,ifail,k,kntmax,n,nmesh C .. C .. Array Arguments .. DOUBLE PRECISION work(24*n**2+12*n),xmesh(0:nmesh) INTEGER lwk(n) C .. C .. Subroutine Arguments .. EXTERNAL hamat,sysbc C .. C .. Local Scalars .. DOUBLE PRECISION dl,du,e1,e2,shift,toloc,twopi INTEGER iflag,ind,knt LOGICAL duas C .. C .. Local Arrays .. DOUBLE PRECISION c1(17) C .. C .. External Functions .. DOUBLE PRECISION dlhsys,x01aaf EXTERNAL dlhsys,x01aaf C .. C .. External Subroutines .. EXTERNAL c05azf C .. C .. Intrinsic Functions .. INTRINSIC abs,max,min,dble C .. twopi = 2.0D0*x01aaf(0.0D0) shift = dble(k) + 5.0D-1 el = elam - eps eu = elam + eps duas = .false. knt = 0 ifail = 0 10 dl = dlhsys(xmesh,nmesh,ic,el,hamat,sysbc,work,lwk,n,ifail) - & shift IF (ifail.NE.0) RETURN C WRITE (6,FMT=*) 'elam,dmiss:',el,dl IF (dl.GT.0.0D0) THEN eu = el du = dl duas = .true. eps = 2.0D0*eps el = el - eps knt = knt + 1 IF (knt.GT.kntmax) THEN ifail = 7 RETURN END IF GO TO 10 END IF IF (duas) GO TO 30 knt = 0 du = dlhsys(xmesh,nmesh,ic,eu,hamat,sysbc,work,lwk,n,ifail) - & shift IF (ifail.NE.0) RETURN C WRITE (6,FMT=*) 'elam,dmiss:',eu,du 20 IF (dl*du.GT.0.0D0) THEN C We have not yet bracketed an eigenvalue knt = knt + 1 IF (knt.GT.kntmax) THEN ifail = 7 RETURN END IF eps = 2.0D0*eps IF (du.LT.0.0D0) THEN el = eu dl = du eu = eu + eps du = dlhsys(xmesh,nmesh,ic,eu,hamat,sysbc,work,lwk,n, & ifail) - shift IF (ifail.NE.0) RETURN C WRITE (6,FMT=*) 'elam,dmiss:',eu,du ELSE C I think that this next chunk of code can never be executed unless C by error in the function values! eu = el du = dl el = el - eps dl = dlhsys(xmesh,nmesh,ic,el,hamat,sysbc,work,lwk,n, & ifail) - shift IF (ifail.NE.0) RETURN C WRITE (6,FMT=*) 'elam,dmiss:',el,dl END IF GO TO 20 END IF C If we are here then it means that we have bracketed an eigenvalue and C we can now think about locating it accurately with c05azf. 30 knt = 0 ind = -1 toloc = tol c1(1) = du iflag = 1 e1 = el e2 = eu 40 CALL c05azf(e1,e2,dl,toloc,1,c1,ind,iflag) IF (ind.EQ.0) GO TO 50 knt = knt + 1 IF (knt.GT.kntmax) THEN ifail = 7 RETURN END IF dl = dlhsys(xmesh,nmesh,ic,e1,hamat,sysbc,work,lwk,n,ifail) - & shift IF (ifail.NE.0) RETURN C WRITE (6,FMT=*) 'elam,dmiss:',e1,dl IF (dl.GT.0.0D0) eu = min(eu,e1) IF (dl.LT.0.0D0) el = max(el,e1) IF (abs(eu-el).LT.toloc*max(1.0D0,abs(el))) GO TO 50 C Modify to a continuous miss distance function IF (abs(dl).LT.0.6D0) THEN IF (dl.LT.0.0D0) THEN dl = work(n) - twopi ELSE dl = work(1) END IF C WRITE (6,FMT=*) 'dmiss (modified):',dl END IF GO TO 40 50 elam = el C WRITE (6,FMT=*) 'el,eu,eu-el:',el,eu,eu - el RETURN END C ---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION dlhsys(xmesh,nmesh,ic,elam,hamat,sysbc, & work,lwk,n,ifail) C This routine computes the Atkinson miss-distance function for a C Hamiltonian system. C .. Parameters .. DOUBLE PRECISION zero,one,two PARAMETER (zero=0.0D0,one=1.0D0,two=2.0D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION elam INTEGER ic,ifail,n,nmesh C .. C .. Array Arguments .. DOUBLE PRECISION work(n* (24*n+12)),xmesh(0:nmesh) INTEGER lwk(n) C .. C .. Subroutine Arguments .. EXTERNAL hamat,sysbc C .. C .. Local Scalars .. DOUBLE PRECISION alfsum,argc,argl,argr,twopi INTEGER i,ir,iu,iuo,iv,ivo,iwork,j LOGICAL altnor,ising,matv C .. C .. External Functions .. DOUBLE PRECISION x01aaf EXTERNAL x01aaf C .. C .. External Subroutines .. EXTERNAL ethtex,f06qff,m01caf,spthta,sumvec,thstp4 C .. ifail = 0 twopi = two*x01aaf(one) iu = 1 iv = n*n + iu iuo = n*n + iv ivo = n*n + iuo ir = n*n + ivo iwork = ir + n*n C Step 1: set up the boundary condition at the left-hand end. CALL sysbc(0,ising,xmesh(0),work(iu),work(iv),n,n,elam) C Step 2: Compute the eigenvalues of the Atkinson theta-matrix associated with C the initial conditions and store in work(1,..,n) altnor = .false. CALL spthta(work(iu),work(iv),n,n,work(iwork),work(ir), & work(iwork+n),lwk,altnor,ifail) IF (ifail.NE.0) THEN ifail = 6 RETURN END IF CALL sumvec(work(iwork),n,argl) alfsum = argl C Step 3: Integrate the system, following the eigenvalues as we go. DO 10 i = 1,ic C Integrate from xmesh(i-1) to xmesh(i): CALL thstp4(xmesh(i-1),xmesh(i),work(iu),work(iv),n,n,elam, & argl,alfsum,work(iwork),lwk,hamat,ifail) IF (ifail.NE.0) RETURN 10 CONTINUE C Completed left-hand leg. Now turn to right hand leg. C Step 4: set up the boundary condition at the right-hand end. CALL sysbc(1,ising,xmesh(nmesh),work(iuo),work(ivo),n,n,elam) C Step 5: Compute the eigenvalues of the Atkinson theta-matrix associated with C the initial conditions. altnor = .true. CALL spthta(work(iuo),work(ivo),n,n,work(iwork),work(ir), & work(iwork+n),lwk,altnor,ifail) IF (ifail.NE.0) THEN ifail = 6 RETURN END IF CALL sumvec(work(iwork),n,argr) alfsum = 0.0D0 DO 20 i = iwork,iwork + n - 1 IF (work(i).NE.twopi) alfsum = alfsum + work(i) 20 CONTINUE C Step 6: Integrate the system, following the eigenvalues as we go. DO 30 i = nmesh,ic + 1,-1 C Integrate from xmesh(i) to xmesh(i-1): CALL thstp4(xmesh(i),xmesh(i-1),work(iuo),work(ivo),n,n,elam, & argr,alfsum,work(iwork),lwk,hamat,ifail) IF (ifail.NE.0) RETURN 30 CONTINUE C We now have the eigenvalues for both left and right-hand theta matrices. C Use these to compute the miss-distance function. matv = .true. CALL ethtex(work(iu),work(iv),work(iuo),work(ivo),n,n,work(iwork), & work(ir),work(iwork+n),n* (2*n+3),lwk,n,matv,ifail) IF (ifail.NE.0) RETURN C Order the eigenvalues: CALL m01caf(work(iwork),1,n,'A',ifail) C WRITE (6,FMT=*) 'Eigenvalues of theta_{miss}:', C & (work(j),j=iwork,iwork+n-1) CALL sumvec(work(iwork),n,argc) C We may now compute the required miss-distance function. dlhsys = (argl-argr-argc)/twopi + n C Also for convenience of the calling program let us copy the phase angles C into the first n elements of the array work. CALL f06qff('G',n,1,work(iwork),n,work(1),n) RETURN END C ----------------------------------------------------------------------- SUBROUTINE thstp4(xo,xend,u,v,ndim,n,elam,argdet,alfsum,work,lwk, & hamat,ifail) C .. Scalar Arguments .. DOUBLE PRECISION alfsum,argdet,elam,xend,xo INTEGER ifail,n,ndim C .. C .. Array Arguments .. DOUBLE PRECISION u(ndim,n),v(ndim,n),work(n* (18*n+12+ndim)) INTEGER lwk(n) C .. C .. Local Scalars .. INTEGER ia,ia0,ia1,ia2,ib,ibet,ibs,ic,id,ievals,im,ir,iscl,ist1, & ist2,iwork,n2 C .. C .. External Subroutines .. EXTERNAL thstp3 C .. C .. Subroutine Arguments .. EXTERNAL hamat C .. n2 = n*n ia = 1 ib = ia + n2 ibs = ib + n2 ic = ibs + n2 ia0 = ic + n2 ia1 = ia0 + n2 ia2 = ia1 + n2 ist1 = ia2 + n2 ist2 = ist1 + n2 im = ist2 + n2 ir = im + n2 ibet = ir + ndim*n id = ibet + n iscl = id + n ievals = iscl + n iwork = ievals + n ifail = 0 CALL thstp3(xo,xend,u,v,ndim,n,elam,argdet,alfsum,work(ia), & work(ib),work(ibs),work(ic),work(ia0),work(ia1), & work(ia2),work(ist1),work(ist2),work(im),work(ir), & work(ibet),work(id),work(iscl),work(ievals), & work(iwork),lwk,hamat,ifail) RETURN END SUBROUTINE thstp3(xo,xend,u,v,ndim,n,elam,argdet,alfsum,a,b,bsqrt, & c,a0,a1,a2,store1,store2,m,r,beta,d,sclrow, & evals,work,lwk,hamat,ifail) C The purpose of this subroutine is to integrate a Hamiltonian C system approximately over an interval [xo,xend] on which the C coefficients are constant. C .. Parameters .. DOUBLE PRECISION zero,qtr,half,one,two PARAMETER (zero=0.0D0,qtr=2.5D-1,half=0.50D0,one=1.0D0,two=2.0D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION alfsum,argdet,elam,xend,xo INTEGER ifail,n,ndim C .. C .. Array Arguments .. DOUBLE PRECISION a(n,n),a0(n,n),a1(n,n),a2(n,n),b(n,n),beta(n), & bsqrt(n,n),c(n,n),d(n),evals(n),m(n,n),r(ndim,n), & sclrow(n),store1(n,n),store2(n,n),u(ndim,n), & v(ndim,n),work(8*n*(n+1)) INTEGER lwk(n) C .. C .. Local Scalars .. DOUBLE PRECISION alfnew,alpha,arg0,argbet,betah,csdk,diff,difmax, & dk,dkssdk,dyold,fac,h,ht,psdk,scmax,scmin,sdk, & sgn,ssdk,t,twopi,yold INTEGER i,imax,j,k,kpjm1n LOGICAL altnor C .. C .. External Functions .. DOUBLE PRECISION armax,armin,phi,psi,scl,spexp,x01aaf INTEGER x02bkf EXTERNAL armax,armin,phi,psi,scl,spexp,x01aaf,x02bkf C .. C .. External Subroutines .. EXTERNAL acomp,f02abf,f04aaf,f06qff,f06yaf,m01caf,matexp,spthta, & sqrmat,sumvec C .. C .. Intrinsic Functions .. INTRINSIC abs,atan2,cos,max,min,sign,sin,sqrt,float C .. C .. Subroutine Arguments .. EXTERNAL hamat C .. ifail = 0 C Evaluate the coefficients and compute any square roots of C matrices which may be required. C CALL hamat(half* (xo+xend),elam,a,b,c,store1,n) CALL hamat(half*(xo+xend),elam,c,store1,a,b,n) DO 7020 j = 1,n DO 7030 i = 1,n c(i,j) = -c(i,j) 7030 CONTINUE 7020 CONTINUE C Get the matrix sqrt(b): CALL sqrmat(b,n,bsqrt,n,n,work,ifail) IF (ifail.NE.0) THEN IF (ifail.EQ.1) THEN ifail = 20 ELSE ifail = 21 END IF END IF C Now form the matrix A0 = inv(bsqrt).a.bsqrt using an appropriate NAG C routine. Note that bsqrt has to be copied into a storage array as C this is overwritten by f04aaf CALL f06qff('G',n,n,bsqrt,n,store1,n) ifail = 1 CALL f04aaf(store1,n,a,n,n,n,store2,n,work,ifail) IF (ifail.NE.0) THEN C WRITE (7,*) 'Matrix bsqrt (1):' C CALL monit(bsqrt,n,n,7) ifail = 21 RETURN END IF C The matrix in store2 is inv(bsqrt).a, now post-multiply by bsqrt CALL f06yaf('N','N',n,n,n,one,store2,n,bsqrt,n,zero,store1,n) C The matrix in store1 is now A0 CALL f06qff('G',n,n,store1,n,a0,n) C Now form the matrix A1 = 0.5*(A0+transpose(A0)) and the matrix C A2 = 0.5*(A0-transpose(A0)). C (This is a pretty crude piece of code but it will do until such time as C I can think of something better DO 20 j = 1,n DO 10 i = j,n a1(i,j) = half* (a0(i,j)+a0(j,i)) a1(j,i) = a1(i,j) a2(i,j) = half* (a0(i,j)-a0(j,i)) a2(j,i) = -a2(i,j) 10 CONTINUE 20 CONTINUE C Form the matrix C0 = bsqrt.c.bsqrt and overwrite it on c: CALL f06yaf('N','N',n,n,n,one,c,n,bsqrt,n,zero,store1,n) CALL f06yaf('N','N',n,n,n,one,bsqrt,n,store1,n,zero,c,n) C C Now we must form the matrix M(x), which is given by C M(x) = exp(A2(x-xmid)) where xmid = (xo+xend)/2. We will compute this C matrix for the case where x = xend; the case x=xo will be the inverse of C the case x=xend. The routine matexp computes the exponential for us. C fac = half* (xend-xo) DO 40 j = 1,n DO 30 i = 1,n store1(i,j) = a2(i,j)*fac store2(i,j) = zero 30 CONTINUE 40 CONTINUE CALL matexp(store1,store2,n,n,2,m,work(1),n,work(n*n+1),ifail) IF (ifail.NE.0) THEN IF (ifail.EQ.1) THEN ifail = 17 ELSE ifail = 18 END IF RETURN END IF C WRITE (6,*) 'M matrix is:' C DO 17 i = 1,n C WRITE (6,*) (real(m(i,j)),j=1,n) C17 CONTINUE C The matrix m has now been computed. Now carry out the initial transformation C from the matrices u and v to the new matrices u2 and v2 given by C u2 = inv(M(xo))*inv(bsqrt)*u, v2 = inv(M(xo))*(bsqrt*v+a1*inv(bsqrt)*u) C Note that inv(M(xo)) = M(xend) = m, which saves us some work. CALL f06qff('G',n,n,bsqrt,n,store1,n) C Use the array A0 to contain inv(bsqrt)*u, since we don't need its current C contents any more. ifail = 1 CALL f04aaf(store1,n,u,ndim,n,n,a0,n,work,ifail) IF (ifail.NE.0) THEN C WRITE (7,*) 'Matrix bsqrt (2):' C CALL monit(bsqrt,n,n,7) ifail = 21 RETURN END IF C Now form u2 (overwrite it on u): CALL f06yaf('N','N',n,n,n,one,m,n,a0,n,zero,u,ndim) C Now form v2 (overwrite it on v): CALL f06yaf('N','N',n,n,n,one,a1,n,a0,n,zero,store1,n) CALL f06yaf('N','N',n,n,n,one,bsqrt,n,v,ndim,one,store1,n) CALL f06yaf('N','N',n,n,n,one,m,n,store1,n,zero,v,ndim) C Now find the phase angles of the matrix (v2+iu2)*inv(v2-iu2) using the C routine spthta. ifail = 0 CALL spthta(u,v,n,ndim,work(1),r,work(n+1),lwk,.false.,ifail) IF (ifail.NE.0) THEN ifail = 16 RETURN END IF CALL sumvec(work(1),n,argbet) C In perparation for the integration step we compute the matrix C which Reid claims is given by C C2 = C0 + A1*A2 -A2*A1 C but is actually given by C C2 = C0 + A1*A1 (which we overwrite on C) C CALL f06qff('G',n,n,a1,n,a2,n) CALL f06yaf('N','N',n,n,n,one,a1,n,a2,n,one,c,n) C CALL f06yaf('N','N',n,n,n,-one,a2,n,a1,n,one,c,n) C Now we are ready to go. The next chunk of code we lift straight from the C old coefficient approximation routine. twopi = two*x01aaf(one) h = xend - xo ht = xend/ (one+abs(xend)) - xo/ (one+abs(xo)) sgn = one C Use a NAG matrix eigenvalue routine to reduce c to diagonal C form. Note: it might be a good idea to do some balancing first. C WRITE (7,*) 'xo,xend:',xo,xend C WRITE (7,*) 'Q matrix after Reid''s transformation:' C CALL monit(c,n,n,7) CALL f02abf(c,n,n,d,r,ndim,work,ifail) IF (ifail.NE.0) THEN ifail = 9 RETURN END IF C WRITE (7,*) 'Eigenvalues of Q matrix:' C DO 17 jj = 1,n C WRITE (7,*) jj,d(jj) C17 CONTINUE C The columns of the array r are the eigenvectors of c. We may put these C into c for safe storage: having diagonalised c we do not need it any more. CALL f06qff('G',n,n,r,ndim,c,n) C Transform dependent variables: CALL f06qff('G',n,n,u,ndim,work(1),n) CALL f06yaf('T','N',n,n,n,one,r,ndim,work,n,zero,u,ndim) CALL f06qff('G',n,n,v,ndim,work(1),n) CALL f06yaf('T','N',n,n,n,one,r,ndim,work,n,zero,v,ndim) C Finished transformation. C Integrate from xo to xend: DO 60 k = 1,n dk = -d(k) sdk = sqrt(abs(dk)) psdk = sdk t = dk* (h**2) IF (t.GE.0.1) THEN csdk = cos(sdk*h) ssdk = sin(sdk*h)/sdk dkssdk = -dk*ssdk sclrow(k) = zero evals(k) = two*scl(sdk*h,one/sdk) ELSE IF (abs(t).LT.0.1) THEN csdk = psi(-t) ssdk = h*phi(-t) dkssdk = -dk*ssdk sclrow(k) = zero evals(k) = two*atan2(ssdk,csdk) ELSE C IF (t.LE.-0.1) THEN sclrow(k) = abs(sdk*h) arg0 = spexp(-two*sclrow(k)) csdk = half* (one+arg0) ssdk = half*sign(one,h)* (one-arg0)/sdk dkssdk = -dk*ssdk evals(k) = two*atan2(ssdk,csdk) END IF kpjm1n = k DO 50 j = 1,n yold = u(k,j) dyold = v(k,j) u(k,j) = csdk*yold + ssdk*dyold v(k,j) = dkssdk*yold + csdk*dyold C Additional matrix needed for technical reasons: work(kpjm1n) = csdk*dyold - ssdk*yold kpjm1n = kpjm1n + n 50 CONTINUE 60 CONTINUE C Finished integration. yold = zero DO 80 j = 1,n DO 70 i = 1,n yold = max(yold,abs(u(i,j)),abs(v(i,j))) 70 CONTINUE 80 CONTINUE scmax = armax(sclrow,n) scmin = armin(sclrow,n) fac = qtr*float(min(127,-x02bkf(one))) C IF (scmax-scmin.GT.fac*abs(ht).OR.yold.GT.toobig.OR. C & .NOT.cheap) THEN C Get argdet((C+iS)*inv(C-iS)): CALL sumvec(evals,n,arg0) C Get eigenvalues of phi (stored in evals): altnor = .false. CALL f06qff('G',n,n,u,ndim,work(n**2+1),n) ifail = 1 CALL spthta(work(n**2+1),work(1),n,n,evals,r,work(2* (n**2)+1), & lwk,altnor,ifail) IF (ifail.NE.0) THEN ifail = 16 RETURN END IF C Form the sum of these phase angles: CALL sumvec(evals,n,betah) C Now get the eigenvalues of theta2 correct including multiples of 2pi. C This is a bit of a pain because we don't actually need these eigenvalues C except to find a well-conditioned normalising matrix. CALL spthta(u,v,n,ndim,d,r,work,lwk,altnor,ifail) IF (ifail.NE.0) THEN ifail = 16 RETURN END IF C Now prepare for possible renormalisation of y and dy. First step is to C select alpha such that cos(alpha/2)*y-sin(alpha/2)*dy is nowhere near C singular. C First order the eigenvalues using m01caf: CALL m01caf(d,1,n,'A',ifail) difmax = twopi + d(1) - d(n) imax = n DO 90 i = 1,n - 1 diff = d(i+1) - d(i) IF (diff.GT.difmax) THEN imax = i difmax = diff END IF 90 CONTINUE C Select phase-shift alpha for next step: IF (imax.LT.n) THEN alpha = half* (d(imax+1)+d(imax)) ELSE alpha = half* (twopi+d(1)+d(n)) END IF C work has to be of dimension at least n*(2*n+1) here altnor = .true. C Re-assign R to be the orthogonal diagonalising matrix: CALL f06qff('G',n,n,c,n,r,ndim) CALL acomp(u,v,ndim,r,ndim,sclrow,alpha,work,n,altnor,sgn,ifail) IF (ifail.NE.0) THEN ifail = 19 RETURN END IF C We have now post-multiplied u2 and v2 by a suitable normalising C matrix. We can transform them back to u and v now. The formulae are C u = bsqrt*M(xend)*u2 and v = inv(bsqrt)*(M(xend)*v2-a1*M(xend)*u2) CALL f06yaf('N','N',n,n,n,one,m,n,u,ndim,zero,store1,n) CALL f06yaf('N','N',n,n,n,one,bsqrt,n,store1,n,zero,u,ndim) CALL f06yaf('N','N',n,n,n,-one,a1,n,store1,ndim,zero,store2,n) CALL f06yaf('N','N',n,n,n,one,m,n,v,ndim,one,store2,n) C Now we premultiply the contents of store2 by inv(bsqrt): CALL f06qff('G',n,n,bsqrt,n,store1,n) ifail = 1 CALL f04aaf(store1,n,store2,n,n,n,v,ndim,work,ifail) IF (ifail.NE.0) THEN C WRITE (7,*) 'Matrix bsqrt (3):' C CALL monit(bsqrt,n,n,7) ifail = 21 RETURN END IF C Now get the eigenvalues of the theta matrix: altnor = .false. ifail = 0 CALL spthta(u,v,n,ndim,d,r,work,lwk,altnor,ifail) IF (ifail.NE.0) THEN ifail = 16 RETURN END IF CALL sumvec(d,n,alfnew) argdet = argdet - alfsum + alfnew + arg0 + argbet - betah alfsum = alfnew RETURN END C ---------------------------------------------------------------------- SUBROUTINE matexp(ar,ai,ia,n,iswtch,expr,expi,iexp,work,ifail) C This subroutine takes in a matrix ar+i*ai, which is either Hermitian or C skew-Hermitian, diagonalises it, and hence computes its exponential. C .. Scalar Arguments .. INTEGER ia,iexp,ifail,iswtch,n C .. C .. Array Arguments .. DOUBLE PRECISION ai(ia,n),ar(ia,n),expi(iexp,n),expr(iexp,n), & work(4*n* (n+1)) C .. C .. Local Scalars .. DOUBLE PRECISION c,s INTEGER i,index,ir,ivi,ivis,ivr,ivrs,iwk1,iwk2,iwk3,j C .. C .. External Subroutines .. EXTERNAL f02axf,f06yaf C .. C .. Intrinsic Functions .. INTRINSIC abs,cos,exp,sin C .. IF (iswtch.EQ.1) THEN C The matrix is Hermitian and has real eigenvalues ir = 1 ivr = n + 1 ivi = ivr + n*n iwk1 = ivi + n*n iwk2 = iwk1 + n iwk3 = iwk2 + n ivrs = iwk3 + n ivis = ivrs + n*n ifail = -1 CALL f02axf(ar,ia,ai,ia,n,work(ir),work(ivr),n,work(ivi),n, & work(iwk1),work(iwk2),work(iwk3),ifail) IF (ifail.NE.0) THEN ifail = 1 RETURN END IF C Now get the exponential of the eigenvalues: index = ir DO 10 j = 1,n work(index) = exp(work(index)) index = index + 1 10 CONTINUE C Now the matrix exponential is given by C C expr + i*expi = (rr+i*ri)*exp(D)*(transpose(rr)-i*transpose(ri)) C index = 0 DO 30 j = 1,n DO 20 i = 1,n work(ivrs+index) = work(ivr+index)*work(ir+j-1) work(ivis+index) = work(ivi+index)*work(ir+j-1) index = index + 1 20 CONTINUE 30 CONTINUE CALL f06yaf('N','T',n,n,n,1.0D0,work(ivr),n,work(ivrs),n, & 0.0D0,expr,iexp) CALL f06yaf('N','T',n,n,n,1.0D0,work(ivi),n,work(ivis),n, & 1.0D0,expr,iexp) CALL f06yaf('N','T',n,n,n,1.0D0,work(ivi),n,work(ivrs),n, & 0.0D0,expi,iexp) CALL f06yaf('N','T',n,n,n,-1.0D0,work(ivr),n,work(ivis),n, & 1.0D0,expi,iexp) C expr + i*expi = (rr+i*ri)*exp(D)*(transpose(rr)-i*transpose(ri)) ELSE C In this case the matrix A is anti-Hermitian. This means that -iA is C Hermitian. Note that -iA = ai - i*ar. We find the eigenvalues of -iA, C which are all real, and translate them back into the eigenvalues of A C which follow from A = i(-iA). C Negate all the elements of ar: DO 50 j = 1,n DO 40 i = 1,n ar(i,j) = -ar(i,j) 40 CONTINUE IF (abs(ar(j,j)).GT.1.0D-6) THEN ifail = 2 RETURN END IF ar(j,j) = 0.0D0 50 CONTINUE C Get eigenvalues: ir = 1 ivr = n + 1 ivi = ivr + n*n iwk1 = ivi + n*n iwk2 = iwk1 + n iwk3 = iwk2 + n ivrs = iwk3 + n ivis = ivrs + n*n ifail = -1 CALL f02axf(ai,ia,ar,ia,n,work(ir),work(ivr),n,work(ivi),n, & work(iwk1),work(iwk2),work(iwk3),ifail) IF (ifail.NE.0) THEN ifail = 1 RETURN END IF C Now we can return the matrix A to its original form: DO 70 j = 1,n DO 60 i = 1,n ar(i,j) = -ar(i,j) 60 CONTINUE 70 CONTINUE C This time we need both sines and cosines of the eigenvalues in order to form C the matrix exponential. index = 0 DO 90 j = 1,n DO 80 i = 1,n c = cos(work(ir+j-1)) s = sin(work(ir+j-1)) work(ivrs+index) = work(ivr+index)*c + & work(ivi+index)*s work(ivis+index) = work(ivr+index)*s - & work(ivi+index)*c index = index + 1 80 CONTINUE 90 CONTINUE CALL f06yaf('N','T',n,n,n,1.0D0,work(ivr),n,work(ivrs),n, & 0.0D0,expr,iexp) CALL f06yaf('N','T',n,n,n,-1.0D0,work(ivi),n,work(ivis),n, & 1.0D0,expr,iexp) CALL f06yaf('N','T',n,n,n,1.0D0,work(ivr),n,work(ivis),n, & 0.0D0,expi,iexp) CALL f06yaf('N','T',n,n,n,1.0D0,work(ivi),n,work(ivrs),n, & 1.0D0,expi,iexp) END IF ifail = 0 RETURN END C ---------------------------------------------------------------------- SUBROUTINE sqrmat(a,ia,as,ias,n,work,ifail) C The purpose of this routine is to compute the square root of positive C definite Hermitian matrix A. C .. Scalar Arguments .. INTEGER ia,ias,ifail,n C .. C .. Array Arguments .. DOUBLE PRECISION a(ia,n),as(ias,n),work(n* (2*n+1)) C .. C .. Local Scalars .. INTEGER i C .. C .. External Subroutines .. EXTERNAL f02abf,rsmble C .. C .. Intrinsic Functions .. INTRINSIC sqrt C .. ifail = 0 C First we obtain the Schur decomposition of the matrix: CALL f02abf(a,ia,n,work(1),work(2*n+1),n,work(n+1),ifail) IF (ifail.NE.0) THEN ifail = 1 RETURN END IF C Now we take the square roots of the eigenvalues DO 10 i = 1,n IF (work(i).LT.0.0D0) THEN ifail = 2 RETURN END IF work(i) = sqrt(work(i)) 10 CONTINUE C Now reassemble the matrix sqrt(A): CALL rsmble(as,ias,n,work(1),work(2*n+1),n) RETURN END C --------------------------------------------------------------------- SUBROUTINE rsmble(as,is,n,d,r,ir) C .. Scalar Arguments .. INTEGER ir,is,n C .. C .. Array Arguments .. DOUBLE PRECISION as(is,n),d(n),r(ir,n) C .. C .. Local Scalars .. DOUBLE PRECISION fac1,fac2 INTEGER i,j,k C .. DO 20 j = 1,n DO 10 i = 1,n as(i,j) = 0.0D0 10 CONTINUE 20 CONTINUE DO 50 k = 1,n fac1 = d(k) DO 40 j = 1,n fac2 = r(j,k) DO 30 i = 1,n as(i,j) = as(i,j) + r(i,k)*fac1*fac2 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END C ----------------------------------------------------------------------- SUBROUTINE monit(a,ia,n,iprn) C .. Scalar Arguments .. INTEGER ia,iprn,n C .. C .. Array Arguments .. DOUBLE PRECISION a(ia,n) C .. C .. Local Scalars .. INTEGER i,j C .. C .. Intrinsic Functions .. INTRINSIC real C .. DO 10 i = 1,n WRITE (iprn,FMT=*) (real(a(i,j)),j=1,n) 10 CONTINUE RETURN END C --------------------------------------------------------------------- C --------------------------------------------------------------------- C C Additional subroutines required by SL12F. C (C) Copyright Marco Marletta, 1991. All rights reserved. C SUBROUTINE spthta(y,dy,n,ndim,omega,r,work,lwk,lnor,ifail) C Transfer y and dy values in order to avoid overwriting. C .. Parameters .. DOUBLE PRECISION zero,one PARAMETER (zero=0.d0,one=1.d0) C .. C .. Scalar Arguments .. INTEGER ifail,n,ndim LOGICAL lnor C .. C .. Array Arguments .. DOUBLE PRECISION dy(ndim,n),omega(n),r(ndim,n),work(n* (2*n+3)), & y(ndim,n) INTEGER lwk(n) C .. C .. Local Scalars .. DOUBLE PRECISION alpha,beta,fac,eps,twopi INTEGER i,knt,n2 LOGICAL matv C .. C .. External Subroutines .. EXTERNAL f06qff,f02bjf C .. C .. Intrinsic Functions .. INTRINSIC atan2,sign C .. C .. External Functions .. DOUBLE PRECISION x01aaf,x02ajf EXTERNAL x01aaf,x02ajf C .. LOGICAL prin COMMON/printb/prin n2 = n**2 CALL f06qff('G',n,n,y,ndim,work(1),n) CALL f06qff('G',n,n,dy,ndim,work(n*n+1),n) ifail = 1 matv = .false. eps = x02ajf(zero) twopi = 2.0D0*x01aaf(zero) knt = 0 100 CALL f02bjf(n,work(1),n,work(n*n+1),n,eps,work(2*n2+1), & work(2*n2+n+1),work(2*(n2+n)+1),matv,r,ndim,lwk, & ifail) IF (ifail.NE.0) THEN knt = knt + 1 ifail = 1 IF (knt.LT.4) GO TO 100 ifail = 8 RETURN END IF DO 30 i = 1,n alpha = work(2*n2+i) beta = work(2*(n2+n)+i) IF ((abs(alpha)+abs(beta)).EQ.zero) THEN ifail = 7 RETURN ENDIF IF (alpha.EQ.zero) THEN omega(i) = 0.0D0 IF (lnor) omega(i) = twopi GO TO 30 END IF fac = sign(one,alpha) alpha = fac*alpha beta = fac*beta omega(i) = 2.d0*atan2(alpha,beta) 30 CONTINUE RETURN END C --------------------------------------------------------------------- SUBROUTINE ethtex(ul,vl,ur,vr,nsys,ndim,omega,vmat,wk,iwk,lwk, & ilwk,matv,ifail) C This routine acts as the acceptable face of SUBROUTINE sptht0. C Its purpose is described in the in-code comment for sptht0. C iwk must be at least 2*(nsys**2) + 3*nsys. C .. Scalar Arguments .. INTEGER ifail,ilwk,iwk,ndim,nsys LOGICAL matv C .. C .. Array Arguments .. DOUBLE PRECISION omega(nsys),vmat(1:ndim,1:nsys),ul(ndim,ndim), & ur(ndim,ndim),vl(ndim,ndim),vr(ndim,ndim),wk(iwk) INTEGER lwk(ilwk) C .. C .. Local Scalars .. INTEGER i,i1,i2,i3,i4,i5,i6 C .. C .. External Subroutines .. EXTERNAL sptht0 C .. i1 = 1 i = nsys**2 i2 = i + 1 i3 = i2 + i i4 = i3 + nsys i5 = i4 + nsys CALL sptht0(ul,vl,ur,vr,nsys,ndim,omega,wk(i1),wk(i2),vmat, & wk(i3),wk(i4),wk(i5),lwk,matv,ifail) RETURN END C --------------------------------------------------------------------- SUBROUTINE sptht0(ul,vl,ur,vr,nsys,ndim,omega,a,b,v,alfar,alfai, & beta,iter,matv,ifail) DOUBLE PRECISION zero,one PARAMETER (zero=0.d0,one=1.d0) C Subroutine to find the eigenvalues of the unitary matrix thetar^{*}.thetal C where thetal and thetar are unitary matrices given by the Cayley transform C thetal(x) = (vl(x)+iul(x))(vl(x)-iul(x))^{-1} C thetar(x) = (vr(x)+iur(x))(vr(x)-iur(x))^{-1} C C The eigenvalues are all of the form cos(w)+isin(w), where w satisfies C C det( cos(w/2)(vr^{*}ul-ur^{*}vl) - sin(w/2)(vr^{*}vl+ur^{*}ul) ) = 0 C We use the QZ algorithm to find cos(w/2) and sin(w/2), and thus recover w. C C .. Scalar Arguments .. INTEGER ifail,ndim,nsys LOGICAL matv C .. C .. Array Arguments .. DOUBLE PRECISION a(nsys,nsys),alfai(nsys),alfar(nsys), & b(nsys,nsys),beta(nsys),omega(nsys), & ul(ndim,nsys),ur(ndim,nsys),v(ndim,nsys), & vl(ndim,nsys),vr(ndim,nsys) INTEGER iter(nsys) C .. C .. Local Scalars .. DOUBLE PRECISION eps1,fac,sum,ulmax,urmax,vlmax,vrmax INTEGER i C .. C .. External Functions .. DOUBLE PRECISION x02aaf,x01aaf EXTERNAL x02aaf,x01aaf C .. C .. External Subroutines .. EXTERNAL f06yaf,f06qff,f06qhf,f02bjf,m01caf C .. C .. Intrinsic Functions .. INTRINSIC abs,atan2 C .. eps1 = 0.d0 ifail = 0 C First set A and B to zero: CALL f06qhf('G',nsys,nsys,zero,zero,a,nsys) CALL f06qhf('G',nsys,nsys,zero,zero,b,nsys) IF (matv) THEN C Form the matrices A and B given by C A = vr^{*}ul - ur^{*}vl C B = vr^{*}vl + ur^{*}ul C (where * denotes transpose) fac = 1.d0 CALL f06yaf('T','N',nsys,nsys,nsys,fac,vr,ndim,ul,ndim,fac, & a,nsys) CALL f06yaf('T','N',nsys,nsys,nsys,-fac,ur,ndim,vl,ndim,fac, & a,nsys) CALL f06yaf('T','N',nsys,nsys,nsys,fac,vr,ndim,vl,ndim,fac, & b,nsys) CALL f06yaf('T','N',nsys,nsys,nsys,fac,ur,ndim,ul,ndim,fac, & b,nsys) ELSE C Set a = ul and b = vl CALL f06qff('G',nsys,nsys,ul,ndim,a,nsys) CALL f06qff('G',nsys,nsys,vl,ndim,b,nsys) ENDIF C Step 2: find the simultaneous eigenvalues of A and B. ifail = 1 CALL f02bjf(nsys,a,nsys,b,nsys,eps1,alfar,alfai,beta,matv,v,ndim, & iter,ifail) IF (ifail.NE.0) THEN IFAIL = 8 RETURN END IF sum = 0.d0 DO 13 i = 1,nsys fac = sign(one,alfar(i)) omega(i) = 2.d0*atan2(fac*alfar(i),fac*beta(i)) sum = sum + abs(alfai(i))/nsys 13 CONTINUE IF (sum.GT.x02aaf(0.d0)) THEN WRITE (6,FMT=*) 'Warning: imaginary eigenvalues detected' write(6,*) 'size of sum:',sum END IF RETURN END C --------------------------------------------------------------------- SUBROUTINE acomp(z,dz,iz,r,ir,sclrow,alpha,work,n,back,sgn,ifail) C .. Scalar Arguments .. DOUBLE PRECISION alpha,sgn INTEGER ifail,ir,iz,n,ndim LOGICAL back C .. C .. Array Arguments .. DOUBLE PRECISION dz(iz,n),r(ir,n),sclrow(n),work(8*n*n+n),z(iz,n) C .. C .. External Subroutines .. EXTERNAL compa C .. CALL compa(z,dz,iz,r,ir,sclrow,alpha,work(1),work(n*n+1), & work(3*n*n+1),work(5*n*n+1),n,back,sgn,ifail) RETURN END SUBROUTINE compa(z,dz,iz,r,ir,sclrow,alpha,astore,bstore,cstore, & work,n,back,sgn,ifail) C .. Parameters .. DOUBLE PRECISION zero,half,one PARAMETER (zero=0.d0,half=5.d-1,one=1.d0) C .. C .. Scalar Arguments .. DOUBLE PRECISION alpha,sgn INTEGER ifail,ir,iz,n,ndim LOGICAL back C .. C .. Array Arguments .. DOUBLE PRECISION astore(n,n),bstore(2*n*n), & cstore(2*n*n),dz(iz,n),r(ir,n),work(n*(3*n+1)), & z(iz,n),sclrow(n) C .. C .. Local Scalars .. DOUBLE PRECISION ca2,dummy,d1,eps,sa2,arg1,arg2,arg3 INTEGER i,id,j,k,p,ipjm1n,jpim1n,lwk(10) LOGICAL doy C .. C .. External Subroutines .. EXTERNAL f06qff,f06qhf,f03aff,f04ahf,getzdz,f04ayf,f01btf C .. C .. Intrinsic Functions .. INTRINSIC cos,sin C .. C .. External Functions .. DOUBLE PRECISION spexp,x01aaf,x02ajf EXTERNAL spexp,x01aaf,x02ajf C .. eps = x02ajf(one) sa2 = sin(half*alpha) ca2 = cos(half*alpha) C Transpose z and dz: CALL transp(z,iz,n) CALL transp(dz,iz,n) C Form matrix A: CALL lncomb(astore,n,z,iz,dz,iz,n,ca2,-sa2) CALL f06qff('G',n,n,astore,n,work(1),n) C Now multiply z and dz by the inverse of a: C Step 1: copy z and dz into a single unit of storage: doy = .true. IF (ca2**2.GE.half) doy = .false. IF (doy) THEN CALL f06qff('G',n,n,z,iz,bstore(1),n) ELSE CALL f06qff('G',n,n,dz,iz,bstore(1),n) END IF C Step 2: solve an appropriate linear system with multiple RHS: ifail = 1 C Get PLU factorisation: CALL f01btf(n,work(1),n,work(n*n+1),d1,ifail) IF (ifail.NE.0) THEN IFAIL = 10 RETURN END IF C CALL f03aff(n,eps,work(1),n,d1,id,work(n*n+1),ifail) sgn = sign(one,d1) index = 1 DO 17 i = 1,n sgn = sgn*sign(one,work(index)) index = index + n + 1 17 CONTINUE C Solve system: ifail = 1 CALL f06qff('G',n,n,bstore(1),n,cstore(1),n) CALL f04ayf(n,n,work(1),n,work(n*n+1),cstore(1),n,ifail) C CALL f04ahf(n,n,astore,n,work(1),n,work(n*n+1),bstore,n,eps, C & cstore,n,work(n*(n+1)+1),n,k,ifail) IF (ifail.NE.0) THEN IFAIL = 10 RETURN END IF C Recover the required solutions: IF (doy) THEN C Recover z: CALL f06qff('G',n,n,cstore(1),n,z,iz) C Store an identity matrix: CALL f06qhf('G',n,n,zero,one,cstore(1),n) C Get dz from z: CALL getzdz(z,dz,iz,cstore(1),n,n,ca2,sa2,1) ELSE C Recover dz: CALL f06qff('G',n,n,cstore(1),n,dz,iz) C Store an identity matrix: CALL f06qhf('G',n,n,zero,one,cstore(1),n) C Get z from dz: CALL getzdz(z,dz,iz,cstore(1),n,n,ca2,sa2,0) END IF C Take account to sclrow scaling: CALL recovr(astore,n,z,iz,sclrow,n,1) CALL recovr(bstore(1),n,dz,iz,sclrow,n,1) C Transpose z and dz again: CALL transp(z,iz,n) CALL transp(dz,iz,n) CALL f06qff('G',n,n,astore,n,z,iz) CALL f06qff('G',n,n,bstore(1),n,dz,iz) C We can set sclrow to zero now. CALL f06qhf('G',n,1,zero,zero,sclrow,n) IF (back) THEN C Apply the transformation z = r.z, dz = r.dz IF (doy) THEN C set astore = 0 CALL f06qhf('G',n,n,zero,zero,astore,n) CALL f06yaf('N','N',n,n,n,one,r,ir,z,iz,one,astore,n) CALL f06qff('G',n,n,astore,n,z,iz) iopt = 1 ELSE CALL f06qhf('G',n,n,zero,zero,astore,n) CALL f06yaf('N','N',n,n,n,one,r,ir,dz,iz,one,astore,n) CALL f06qff('G',n,n,astore,n,dz,iz) iopt = 0 END IF CALL getzdz(z,dz,iz,r,ir,n,ca2,sa2,iopt) END IF RETURN END C ------------------------------------------------------------------------- SUBROUTINE getzdz(z,dz,iz,r,ir,n,ca2,sa2,iopt) DOUBLE PRECISION zero,one PARAMETER (zero=0.d0,one=1.d0) INTEGER iz,n,iopt DOUBLE PRECISION z(iz,n),dz(iz,n),r(ir,n),ca2,sa2,arg1,arg2 EXTERNAL lncomb IF (iopt.EQ.1) THEN C Get dz from z arg1 = -one/sa2 arg2 = ca2/sa2 CALL f06qhf('G',n,n,zero,zero,dz,iz) CALL lncomb(dz,iz,r,ir,z,iz,n,arg1,arg2) ELSE C Get z from dz arg1 = one/ca2 arg2 = sa2/ca2 CALL f06qhf('G',n,n,zero,zero,z,iz) CALL lncomb(z,iz,r,ir,dz,iz,n,arg1,arg2) END IF RETURN END C --------------------------------------------------------------------- C Subroutine to compute the sum of the elements of a vector SUBROUTINE sumvec(vec,nvec,sum) C .. Scalar Arguments .. DOUBLE PRECISION sum INTEGER nvec C .. C .. Array Arguments .. DOUBLE PRECISION vec(nvec) C .. C .. Local Scalars .. INTEGER i C .. sum = 0.d0 DO 10 i = 1,nvec sum = sum + vec(i) 10 CONTINUE RETURN END C ---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION armax(v,n) C .. Scalar Arguments .. INTEGER n C .. C .. Array Arguments .. DOUBLE PRECISION v(1:n) C .. C .. Local Scalars .. INTEGER i C .. armax = v(1) DO 10 i = 2,n IF (armax.LT.v(i)) armax = v(i) 10 CONTINUE RETURN END C ---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION armin(v,n) C .. Scalar Arguments .. INTEGER n C .. C .. Array Arguments .. DOUBLE PRECISION v(1:n) C .. C .. Local Scalars .. INTEGER i C .. armin = v(1) DO 10 i = 2,n IF (armin.GT.v(i)) armin = v(i) 10 CONTINUE RETURN END C -------------------------------------------------------------------- SUBROUTINE lindep(v1,v2,n,tf) C This routine examines two vectors v1 and v2 for linear dependence C .. Parameters .. DOUBLE PRECISION zero,one PARAMETER (zero=0.d0,one=1.d0) C .. C .. Scalar Arguments .. INTEGER n LOGICAL tf C .. C .. Array Arguments .. DOUBLE PRECISION v1(n),v2(n) C .. C .. Local Scalars .. DOUBLE PRECISION v1norm,v1v2ip,v2norm INTEGER i C .. C .. External Functions .. DOUBLE PRECISION x02ajf EXTERNAL x02ajf C .. C .. Intrinsic Functions .. INTRINSIC abs,sqrt C .. v1norm = zero v2norm = zero v1v2ip = zero DO 10 i = 1,n v1norm = v1norm + v1(i)**2 v2norm = v2norm + v2(i)**2 v1v2ip = v1v2ip + v1(i)*v2(i) 10 CONTINUE IF (v1norm.LE.x02ajf(one) .OR. v2norm.LE.x02ajf(one)) THEN tf = .true. RETURN END IF IF (abs(abs(v1v2ip)/sqrt(v1norm*v2norm)-one).LT. & sqrt(x02ajf(one))) THEN tf = .true. ELSE tf = .false. END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE recovr(y,iy,a,ia,d,n,iopt) DOUBLE PRECISION zero PARAMETER (zero=0.d0) INTEGER iy,ia,n,iopt DOUBLE PRECISION y(iy,n),a(ia,n),d(n) DOUBLE PRECISION diff DOUBLE PRECISION spexp EXTERNAL spexp GOTO (100,200,300,400) iopt 100 DO 10 j = 1,n DO 20 i = 1,n diff = d(i)-d(j) IF (diff.LE.zero) THEN y(i,j) = spexp(diff)*a(j,i) y(j,i) = y(i,j) ELSE y(j,i) = spexp(-diff)*a(i,j) y(i,j) = y(j,i) END IF 20 CONTINUE 10 CONTINUE RETURN 200 DO 30 j = 1,n DO 40 i = 1,n y(i,j) = a(j,i) 40 CONTINUE 30 CONTINUE RETURN 300 DO 50 j = 1,n DO 60 i = 1,n diff = -d(i)-d(j) y(i,j) = spexp(diff)*a(j,i) 60 CONTINUE 50 CONTINUE RETURN 400 DO 70 j = 1,n diff = spexp(-d(j)) DO 80 i = 1,n y(i,j) = diff*a(j,i) 80 CONTINUE 70 CONTINUE RETURN END C --------------------------------------------------------------------- SUBROUTINE transp(a,ia,n) INTEGER ia,n DOUBLE PRECISION a(ia,n) DOUBLE PRECISION dummy DO 10 i = 1,n DO 20 j = i+1,n dummy = a(i,j) a(i,j) = a(j,i) a(j,i) = dummy 20 CONTINUE 10 CONTINUE RETURN END C ----------------------------------------------------------------- SUBROUTINE lncomb(a,ia,b,ib,c,ic,n,facb,facc) INTEGER ia,ib,ic,n DOUBLE PRECISION a(ia,n),b(ib,n),c(ic,n) DOUBLE PRECISION facb,facc DO 10 i = 1,n DO 20 j = 1,n a(j,i) = facb*b(j,i)+facc*c(j,i) 20 CONTINUE 10 CONTINUE RETURN END C ------------------------------------------------------------------ DOUBLE PRECISION FUNCTION spexp(x) C .. Scalar Arguments .. DOUBLE PRECISION x C .. C .. Local Scalars .. DOUBLE PRECISION eneg LOGICAL first C .. C .. External Functions .. DOUBLE PRECISION x02amf EXTERNAL x02amf C .. C .. Intrinsic Functions .. INTRINSIC exp,log C .. C .. Save statement .. SAVE first,eneg C .. C .. Data statements .. DATA first/.true./ C .. IF (first) THEN eneg = log(x02amf(x)) first = .false. END IF spexp = x02amf(x) IF (x.GT.eneg) spexp = exp(x) RETURN END C --------------------------------------------------------------------- DOUBLE PRECISION FUNCTION phi(v) C .. Parameters .. DOUBLE PRECISION a0,a1,a2,a3,a4,a5 PARAMETER (a0=1.D0,a1=a0/6.D0,a2=a1/20.D0,a3=a2/42.D0,a4=a3/72.D0, & a5=a4/110.D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION v C .. phi = a0 + v* (a1+v* (a2+v* (a3+v* (a4+v*a5)))) RETURN END C --------------------------------------------------------------------- DOUBLE PRECISION FUNCTION psi(v) C .. Parameters .. DOUBLE PRECISION b0,b1,b2,b3,b4,b5 PARAMETER (b0=1.D0,b1=b0/2.D0,b2=b1/12.D0,b3=b2/30.D0, & b4=b3/56.D0,b5=b4/90.D0) C .. C .. Scalar Arguments .. DOUBLE PRECISION v C .. psi = b0 + v* (b1+v* (b2+v* (b3+v* (b4+v*b5)))) RETURN END C --------------------------------------------------------------------- DOUBLE PRECISION FUNCTION scl(th,u) C .. Parameters .. DOUBLE PRECISION one PARAMETER (one=1.d0) C .. C .. Scalar Arguments .. DOUBLE PRECISION th,u C .. C .. Intrinsic Functions .. INTRINSIC atan2,cos,sin C .. scl = th + atan2((u-one)*sin(th)*cos(th),one+ (u-one)*sin(th)**2) RETURN END