SUBROUTINE SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,WORK, & IWORK,LWK,ILWK,IFAIL) C Specification: SL11F solves Hamiltonian eigenvalue problems which are C reducible to the form of a first order matrix differential equation C Theta' = i Theta Omega(x,lambda,Theta), C where Theta is a unitary matrix prescribed by boundary conditions at C x=a and x=b. C C Variables: C elam: double precision. On entry elam should be set to an approximation C to the eigenvalue to be found. The value of elam is not important C in determining the accuracy of the final result but it may affect C run times. On exit, elam specifies the approximation to the C eigenvalue computed by SL11F. 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 C lies within [elam-C*eps,elam+C*eps] where by elam we mean the C initial approximation provided by the user, and C is a constant C of modest size. Again the value of eps is not critical. In cases C where absolutely no prior information is available it would be C acceptable in almost all cases to set elam = 0 and eps = 10 on C entry. C a,b: double precision. On entry a and b specify the finite regular endpoints C of the interval on which the problem is posed. Unchanged 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 C of 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. It is C intended that the smaller tol is, the smaller will be the error in C the final approximation elam returned by the code. The value of C tol is unchanged on exit. C knobs: integer. Used to control the maximum number of miss-distance C evaluations carried out by the code. The number of miss-distance C evaluations will be max(60,knobs). knobs is unchanged on exit. C work: double precision array of dimension iwork. Used as workspace. C iwork: integer. On entry iwork must be at least n*(16+n*(33+2*nsampl))+25+ C 2*nsampl. Unchanged on exit. C lwk: integer array of dimension ilwk. Used as workspace. C ilwk: integer. On entry ilwk must be at least n. Unchanged on exit. 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 USUB1 and USUB2. C C USUB1: C C SUBROUTINE usub1(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 USUB2: C C SUBROUTINE usub2(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 USUB2 must specify the boundary condition at C x=a, while if iend=1 then USUB2 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 USUB2. 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 NSAMPL PARAMETER (NSAMPL=20) C .. C .. Local Scalars .. DOUBLE PRECISION C,EL,EU,TOLC05 INTEGER IFLAG,IFO,IREST,ISAMPL,ISIGIL,ISIGIR,ISIGRL,ISIGRR,KNTMAX, & NREC CHARACTER*6 SRNAME C .. C .. Scalar Arguments .. DOUBLE PRECISION A,B,ELAM,EPS,TOL INTEGER IFAIL,ILWK,IWORK,K,KNOBS,N C .. C .. Array Arguments .. C IWORK must be at least n*(16+(33+2*nsampl)*n)+25+2*nsampl. With nsampl=10 C this gives n*(16+53*n)+45, with nsampl=20 we get n*(16+73*n)+65 DOUBLE PRECISION WORK(IWORK) INTEGER LWK(ILWK) C .. C .. Subroutine Arguments .. EXTERNAL USUB1,USUB2 C .. C .. External Subroutines .. EXTERNAL CHOOSC,SOLVSY C .. C .. Intrinsic Functions .. INTRINSIC MAX,SQRT C .. C .. External Functions .. DOUBLE PRECISION X02AJF INTEGER P01ABF EXTERNAL X02AJF,P01ABF C .. C .. Local Arrays .. CHARACTER*80 REC(2) C .. SRNAME = 'SL11F' 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* (16+ (33+2*NSAMPL)*N)+25+2*NSAMPL .OR. & ILWK.LT.N .OR. EPS.LE.0.0D0) THEN C We have an input error IFLAG = 1 IF (A.GE.B) THEN WRITE (REC,FMT=9000) A,B NREC = 2 GO TO 10 ELSE IF (N.LT.2) THEN WRITE (REC,FMT=9010) N NREC = 2 GO TO 10 ELSE IF (TOL.LE.0.0D0) THEN WRITE (REC,FMT=9020) TOL NREC = 2 GO TO 10 ELSE IF (IWORK.LT.N* (16+ (33+2*NSAMPL)*N)+25+2*NSAMPL) THEN WRITE (REC,FMT=9030) N* (16+ (33+2*NSAMPL)*N) + 25 + & 2*NSAMPL,IWORK NREC = 2 GO TO 10 ELSE IF (ILWK.LT.N) THEN WRITE (REC,FMT=9040) N,ILWK NREC = 2 GO TO 10 ELSE IF (EPS.LE.0.0D0) THEN WRITE (REC,FMT=9050) EPS NREC = 2 GO TO 10 ELSE END IF END IF ISIGRL = 1 ISIGIL = N*N + 1 ISIGRR = 2*N*N + 1 ISIGIR = 3*N*N + 1 IREST = 4*N*N + 1 ISAMPL = 40*N*N + 4*N + 1 KNTMAX = MAX(60,KNOBS) TOLC05 = 1.0D-2 IFAIL = 0 C = 0.50D0* (A+B) CALL SOLVSY(A,B,C,K,ELAM,EPS,EL,EU,USUB1,USUB2,TOL,TOLC05, & WORK(IREST),WORK(ISIGRL),WORK(ISIGIL),N,WORK(ISIGRR), & WORK(ISIGIR),N,LWK,N,KNTMAX,IFAIL) IF (IFAIL.EQ.0) THEN C Choose a better matching point and repeat CALL CHOOSC(A,B,C,EL,EU,TOL,USUB1,USUB2,WORK,LWK,N,K,NSAMPL, & WORK(N* (16+33*N)+25+2*NSAMPL+1),IFAIL) IF (IFAIL.EQ.0) THEN TOLC05 = MAX(1.0D-1*TOL,SQRT(X02AJF(0.0D0))) ELAM = 0.50D0* (EU+EL) EPS = 0.50D0* (EU-EL) CALL SOLVSY(A,B,C,K,ELAM,EPS,EL,EU,USUB1,USUB2,TOL,TOLC05, & WORK(IREST),WORK(ISIGRL),WORK(ISIGIL),N, & WORK(ISIGRR),WORK(ISIGIR),N,LWK,N,KNTMAX, & IFAIL) IF (IFAIL.EQ.0) THEN RETURN END IF END IF END IF IFLAG = IFAIL IF (IFAIL.EQ.7) THEN WRITE (REC,FMT=9060) NREC = 2 IFLAG = 2 ELSE IF (IFAIL.EQ.12 .OR. IFAIL.EQ.14) THEN WRITE (REC,FMT=9070) NREC = 2 IFLAG = 3 ELSE IF (IFAIL.EQ.13 .OR. IFAIL.EQ.15) THEN WRITE (REC,FMT=9080) NREC = 2 IFLAG = 4 ELSE IF (IFAIL.EQ.10 .OR. IFAIL.EQ.11 .OR. IFAIL.EQ.16 .OR. & IFAIL.EQ.17) THEN WRITE (REC,FMT=9090) NREC = 2 IFLAG = 5 ELSE IF (IFAIL.EQ.9 .OR. IFAIL.EQ.18) THEN WRITE (REC,FMT=9100) NREC = 2 IFLAG = 6 ELSE IF (IFAIL.EQ.20) THEN WRITE (REC,FMT=9110) NREC = 2 IFLAG = 7 ELSE IF (IFAIL.EQ.21) THEN WRITE (REC,FMT=9120) NREC = 2 IFLAG = 8 ELSE IF (IFAIL.EQ.19) THEN WRITE (REC,FMT=9130) NREC = 2 IFLAG = 9 END IF 10 IFAIL = P01ABF(IFO,IFLAG,SRNAME,NREC,REC) 9000 FORMAT (' ** A must be less than B, you have',/,' ** A = ',g18.8, & ' B = ',g18.8) 9010 FORMAT (' ** N must be at least 2, you have',/,' ** N = ',i8) 9020 FORMAT (' ** TOL must be strictly positive, you have',/, & ' ** TOL = ',g18.8) 9030 FORMAT (' ** IWORK must be at least ',i8,' you have',/, & ' ** IWORK = ',i8) 9040 FORMAT (' ** ILWK must be at least',i8,' you have',/, & ' ** ILWK = ',i8) 9050 FORMAT (' ** EPS must be strictly positive, you have',/, & ' ** EPS = ',g18.8) 9060 FORMAT (' ** Too many attempts have been made to find',/, & ' ** required eigenvalue; try increasing EPS') 9070 FORMAT (' ** Failure in F02BJF when trying to convert boundary',/, & ' ** conditions from user variables to H-matrix variables') 9080 FORMAT (' ** The boundary conditions supplied do not',/, & ' ** satisfy the rank condition') 9090 FORMAT (' ** An H-matrix has occurred whose exponential',/, & ' ** cannot be accurately computed') 9100 FORMAT (' ** A failure has occurred on a call to F02AJF;',/, & ' ** cannot find eigenvalues of central miss-matrix') 9110 FORMAT (' ** An H-matrix has arisen which cannot be',/, & ' ** diagonalised; try reducing TOL') 9120 FORMAT (' ** Too many unsuccessful attempts have been made',/, & ' ** to find an initial stepsize') 9130 FORMAT (' ** Too many unsuccessful attempts have been made',/, & ' ** to find a suitable next stepsize') END C ---------------------------------------------------------------------- SUBROUTINE SOLVSY(A,B,C,K,ELAM,EPS,EL,EU,COFMAT,SYSBC,TOL,TOLC05, & WORK,SIGRL,SIGIL,ISIGL,SIGRR,SIGIR,ISIGR,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 A,B,C,EL,ELAM,EPS,EU,TOL,TOLC05 INTEGER IFAIL,ISIGL,ISIGR,K,KNTMAX,N C .. C .. Array Arguments .. DOUBLE PRECISION SIGIL(ISIGL,N),SIGIR(ISIGR,N),SIGRL(ISIGL,N), & SIGRR(ISIGR,N),WORK(N* (13+25*N)+25) INTEGER LWK(N) C .. C .. Subroutine Arguments .. EXTERNAL COFMAT,SYSBC C .. C .. Local Scalars .. DOUBLE PRECISION DL,DU,E1,E2,SHIFT,TOLOC,E3,F3 INTEGER IFLAG,IND,KNT LOGICAL DUAS,TREVAL C .. C .. Local Arrays .. DOUBLE PRECISION C1(17) C .. C .. External Functions .. DOUBLE PRECISION DLH13F EXTERNAL DLH13F C .. C .. External Subroutines .. C EXTERNAL C05AZF EXTERNAL C05TST C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,MAX,MIN C .. SHIFT = DBLE(K) + 5.0D-1 EL = ELAM - EPS EU = ELAM + EPS DUAS = .false. KNT = 0 IFAIL = 0 10 CONTINUE C WRITE (6,FMT=*) 'el=',el DL = DLH13F(A,B,C,EL,N,TOL,SIGRL,ISIGL,SIGRR,ISIGR,WORK(1), & WORK(N*N+1),WORK(2*N*N+1),WORK(3*N*N+1),WORK(4*N*N+1),LWK, & COFMAT,SYSBC,SHIFT,IFAIL) IF (IFAIL.NE.0) THEN RETURN C WRITE (6,FMT=*) 'elam,dmiss:',el,dl ELSE IF (DL.GT.0.0D0) THEN EU = EL DU = DL DUAS = .true. EPS = 2.0D0*EPS EL = EL - EPS GO TO 10 END IF IF (.NOT.DUAS) THEN DU = DLH13F(A,B,C,EU,N,TOL,SIGRL,ISIGL,SIGRR,ISIGR,WORK(1), & WORK(N*N+1),WORK(2*N*N+1),WORK(3*N*N+1),WORK(4*N*N+1), & LWK,COFMAT,SYSBC,SHIFT,IFAIL) IF (IFAIL.NE.0) THEN RETURN ELSE 20 CONTINUE C WRITE (6,FMT=*) 'elam,dmiss:',eu,du IF (DL*DU.GT.0.0D0) THEN C We have not yet bracketed an eigenvalue KNT = KNT + 1 IF (KNT.LE.KNTMAX) THEN EPS = 2.0D0*EPS IF (DU.LT.0.0D0) THEN EL = EU DL = DU EU = EU + EPS DU = DLH13F(A,B,C,EU,N,TOL,SIGRL,ISIGL,SIGRR, & ISIGR,WORK(1),WORK(N*N+1),WORK(2*N*N+1), & WORK(3*N*N+1),WORK(4*N*N+1),LWK,COFMAT, & SYSBC,SHIFT,IFAIL) IF (IFAIL.NE.0) RETURN ELSE C WRITE (6,FMT=*) 'elam,dmiss:',eu,du EU = EL DU = DL EL = EL - EPS DL = DLH13F(A,B,C,EL,N,TOL,SIGRL,ISIGL,SIGRR, & ISIGR,WORK(1),WORK(N*N+1),WORK(2*N*N+1), & WORK(3*N*N+1),WORK(4*N*N+1),LWK,COFMAT, & SYSBC,SHIFT,IFAIL) IF (IFAIL.NE.0) THEN RETURN ELSE C WRITE (6,FMT=*) 'elam,dmiss:',el,dl END IF END IF GO TO 20 END IF ELSE GO TO 30 END IF IFAIL = 7 RETURN END IF 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 = TOLC05 C1(1) = DU IFLAG = 1 E1 = EL E2 = EU TREVAL = .FALSE. 40 CONTINUE C CALL C05AZF(E1,E2,DL,TOLOC,1,C1,IND,IFLAG) CALL C05TST(E1,E2,E3,F3,TREVAL,DL,TOLOC,1,C1,IND,IFLAG) IF (IND.EQ.0) THEN GO TO 50 ELSE KNT = KNT + 1 IF (KNT.LE.KNTMAX) THEN DL = DLH13F(A,B,C,E1,N,TOL,SIGRL,ISIGL,SIGRR,ISIGR, & WORK(1),WORK(N*N+1),WORK(2*N*N+1),WORK(3*N*N+1), & WORK(4*N*N+1),LWK,COFMAT,SYSBC,SHIFT,IFAIL) C WRITE (6,FMT=*) 'elam,dmiss:',E1,DL IF (IFAIL.NE.0) THEN RETURN ELSE 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))) THEN GO TO 50 ELSE C WRITE (6,FMT=*) 'elam,dmiss:',E1,DL GO TO 40 END IF END IF END IF END IF IFAIL = 7 RETURN 50 ELAM = E1 END C ------------------------------------------------------------------------ SUBROUTINE c05tst(x,y,z,fz,treval,fx,tolx,ir,c,ind,ifail) C First written by Chris Graves (ESF Student, Leicester University) C on 13/10/92. Modified to cope with `nasty' functions by switching C to bisection, by Marco Marletta on 7th June 1993. Specification C similar to NAG routine C05AZF. DOUBLE PRECISION x,y,z,fx,fz,tolx,c(17) INTEGER ir,ind,ifail LOGICAL treval DOUBLE PRECISION x0,x1,valx0,valx1,new,new2,x1old SAVE x0,x1,valx0,valx1,new GOTO (100,110,120,130,140,150) ind+2 C (IND=-1) 100 x0=x x1=y z = x0 valx0=fx valx1=c(1) fz = valx0 GOTO 150 c (IND=0) 110 RETURN c (IND=1) 120 x0=x x1=y ind=2 RETURN c (IND=2) 130 valx0=fx z = x fz = fx x=x1 ind=3 RETURN c (IND=3) 140 valx1=fx c (IND=4) 150 new=x1-(valx1*(x1-x0)/(valx1-valx0)) IF (treval) THEN new2=x1-(valx1*(x1-z)/(valx1-fz)) ELSE new2 = new END IF IF (sign(1.D0,valx0)*valx1.LT.0.D0) THEN z = x0 fz = valx0 END IF treval = .true. x0 = x1 valx0=valx1 IF (abs(new-new2).LT.0.25D0*abs(x1-z)) THEN x1=new ELSE x1 = 0.5D0*(x1+z) END IF c ..CHECK TOLERANCE.. IF (ABS(x1-z).le.tolx*max(1.D0,x0,x1)) THEN ind=0 RETURN END IF x=x1 ind=3 RETURN END C ------------------------------------------------------------------------ SUBROUTINE CHOOSC(A,B,C,EL,EU,TOL,COFMAT,SYSBC,WORK,LWK,N,K, & NSAMPL,WSAMPL,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 .. C .. C .. Scalar Arguments .. DOUBLE PRECISION A,B,C,EL,EU,TOL INTEGER IFAIL,K,N,NSAMPL C .. C .. Array Arguments .. DOUBLE PRECISION WORK(N* (16+33*N)+25+2*NSAMPL), & WSAMPL(1:N,1:N,1:2*NSAMPL) INTEGER LWK(N) C .. C .. Subroutine Arguments .. EXTERNAL COFMAT,SYSBC C .. C .. Local Scalars .. DOUBLE PRECISION ADJUST,ALFMAX,ALFMIN,ALPHA,ARGC,ARGL,ARGR,BETA, & DIFF,DLA13F,ELAM,FAC,HSAMPL,OMEGA,SHIFT,TWOPI,X1, & X2 INTEGER I,IALFI,IALFR,IASMP,IBETA,IDLSMP,IHL,II,INDEX,ISAMPI, & ISAMPR,IUL,IULO,IUR,IV,IVL,IVLO,IVR,J,JMATCH,JSAMPL,KNT,N2 LOGICAL ISING C .. C .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF C .. C .. External Subroutines .. EXTERNAL F02AJF,F02BJF,F06QFF,F06QHF,F06YAF,HINT,MATEXP,MVN,UDUT C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,DBLE,SIGN C .. IFAIL = 0 SHIFT = DBLE(K) + 5.0D-1 TWOPI = 2.0D0*X01AAF(0.0D0) N2 = N*N C Set up pointers for the workspace: IUL = N* (13+25*N) + 25 + 1 IVL = IUL + N2 IUR = IVL + N2 IVR = IUR + N2 IULO = IVR + N2 IVLO = IULO + N2 IALFR = IVLO + N2 IALFI = IALFR + N IBETA = IALFI + N IV = IBETA + N IHL = IV + N2 IASMP = IHL + N2 IDLSMP = IASMP + NSAMPL C Need workspace of length n*(16+33*n)+25+2*nsampl ELAM = EL DIFF = 100.0D0 10 CONTINUE C Set up initial condition at x=a: CALL SYSBC(0,ISING,A,WORK(IUL),WORK(IVL),N,N,ELAM) C Now we convert the boundary condition information into information C about the H-matrix. CALL F06QFF('G',N,N,WORK(IUL),N,WORK(IULO),N) CALL F06QFF('G',N,N,WORK(IVL),N,WORK(IVLO),N) 20 CONTINUE CALL F02BJF(N,WORK(IULO),N,WORK(IVLO),N,0.0D0,WORK(IALFR), & WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,IFAIL) IF (IFAIL.NE.0) THEN KNT = KNT + 1 IFAIL = 1 IF (KNT.LT.4) THEN GO TO 20 ELSE GO TO 230 END IF END IF C Now retrieve the eigenvectors. DO 30 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN GO TO 220 ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN C Form the eigenvector as Uz and normalise CALL MVN(WORK(IUL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N) ELSE C Form the eigenvector Vz and normalise CALL MVN(WORK(IVL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N) END IF 30 CONTINUE C Finished retrieval of eigenvectors; now get phase angles. Remember the C funny things that ATAN2 can do in IEEE arithmetic. DO 40 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ALPHA.EQ.0.0D0) THEN WORK(IALFR+J-1) = 0.0D0 ELSE FAC = SIGN(1.0D0,ALPHA) ALPHA = FAC*ALPHA BETA = FAC*BETA WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA) END IF 40 CONTINUE C C Now we are ready to form the initial matrix H. Remember that: C the eigenvectors are in the matrix work(iv+j),j=0,..,n*n-1; C the eigenvalues are in the elements work(ialfr+j),j=0,..,n. C The following routine call will overwrite the eigenvectors but C not the eigenvalues. Also, notice that work(iul+j), j=0,..,n*n-1 C is used here purely for storage. C CALL UDUT(WORK(IV),N,WORK(IALFR),WORK(IHL),N,WORK(IULO),N) C Do case of lower eigenvalue approximation first: C Integrate forward, 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: HSAMPL = (B-A)/DBLE(NSAMPL) X1 = A X2 = X1 + HSAMPL ISAMPR = 1 ISAMPI = N*N + 1 DO 60 JSAMPL = 1,NSAMPL - 1 IFAIL = 0 CALL HINT(X1,X2,ELAM,WORK(IHL),N,N,ADJUST,TOL,WORK(1),COFMAT, & IFAIL) IF (IFAIL.NE.0) THEN RETURN ELSE C Store the outcome of the integration: CALL F06QFF('G',N,N,WORK(IHL),N,WSAMPL(1,1,JSAMPL),N) ARGL = 0.0D0 INDEX = IHL DO 50 II = 1,N ARGL = ARGL + WORK(INDEX) INDEX = INDEX + N + 1 50 CONTINUE ARGL = ARGL + ADJUST WORK(IASMP+JSAMPL-1) = ARGL ISAMPR = ISAMPR + N*N X1 = X2 X2 = X2 + HSAMPL END IF 60 CONTINUE C Set up initial condition at x=b: CALL SYSBC(1,ISING,B,WORK(IUL),WORK(IVL),N,N,ELAM) CALL F06QFF('G',N,N,WORK(IUL),N,WORK(IULO),N) CALL F06QFF('G',N,N,WORK(IVL),N,WORK(IVLO),N) 70 CONTINUE CALL F02BJF(N,WORK(IULO),N,WORK(IVLO),N,0.0D0,WORK(IALFR), & WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,IFAIL) IF (IFAIL.NE.0) THEN KNT = KNT + 1 IFAIL = 1 IF (KNT.LT.4) THEN GO TO 70 ELSE GO TO 210 END IF END IF C Now retrieve the eigenvectors. DO 80 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN GO TO 200 ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN C Form the eigenvector as Uz and normalise CALL MVN(WORK(IUL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N) ELSE C Form the eigenvector Vz and normalise CALL MVN(WORK(IVL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N) END IF 80 CONTINUE C Finished retrieval of eigenvectors; now get phase angles. Remember the C funny things that ATAN2 can do in IEEE arithmetic. DO 90 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ALPHA.EQ.0.0D0) THEN WORK(IALFR+J-1) = 0.0D0 ELSE FAC = SIGN(1.0D0,ALPHA) ALPHA = FAC*ALPHA BETA = FAC*BETA WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA) END IF 90 CONTINUE C CALL UDUT(WORK(IV),N,WORK(IALFR),WORK(IHL),N,WORK(IULO),N) C Integrate backwards: X1 = B HSAMPL = (A-B)/DBLE(NSAMPL) X2 = X1 + HSAMPL JMATCH = NSAMPL DO 160 JSAMPL = NSAMPL - 1,1,-1 IFAIL = 0 CALL HINT(X1,X2,ELAM,WORK(IHL),N,N,ADJUST,TOL,WORK(1),COFMAT, & IFAIL) IF (IFAIL.NE.0) THEN RETURN ELSE INDEX = IHL ARGR = 0.0D0 DO 100 II = 1,N ARGR = ARGR + WORK(INDEX) INDEX = INDEX + N + 1 100 CONTINUE ARGR = ARGR + ADJUST C C Form the miss matrix at the current point: C CALL F06QFF('G',N,N,WSAMPL(1,1,JSAMPL),N,WORK(IUL),N) ISAMPR = JSAMPL*N*N + 1 CALL F06QHF('G',N,N,0.0D0,0.0D0,WORK(IVL),N) CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,WORK(IULO), & WORK(IVLO),N,WORK(1),IFAIL) IF (IFAIL.NE.0) THEN GO TO 190 ELSE C That does exp(i.hl) INDEX = IHL DO 120 J = 1,N DO 110 I = 1,N WORK(INDEX) = -WORK(INDEX) INDEX = INDEX + 1 110 CONTINUE 120 CONTINUE CALL F06QFF('G',N,N,WORK(IHL),N,WORK(IUL),N) CALL F06QHF('G',N,N,0.0D0,0.0D0,WORK(IVL),N) CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,WORK(IUR), & WORK(IVR),N,WORK(1),IFAIL) IF (IFAIL.NE.0) THEN GO TO 180 ELSE C Now form C exp(-i.hr).exp(i.hl) = (ur + i.vr)*(ul + i.vl) CALL F06YAF('N','N',N,N,N,1.0D0,WORK(IUR),N, & WORK(IULO),N,0.0D0,WORK(IUL),N) CALL F06YAF('N','N',N,N,N,-1.0D0,WORK(IVR),N, & WORK(IVLO),N,1.0D0,WORK(IUL),N) CALL F06YAF('N','N',N,N,N,1.0D0,WORK(IVR),N, & WORK(IULO),N,0.0D0,WORK(IVL),N) CALL F06YAF('N','N',N,N,N,1.0D0,WORK(IUR),N, & WORK(IVLO),N,1.0D0,WORK(IVL),N) IFAIL = 1 CALL F02AJF(WORK(IUL),N,WORK(IVL),N,N,WORK(1), & WORK(1+N),LWK,IFAIL) IF (IFAIL.NE.0) THEN GO TO 170 ELSE C Now we get the phase angles and store them in work(ialfr+j) for C j=0,..,n-1. ARGC = 0.0D0 ALFMAX = -100.0D0 ALFMIN = 100.0D0 DO 130 J = 1,N ALPHA = WORK(N+J) BETA = WORK(J) OMEGA = ATAN2(ALPHA,BETA) IF (OMEGA.LT.0.0D0) OMEGA = OMEGA + TWOPI WORK(J) = OMEGA ARGC = ARGC + OMEGA IF (OMEGA.GT.ALFMAX) ALFMAX = OMEGA IF (OMEGA.LT.ALFMIN) ALFMIN = OMEGA 130 CONTINUE IF (ELAM.EQ.EU) DLA13F = ALFMIN/TWOPI IF (ELAM.EQ.EL) DLA13F = ALFMAX/TWOPI - 1.0D0 IF (ELAM.EQ.EL) THEN WORK(IDLSMP+JSAMPL-1) = DLA13F ELSE C WRITE (11,FMT=*) 'Stored miss distance:',dla13f C WRITE (11,FMT=*) 'elam=',elam 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:',dble(jsampl)* (b-a)/ C & dble(nsampl) + a,dla13f - work(idlsmp+jsampl-1) IF (ABS(DLA13F-WORK(IDLSMP+JSAMPL-1)).LE. & DIFF) THEN JMATCH = JSAMPL DIFF = ABS(DLA13F- & WORK(IDLSMP+JSAMPL-1)) END IF END IF X1 = X2 X2 = X2 + HSAMPL INDEX = IHL DO 150 J = 1,N DO 140 I = 1,N WORK(INDEX) = -WORK(INDEX) INDEX = INDEX + 1 140 CONTINUE 150 CONTINUE END IF END IF END IF END IF 160 CONTINUE IF (ELAM.EQ.EL) THEN ELAM = EU GO TO 10 END IF C = A + DBLE(JMATCH)* (B-A)/DBLE(NSAMPL) C WRITE (6,FMT=*) 'Match-point chosen=',c RETURN 170 IFAIL = 9 RETURN 180 IFAIL = 10 RETURN 190 IFAIL = 11 RETURN 200 IFAIL = 15 RETURN 210 IFAIL = 14 RETURN 220 IFAIL = 13 RETURN 230 IFAIL = 12 C$st$ Unreachable comments ... END C ------------------------------------------------------------------------ SUBROUTINE MATSUM(TYPE,A,IA,FACA,B,IB,FACB,C,IC,FACC,N) C .. Scalar Arguments .. DOUBLE PRECISION FACA,FACB,FACC INTEGER IA,IB,IC,N CHARACTER TYPE C .. C .. Array Arguments .. DOUBLE PRECISION A(IA,N),B(IB,N),C(IC,N) C .. C .. Local Scalars .. INTEGER I,J C .. IF (FACC.NE.0.0D0) THEN IF (TYPE.EQ.'S') THEN DO 20 J = 1,N DO 10 I = J,N C(I,J) = FACC*C(I,J) + FACA*A(I,J) + FACB*B(I,J) C(J,I) = C(I,J) 10 CONTINUE 20 CONTINUE ELSE IF (TYPE.EQ.'A') THEN DO 40 J = 1,N DO 30 I = J + 1,N C(I,J) = FACC*C(I,J) + FACA*A(I,J) + FACB*B(I,J) C(J,I) = -C(I,J) 30 CONTINUE C(J,J) = 0.D0 40 CONTINUE ELSE DO 60 J = 1,N DO 50 I = 1,N C(I,J) = FACC*C(I,J) + FACA*A(I,J) + FACB*B(I,J) 50 CONTINUE 60 CONTINUE END IF C Special coding to take account of the fact that the array c might be C un-initialised in this case. ELSE IF (TYPE.EQ.'S') THEN DO 80 J = 1,N DO 70 I = J,N C(I,J) = FACA*A(I,J) + FACB*B(I,J) C(J,I) = C(I,J) 70 CONTINUE 80 CONTINUE ELSE IF (TYPE.EQ.'A') THEN DO 100 J = 1,N DO 90 I = J + 1,N C(I,J) = FACA*A(I,J) + FACB*B(I,J) C(J,I) = -C(I,J) 90 CONTINUE C(J,J) = 0.D0 100 CONTINUE ELSE DO 120 J = 1,N DO 110 I = 1,N C(I,J) = FACA*A(I,J) + FACB*B(I,J) 110 CONTINUE 120 CONTINUE END IF END C ------------------------------------------------------------------------ SUBROUTINE HSTEP(X1,X2,TOL,FSAL,ELAM,H,IH,N,COFMAT,ERR,WORK,PRN, & IFAIL) C .. Parameters .. C .. C .. Scalar Arguments .. DOUBLE PRECISION ELAM,ERR,TOL,X1,X2 INTEGER IFAIL,IH,N LOGICAL FSAL,PRN C .. C .. Array Arguments .. DOUBLE PRECISION H(IH,N),WORK(N* (2+18*N)) C .. C .. Subroutine Arguments .. EXTERNAL COFMAT C .. C .. Local Scalars .. DOUBLE PRECISION ELMAX,FAC1,FAC10,FAC11,FAC12,FAC13,FAC14,FAC15, & FAC16,FAC17,FAC18,FAC19,FAC2,FAC20,FAC21,FAC22, & FAC23,FAC3,FAC4,FAC5,FAC6,FAC7,FAC8,FAC9,HSTP,HN INTEGER COFF11,COFF12,COFF21,COFF22,HDOT1,HDOT2,HDOT3,HDOT4,HDOT5, & HDOT6,HINDEX,IE,IR,N2,S11,S12,S21,S22,STORE,V C .. C .. External Subroutines .. EXTERNAL F06QFF,HDER,MATSUM,MONIT C .. C .. External Functions .. DOUBLE PRECISION F06QGF EXTERNAL F06QGF C .. C .. Intrinsic Functions .. INTRINSIC ABS,REAL C .. N2 = N*N HINDEX = 1 HDOT1 = HINDEX + N2 HDOT2 = HDOT1 + N2 HDOT3 = HDOT2 + N2 HDOT4 = HDOT3 + N2 HDOT5 = HDOT4 + N2 HDOT6 = HDOT5 + N2 S11 = HDOT6 + N2 S12 = S11 + N2 S21 = S12 + N2 S22 = S21 + N2 COFF11 = S22 + N2 COFF12 = COFF11 + N2 COFF21 = COFF12 + N2 COFF22 = COFF21 + N2 STORE = COFF22 + N2 V = STORE + N2 IR = V + N2 IE = IR + N IEND = N* (2+17*N) + 1 C Need workspace of length (18*n+2)*n HN = F06QGF('M','S',N,N,H,IH) C Store current H: CALL F06QFF('G',N,N,H,IH,WORK(HINDEX),N) C Get hdot1: IF (FSAL) THEN CALL F06QFF('G',N,N,WORK(IEND),N,WORK(HDOT1),N) ELSE CALL HDER(X1,ELAM,WORK(HINDEX),N,WORK(HDOT1),N,N,WORK(S11), & WORK(S12),WORK(S21),WORK(S22),WORK(COFF11), & WORK(COFF12),WORK(COFF21),WORK(COFF22),WORK(STORE), & WORK(V),N,WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN END IF C Force a stepsize reduction if hdot is too large. IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT1),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN ENDIF C Get hdot2: HSTP = (X2-X1)/5.D0 CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,HSTP,WORK(HINDEX),N, & 0.0D0,N) CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT2),N,N, & WORK(S11),WORK(S12),WORK(S21),WORK(S22), & WORK(COFF11),WORK(COFF12),WORK(COFF21),WORK(COFF22), & WORK(STORE),WORK(V),N,WORK(IR),WORK(IE),COFMAT, & ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT2),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN END IF C Get hdot3: HSTP = 3.D0* (X2-X1)/10.D0 CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,0.25D0*HSTP,WORK(STORE), & N,0.0D0,N) CALL MATSUM('S',WORK(STORE),N,1.0D0,WORK(HDOT2),N,HSTP*0.75D0, & WORK(HINDEX),N,0.0D0,N) CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT3),N,N,WORK(S11), & WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12), & WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N, & WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT3),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN END IF C Get hdot4: HSTP = (X2-X1) FAC1 = 44.D0/45.D0 FAC2 = -56.D0/15.D0 FAC3 = 32.D0/9.D0 CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,FAC1*HSTP,WORK(HINDEX), & N,0.0D0,N) CALL MATSUM('S',WORK(HDOT2),N,FAC2*HSTP,WORK(HDOT3),N,FAC3*HSTP, & WORK(HINDEX),N,1.D0,N) HSTP = 0.8D0*HSTP CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT4),N,N,WORK(S11), & WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12), & WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N, & WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT4),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN END IF C Get hdot5: HSTP = X2 - X1 FAC4 = 19372.D0/6561.D0 FAC5 = -25360.D0/2187.D0 FAC6 = 64448.D0/6561.D0 FAC7 = -212.D0/729.D0 CALL F06QFF('G',N,N,H,IH,WORK(HINDEX),N) CALL MATSUM('S',WORK(HDOT1),N,FAC4*HSTP,WORK(HDOT2),N,FAC5*HSTP, & WORK(HINDEX),N,1.D0,N) CALL MATSUM('S',WORK(HDOT3),N,FAC6*HSTP,WORK(HDOT4),N,FAC7*HSTP, & WORK(HINDEX),N,1.D0,N) HSTP = HSTP*8.D0/9.D0 CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT5),N,N,WORK(S11), & WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12), & WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N, & WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT5),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN END IF C Get hdot6: HSTP = X2 - X1 FAC8 = 9017.D0/3168.D0 FAC9 = -355.D0/33.D0 FAC10 = 46732.D0/5247.D0 FAC11 = 49.D0/176.D0 FAC12 = -5103.D0/18656.D0 CALL F06QFF('G',N,N,H,IH,WORK(HINDEX),N) CALL MATSUM('S',WORK(HDOT1),N,FAC8*HSTP,WORK(HDOT2),N,FAC9*HSTP, & WORK(HINDEX),N,1.D0,N) CALL MATSUM('S',WORK(HDOT3),N,FAC10*HSTP,WORK(HDOT4),N,FAC11*HSTP, & WORK(HINDEX),N,1.D0,N) CALL MATSUM('S',WORK(HINDEX),N,1.D0,WORK(HDOT5),N,FAC12*HSTP, & WORK(S11),N,0.D0,N) CALL F06QFF('G',N,N,WORK(S11),N,WORK(HINDEX),N) CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT6),N,N,WORK(S11), & WORK(S12),WORK(S21),WORK(S22),WORK(COFF11), & WORK(COFF12),WORK(COFF21),WORK(COFF22),WORK(STORE), & WORK(V),N,WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT5),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN END IF C Get hdot7 (store it in WORK(HDOT2) which is no longer required C at this point): HSTP = X2 - X1 FAC13 = 35.D0/384.D0 FAC14 = 500.D0/1113.D0 FAC15 = 125.D0/192.D0 FAC16 = -2187.D0/6784.D0 FAC17 = 11.D0/84.D0 CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,FAC13*HSTP,WORK(HINDEX), & N,0.0D0,N) CALL MATSUM('S',WORK(HDOT3),N,FAC14*HSTP,WORK(HDOT4),N,FAC15*HSTP, & WORK(HINDEX),N,1.0D0,N) CALL MATSUM('S',WORK(HDOT5),N,FAC16*HSTP,WORK(HDOT6),N,FAC17*HSTP, & WORK(HINDEX),N,1.0D0,N) CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT2),N,N,WORK(S11), & WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12), & WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N, & WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN IF (ABS((X2-X1)*F06QGF('M','S',N,N,WORK(HDOT5),N)).GT.1.0D1) THEN ERR = TOL*32.D0 RETURN END IF C Now store the old value of H: CALL F06QFF('G',N,N,H,IH,WORK(S11),N) C Now update h: CALL F06QFF('G',N,N,WORK(HINDEX),N,H,IH) C Now store the old value of H: CALL F06QFF('G',N,N,WORK(S11),N,WORK(HINDEX),N) C Now form the error estimate HSTP = (X2-X1) FAC18 = 5179.D0/57600.D0 FAC19 = 7571.D0/16695.D0 FAC20 = 393.D0/640.D0 FAC21 = -92097.D0/339200.D0 FAC22 = 187.D0/2100.D0 FAC23 = 1.D0/40.D0 CALL MATSUM('S',WORK(HDOT1),N,FAC18*HSTP,WORK(HDOT3),N,FAC19*HSTP, & WORK(HINDEX),N,1.D0,N) CALL MATSUM('S',WORK(HDOT4),N,FAC20*HSTP,WORK(HDOT5),N,FAC21*HSTP, & WORK(HINDEX),N,1.D0,N) CALL MATSUM('S',WORK(HDOT6),N,FAC22*HSTP,WORK(HDOT2),N,FAC23*HSTP, & WORK(HINDEX),N,1.D0,N) CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT3),N,N, & WORK(S11),WORK(S12),WORK(S21),WORK(S22), & WORK(COFF11),WORK(COFF12),WORK(COFF21), & WORK(COFF22),WORK(STORE),WORK(V),N, & WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL) IF (IFAIL.NE.0) RETURN CALL MATSUM('S',WORK(HINDEX),N,1.D0,H,IH,-1.D0,WORK(S12),N,0.D0,N) CALL MATSUM('S',WORK(HDOT2),N,1.D0,WORK(HDOT3),IH,-1.D0, & WORK(HDOT1),N,0.D0,N) ERR = SQRT( (F06QGF('M','S',N,N,WORK(S12),N)**2 & + F06QGF('M','S',N,N,WORK(HDOT1),N)**2)/(1.D0+HN**2) ) CALL F06QFF('G',N,N,WORK(HDOT2),N,WORK(IEND),N) C Now return to the calling routine 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,MAXEL,S INTEGER I,INDEX,IR,IVI,IVIS,IVR,IVRS,IWK1,IWK2,IWK3,J C .. C .. External Subroutines .. EXTERNAL F02AXF,F06YAF C .. C .. External Functions .. DOUBLE PRECISION F06QGF EXTERNAL F06QGF C .. C .. Intrinsic Functions .. INTRINSIC COS,EXP,MAX,SIN C .. MAXEL = MAX(1.0D0,F06QGF('M','G',N,N,AR,IA)) MAXEL = MAX(MAXEL,F06QGF('M','G',N,N,AI,IA)) 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 DO 20 J = 1,N DO 10 I = 1,N EXPR(I,J) = AR(I,J)/MAXEL EXPI(I,J) = AI(I,J)/MAXEL 10 CONTINUE 20 CONTINUE IFAIL = -1 CALL F02AXF(EXPR,IEXP,EXPI,IEXP,N,WORK(IR),WORK(IVR),N, & WORK(IVI),N,WORK(IWK1),WORK(IWK2),WORK(IWK3), & IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 1 RETURN ELSE C Now get the exponential of the eigenvalues: INDEX = IR DO 30 J = 1,N WORK(INDEX) = EXP(MAXEL*WORK(INDEX)) INDEX = INDEX + 1 30 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 50 J = 1,N DO 40 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 40 CONTINUE 50 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) END IF ELSE C expr + i*expi = (rr+i*ri)*exp(D)*(transpose(rr)-i*transpose(ri)) 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 70 J = 1,N DO 60 I = 1,N EXPR(I,J) = -AR(I,J)/MAXEL EXPI(I,J) = AI(I,J)/MAXEL 60 CONTINUE AR(J,J) = 0.0D0 EXPR(J,J) = 0.0D0 70 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(EXPI,IEXP,EXPR,IEXP,N,WORK(IR),WORK(IVR),N, & WORK(IVI),N,WORK(IWK1),WORK(IWK2),WORK(IWK3), & IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 1 RETURN ELSE CC Now we can return the matrix A to its original form: CC DO 70 j = 1,n C DO 60 i = 1,n C ar(i,j) = -ar(i,j) C60 CONTINUE C70 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(MAXEL*WORK(IR+J-1)) S = SIN(MAXEL*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 END IF IFAIL = 0 C$st$ Unreachable comments ... END C ------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DLH13F(A,B,C,ELAM,N,TOL,HL,IHL,HR,IHR, & UL,VL,UR,VR,WORK,LWK,COFMAT,SYSBC,SHIFT,IFAIL) C This subroutine computes the M(lambda) miss distance function C .. Scalar Arguments .. DOUBLE PRECISION A,B,C,ELAM,SHIFT,TOL INTEGER IFAIL,IHL,IHR,N C .. C .. Array Arguments .. DOUBLE PRECISION HL(IHL,N),HR(IHR,N),UL(N,N),UR(N,N),VL(N,N), & VR(N,N),WORK(N* (13+21*N)+25) INTEGER LWK(N) C .. C .. Subroutine Arguments .. EXTERNAL COFMAT,SYSBC C .. C .. Local Scalars .. DOUBLE PRECISION ADJL,ADJR,ALFMAX,ALFMIN,ALPHA,ARGL,ARGR,BETA,FAC, & OMEGA,SUM,TWOPI INTEGER I,IALFI,IALFR,IBETA,IUL,IV,IVL,IWK,J,KNT,N2 LOGICAL ISING C .. C .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF C .. C .. External Subroutines .. EXTERNAL F02AJF,F02BJF,F06QFF,F06QHF,F06YAF,HINT,MATEXP,MVN,UDUT C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,DBLE,SIGN C .. TWOPI = 2.0D0*X01AAF(0.0D0) N2 = N*N C WRITE (6,FMT=*) '**** In dlh13f, elam=',elam C Set pointers for workspace array: IUL = 1 IVL = IUL + N2 IALFR = IVL + N2 IALFI = IALFR + N IBETA = IALFI + N IV = IBETA + N IWK = IV + N2 C Set up boundary condition at x=a: CALL SYSBC(0,ISING,A,UL,VL,N,N,ELAM) C Now we must transform this into a boundary condition on H by getting C the eigenvalues and eigenvectors of the Theta-matrix associated with C these boundary conditions. We do this by solving the simultaneous C eigenvalue problem C (cos(omega/2)U-sin(omega/2)V)z = 0, C then recovering the eigenvectors w from a formula such as C w = (V-iU)z. C CALL F06QFF('G',N,N,UL,N,WORK(IUL),N) CALL F06QFF('G',N,N,VL,N,WORK(IVL),N) IFAIL = 1 KNT = 0 10 CONTINUE CALL F02BJF(N,WORK(IUL),N,WORK(IVL),N,0.0D0,WORK(IALFR), & WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,IFAIL) IF (IFAIL.NE.0) THEN KNT = KNT + 1 IFAIL = 1 IF (KNT.LT.4) THEN GO TO 10 ELSE GO TO 140 END IF END IF C Now retrieve the eigenvectors. DO 20 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN GO TO 130 ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN C Form the eigenvector as Uz and normalise CALL MVN(UL,N,WORK(IV+ (J-1)*N),WORK(IALFI),N) ELSE C Form the eigenvector Vz and normalise CALL MVN(VL,N,WORK(IV+ (J-1)*N),WORK(IALFI),N) END IF 20 CONTINUE C Finished retrieval of eigenvectors; now get phase angles. Remember the C funny things that ATAN2 can do in IEEE arithmetic. DO 30 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ALPHA.EQ.0.0D0) THEN WORK(IALFR+J-1) = 0.0D0 ELSE FAC = SIGN(1.0D0,ALPHA) ALPHA = FAC*ALPHA BETA = FAC*BETA WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA) END IF 30 CONTINUE C C Now we are ready to form the initial matrix H. Remember that: C the eigenvectors are in the matrix work(iv+j),j=0,..,n*n-1; C the eigenvalues are in the elements work(ialfr+j),j=0,..,n. C The following routine call will overwrite the eigenvectors but C not the eigenvalues. Also, notice that work(iul+j), j=0,..,n*n-1 C is used here purely for storage. C CALL UDUT(WORK(IV),N,WORK(IALFR),HL,IHL,WORK(IUL),N) C Now do the integration from x=a to x=c: CALL HINT(A,C,ELAM,HL,IHL,N,ADJL,TOL,WORK(1),COFMAT,IFAIL) IF (IFAIL.EQ.0) THEN C That does the left-hand leg of the integration. Now do the right C hand leg. C Set up boundary condition at x=b: CALL SYSBC(1,ISING,B,UR,VR,N,N,ELAM) C Now we must transform this into a boundary condition on H by getting C the eigenvalues and eigenvectors of the Theta-matrix associated with C these boundary conditions. We do this by solving the simultaneous C eigenvalue problem C (cos(omega/2)U-sin(omega/2)V)z = 0, C then recovering the eigenvectors w from a formula such as C w = (V-iU)z. C IFAIL = 1 KNT = 0 CALL F06QFF('G',N,N,UR,N,WORK(IUL),N) CALL F06QFF('G',N,N,VR,N,WORK(IVL),N) 40 CONTINUE CALL F02BJF(N,WORK(IUL),N,WORK(IVL),N,0.0D0,WORK(IALFR), & WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK, & IFAIL) IF (IFAIL.NE.0) THEN KNT = KNT + 1 IFAIL = 1 IF (KNT.LT.4) THEN GO TO 40 ELSE GO TO 120 END IF END IF C Now retrieve the eigenvectors. DO 50 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN GO TO 110 ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN C Form the eigenvector as Uz and normalise CALL MVN(UR,N,WORK(IV+ (J-1)*N),WORK(IALFI),N) ELSE C Form the eigenvector Vz and normalise CALL MVN(VR,N,WORK(IV+ (J-1)*N),WORK(IALFI),N) END IF 50 CONTINUE C Finished retrieval of eigenvectors; now get phase angles. Remember the C funny things that ATAN2 can do in IEEE arithmetic. DO 60 J = 1,N ALPHA = WORK(IALFR+J-1) BETA = WORK(IBETA+J-1) IF (ALPHA.EQ.0.0D0) THEN WORK(IALFR+J-1) = TWOPI ELSE FAC = SIGN(1.0D0,ALPHA) ALPHA = FAC*ALPHA BETA = FAC*BETA WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA) END IF 60 CONTINUE C C Now we are ready to form the initial matrix H. Remember that: C the eigenvectors are in the matrix work(iv+j),j=0,..,n*n-1; C the eigenvalues are in the elements work(ialfr+j),j=0,..,n. C The following routine call will overwrite the eigenvectors but C not the eigenvalues. Also, notice that work(iul+j), j=0,..,n*n-1 C is used here purely for storage. C CALL UDUT(WORK(IV),N,WORK(IALFR),HR,IHR,WORK(IUL),N) C Now do the integration from x=b to x=c: CALL HINT(B,C,ELAM,HR,IHR,N,ADJR,TOL,WORK(1),COFMAT,IFAIL) IF (IFAIL.EQ.0) THEN C C Now we have the left-hand and right-hand H-matrices at the centre. C We must compute the phase-angles of the matrix exp(-i.hr).exp(i.hl). C First we must compute the two matrices concerned, since the eigenvalues C are not those of the matrix hl-hr. C First get the traces since we will probably overwrite the matrices C hl and hr and we need to know their traces. ARGL = ADJL ARGR = ADJR DO 70 I = 1,N ARGL = ARGL + HL(I,I) ARGR = ARGR + HR(I,I) 70 CONTINUE CALL F06QFF('G',N,N,HL,IHL,WORK(IUL),N) CALL F06QHF('G',N,N,0.0D0,0.0D0,WORK(IVL),N) CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,UL,VL,N,WORK(IWK), & IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 16 ELSE C That does exp(i.hl) DO 90 J = 1,N DO 80 I = 1,N HR(I,J) = -HR(I,J) 80 CONTINUE 90 CONTINUE CALL F06QFF('G',N,N,HR,IHR,WORK(IUL),N) CALL F06QHF('G',N,N,0.0D0,0.0D0,WORK(IVL),N) CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,UR,VR,N, & WORK(IWK),IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 17 ELSE C Now form C exp(-i.hr).exp(i.hl) = (ur + i.vr)*(ul + i.vl) CALL F06YAF('N','N',N,N,N,1.0D0,UR,N,UL,N,0.0D0, & WORK(IUL),N) CALL F06YAF('N','N',N,N,N,-1.0D0,VR,N,VL,N,1.0D0, & WORK(IUL),N) CALL F06YAF('N','N',N,N,N,1.0D0,VR,N,UL,N,0.0D0, & WORK(IVL),N) CALL F06YAF('N','N',N,N,N,1.0D0,UR,N,VL,N,1.0D0, & WORK(IVL),N) C We have got the real and imaginary parts; now we can find the C eigenvalues. For this we use the NAG routine f02ajf IFAIL = 1 CALL F02AJF(WORK(IUL),N,WORK(IVL),N,N,WORK(IWK), & WORK(IWK+N),LWK,IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 18 ELSE C Now we get the phase angles and store them in work(ialfr+j) for C j=0,..,n-1. SUM = 0.0D0 ALFMAX = -100.0D0 ALFMIN = 100.0D0 DO 100 J = 1,N ALPHA = WORK(IWK+N+J-1) BETA = WORK(IWK+J-1) OMEGA = ATAN2(ALPHA,BETA) IF (OMEGA.LT.0.0D0) OMEGA = OMEGA + TWOPI WORK(IALFR+J-1) = OMEGA SUM = SUM + OMEGA IF (OMEGA.GT.ALFMAX) ALFMAX = OMEGA IF (OMEGA.LT.ALFMIN) ALFMIN = OMEGA 100 CONTINUE C WRITE (6,FMT=*) 'argl,argr,argc:',argl,argr,sum DLH13F = (ARGL-ARGR-SUM)/TWOPI + DBLE(N) - & SHIFT C WRITE (6,FMT=*) ' From DLH13F: ELAM,DMISS:',ELAM,DLH13F IF (ABS(DLH13F).LE.0.6D0) THEN IF (DLH13F.GT.0.0D0) THEN DLH13F = ALFMIN/TWOPI ELSE DLH13F = ALFMAX/TWOPI - 1.0D0 END IF END IF END IF END IF END IF END IF RETURN 110 IFAIL = 15 RETURN 120 IFAIL = 14 END IF RETURN 130 IFAIL = 13 RETURN 140 IFAIL = 12 C$st$ Unreachable comments ... END C ------------------------------------------------------------------------ SUBROUTINE UDUT(U,IU,D,OUT,IOUT,STORE,N) C .. Scalar Arguments .. INTEGER IOUT,IU,N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),OUT(IOUT,N),STORE(N,N),U(IU,N) C .. C .. Local Scalars .. INTEGER I,J C .. C .. External Subroutines .. EXTERNAL F06YAF C .. DO 20 J = 1,N DO 10 I = 1,N STORE(I,J) = U(I,J)*D(J) 10 CONTINUE 20 CONTINUE CALL F06YAF('N','T',N,N,N,1.0D0,STORE,N,U,IU,0.0D0,OUT,IOUT) END C ------------------------------------------------------------------------- SUBROUTINE MVN(A,IA,V,STORE,N) C .. Scalar Arguments .. INTEGER IA,N C .. C .. Array Arguments .. DOUBLE PRECISION A(IA,N),STORE(N),V(N) C .. C .. Local Scalars .. DOUBLE PRECISION ELMAX,SUM INTEGER I,J C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,SQRT C .. DO 10 I = 1,N STORE(I) = 0.0D0 10 CONTINUE DO 30 J = 1,N DO 20 I = 1,N STORE(I) = STORE(I) + A(I,J)*V(J) 20 CONTINUE 30 CONTINUE ELMAX = 0.0D0 DO 40 I = 1,N ELMAX = MAX(ELMAX,ABS(STORE(I))) 40 CONTINUE SUM = 0.0D0 DO 50 I = 1,N STORE(I) = STORE(I)/ELMAX SUM = SUM + STORE(I)**2 50 CONTINUE SUM = SQRT(SUM) DO 60 I = 1,N V(I) = STORE(I)/SUM 60 CONTINUE END C ------------------------------------------------------------------------- SUBROUTINE HINT(XO,XEND,ELAM,H,IH,N,ADJUST,TOL,WORK,COFMAT,IFAIL) C .. Scalar Arguments .. DOUBLE PRECISION ADJUST,ELAM,TOL,XEND,XO INTEGER IFAIL,IH,N C .. C .. Array Arguments .. DOUBLE PRECISION H(IH,N),WORK(N* (21*N+13)+25) C .. C .. Subroutine Arguments .. EXTERNAL COFMAT C .. C .. Local Scalars .. DOUBLE PRECISION ALFA,BETA,BND,ERR,GAP,GAPL,HL,HMAX,PI,RATIO, & TOLOC,TWOPI,U,X1,X2 INTEGER I,IE,IERR,IFLAG,IGAP,IGOOD,IHLOC,INDEX,IR,ITRY,IV,IWORK, & NSTEPS LOGICAL PHASE1,PRN,FSAL C .. C .. External Functions .. DOUBLE PRECISION X01AAF,X02AJF EXTERNAL X01AAF,X02AJF C .. C .. Parameters .. DOUBLE PRECISION SAFE,SMALL,FIFTH PARAMETER (SAFE=0.80D0,SMALL=0.70D0,FIFTH=0.20D0) C .. C .. External Subroutines .. EXTERNAL F02ABF,F06QFF,HSTEP,M01CAF,UDUT C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,MAX,MIN,NINT,SIGN C .. PRN = .false. FSAL = .FALSE. PI = X01AAF(0.0D0) TWOPI = 2.0D0*X01AAF(0.0D0) HL = TOL*SIGN(1.0D0,XEND-XO) HMAX = ABS(XEND-XO)/4.0D0 ADJUST = 0.0D0 TOLOC = MIN(0.1D0,TOL) NSTEPS = 0 ITRY = 0 IERR = 0 IGOOD = 4 IHLOC = 1 IE = IHLOC + N*N IV = IE + N IR = IV + N*N IWORK = IHLOC + N*N PHASE1 = .true. C Phase1 controls the special step-control algorithm which is used only C on the first step to get the meshing `on scale'. X1 = XO 10 CONTINUE CALL F06QFF('G',N,N,H,IH,WORK(IHLOC),N) C Before doing any stepping we subtract off a matrix M which commutes with C H and is such that exp(iM) = I. This will avoid singularities later. IFLAG = 1 CALL F02ABF(WORK(IHLOC),N,N,WORK(IE),WORK(IV),N,WORK(IR),IFLAG) IF (IFLAG.NE.0) THEN GO TO 70 ELSE CALL F06QFF('G',N,1,WORK(IE),N,WORK(IR),N) INDEX = IR DO 20 I = 1,N ALFA = TWOPI*DBLE(NINT(WORK(INDEX)/TWOPI)) WORK(INDEX) = WORK(INDEX) - ALFA IF (WORK(INDEX).LT.0.0D0) WORK(INDEX) = WORK(INDEX) + & TWOPI INDEX = INDEX + 1 20 CONTINUE C Now order the eigenvalues IFLAG = 0 CALL M01CAF(WORK(IR),1,N,'A',IFLAG) C Now find the greatest gap GAP = -WORK(IR+N-1) + WORK(IR) + TWOPI IGAP = N DO 30 I = 1,N - 1 GAPL = WORK(IR+I) - WORK(IR+I-1) IF (GAPL.GT.GAP) THEN GAP = GAPL IGAP = I END IF 30 CONTINUE C Have now found the greatest gap. IF (IGAP.LT.N) THEN BETA = 0.50D0* (WORK(IR+IGAP)+WORK(IR+IGAP-1)) ELSE BETA = 0.50D0* (WORK(IR)+WORK(IR+N-1)-TWOPI) END IF IF (BETA.GE.PI) BETA = BETA - TWOPI C Now normalise eigenvalues so that those which correspond to close points on the C unit circle are indeed close DO 40 I = 1,N ALFA = TWOPI*DBLE(NINT(WORK(IE+I-1)/TWOPI)) U = WORK(IE+I-1) - ALFA IF (U.GE.TWOPI+BETA) THEN U = U - TWOPI ELSE IF (U.LT.BETA) THEN U = U + TWOPI END IF ADJUST = ADJUST + WORK(IE+I-1) - U WORK(IE+I-1) = U 40 CONTINUE C WRITE (6,FMT=*) 'Eigenvalue:',I,WORK(IE+I-1) C Now construct the new matrix H: CALL UDUT(WORK(IV),N,WORK(IE),WORK(IHLOC),N,WORK(IR),N) CALL F06QFF('G',N,N,WORK(IHLOC),N,H,IH) C Use the routine hstep to do the stepping X2 = X1 + HL C WRITE (6,FMT=*) 'H X1 X2 ELAM',X2-X1,X1,X2,ELAM IF (ABS(XEND-X1).LT.1.D1*X02AJF(0.0D0)) THEN GO TO 60 ELSE C C WARNING: The last N**2 elements of the array WORK must be C preserved between one call to HSTEP and the next. C CALL HSTEP(X1,X2,TOLOC,FSAL,ELAM,WORK(IHLOC),N,N,COFMAT, & ERR,WORK(IWORK),PRN,IFAIL) IF (IFAIL.NE.0) THEN RETURN ELSE FSAL = .FALSE. IF (ERR.GT.TOLOC) THEN C WRITE (7,FMT=*) 'Failed step' IERR = IERR + 1 IGOOD = 0 IF (IERR.GT.16 .AND. .NOT.PHASE1) THEN GO TO 50 ELSE IF (.NOT.PHASE1 .OR. IERR.LE.64) THEN RATIO = MIN((TOLOC/ERR)**FIFTH,SAFE) HL = HL*RATIO GO TO 10 END IF ELSE IF (ERR.GT.TOLOC*SAFE) THEN RATIO = (TOLOC*SAFE/ERR)**FIFTH HL = HL*RATIO ELSE IF (ERR.LT.SMALL*TOLOC .AND. IERR.EQ.0) THEN IGOOD = IGOOD + 1 RATIO = (ERR/TOLOC)**FIFTH BND = MAX(FIFTH,1.0D0/ (DBLE(IGOOD+1))) HL = HL/MAX(BND,RATIO) IF (ABS(HL).GT.HMAX) HL = SIGN(1.0D0,HL)*HMAX IF (PHASE1 .AND. ITRY.LT.8 .AND. & ABS(HL).LT.HMAX) THEN ITRY = ITRY + 1 C WRITE (7,FMT=*) 'Failed step' GO TO 10 END IF END IF FSAL = .TRUE. C WRITE (7,FMT=*) 'Successful step' IERR = 0 PHASE1 = .false. IF (ABS(X2-X1).GE.1.0D-3) PRN = .false. X1 = X2 IF ((XEND-X2-HL)*SIGN(1.0D0,HL).LE. & 0.0D0) HL = XEND - X1 C Update h CALL F06QFF('G',N,N,WORK(IHLOC),N,H,IH) NSTEPS = NSTEPS + 1 GO TO 10 END IF END IF END IF END IF IFAIL = 21 RETURN 50 IFAIL = 19 RETURN C WRITE (6,FMT=*) 'No. of steps:',nsteps 60 IFAIL = 0 RETURN 70 IFAIL = 20 END C ---------------------------------------------------------------------- SUBROUTINE HDER(X,ELAM,H,IH,HDOT,IHDOT,N,S11,S12,S21,S22,COFF11, & COFF12,COFF21,COFF22,STORE,V,IV,R,E,COFMAT,ELMAX, & IFAIL) C .. Scalar Arguments .. DOUBLE PRECISION ELAM,ELMAX,X INTEGER IFAIL,IH,IHDOT,IV,N C .. C .. Array Arguments .. DOUBLE PRECISION COFF11(N,N),COFF12(N,N),COFF21(N,N),COFF22(N,N), & E(N),H(IH,N),HDOT(IHDOT,N),R(N),S11(N,N), & S12(N,N),S21(N,N),S22(N,N),STORE(N,N),V(IV,N) C .. C .. Subroutine Arguments .. EXTERNAL COFMAT C .. C .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF C .. C .. External Subroutines .. EXTERNAL F02ABF,F06QFF,F06YAF,FGH C .. C Step 1: find the eigenvalues and eigenvectors of H: C .. Local Scalars .. DOUBLE PRECISION ELOC,TWOPI,F,G,GT,HL INTEGER I,IFLAG,J,P,Q C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX C .. IFLAG = 1 CALL F02ABF(H,IH,N,E,V,IV,R,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = 20 ELSE TWOPI = 2.0D0*X01AAF(0.0D0) C Eigenvalues and eigenvectors now found. Now suppose that Theta = C exp(iH), and Theta' = i.Theta.Omega. From the known expression for Omega C we can compute dH/dt. CALL COFMAT(X,ELAM,S11,S12,S21,S22,N) ELOC = 0.0D0 DO 20 J = 1,N DO 10 I = 1,N ELOC = MAX(ELOC,ABS(S11(I,J)),ABS(S12(I,J)), & ABS(S21(I,J)),ABS(S22(I,J))) 10 CONTINUE 20 CONTINUE ELMAX = MAX(ELOC,ELMAX) CALL F06YAF('T','N',N,N,N,1.0D0,V,IV,S11,N,0.0D0,STORE,N) CALL F06YAF('N','N',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,COFF11,N) CALL F06YAF('T','N',N,N,N,1.0D0,V,IV,S12,N,0.0D0,STORE,N) CALL F06YAF('N','N',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,COFF12,N) CALL F06YAF('T','N',N,N,N,1.0D0,V,IV,S22,N,0.0D0,STORE,N) CALL F06YAF('N','N',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,COFF22,N) DO 40 P = 1,N DO 30 Q = 1,N COFF21(Q,P) = COFF12(P,Q) 30 CONTINUE 40 CONTINUE DO 60 Q = 1,N DO 50 P = 1,Q CALL FGH(E(P),E(Q),F,G,GT,HL) STORE(P,Q) = 0.5D0*( F*COFF11(P,Q)+G*COFF12(P,Q)+ & GT*COFF21(P,Q)+HL*COFF22(P,Q) ) STORE(Q,P) = STORE(P,Q) 50 CONTINUE 60 CONTINUE CALL F06QFF('G',N,N,STORE,N,COFF11,N) CALL F06YAF('N','N',N,N,N,1.0D0,V,IV,COFF11,N,0.0D0,STORE,N) CALL F06YAF('N','T',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,HDOT,IHDOT) C We have now computed the derivative required. END IF END C ---------------------------------------------------------------------- SUBROUTINE ANSYM(A,IA,B,IB,N,C,IC) C .. Scalar Arguments .. INTEGER IA,IB,IC,N C .. C .. Array Arguments .. DOUBLE PRECISION A(IA,N),B(IB,N),C(IC,N) C .. C .. Local Scalars .. INTEGER I,J,K C .. DO 20 J = 1,N DO 10 I = 1,N C(I,J) = 0.D0 10 CONTINUE 20 CONTINUE DO 30 J = 1,N DO 40 I = 1,J DO 50 K = 1,N C(I,J) = C(I,J) + A(I,K)*B(K,J) 50 CONTINUE 40 CONTINUE 30 CONTINUE C DO 50 K = 1,N C DO 40 J = 1,N C DO 30 I = 1,J C C(I,J) = C(I,J) + A(I,K)*B(K,J) C 30 CONTINUE C 40 CONTINUE C 50 CONTINUE DO 70 J = 1,N DO 60 I = 1,J - 1 C(J,I) = C(I,J) 60 CONTINUE 70 CONTINUE END C ---------------------------------------------------------------------- SUBROUTINE FGH(X,Y,F,G,GT,H) DOUBLE PRECISION X,Y,F,G,GT,H DOUBLE PRECISION T,PHI,RATIO,DIFF,SUM,COSDIF,COSSUM,SINDIF,SINSUM EXTERNAL PHI DIFF = (X-Y)*0.5D0 SUM = DIFF + Y T = DIFF**2 SINDIF = SIN(DIFF) IF (T.LT.0.1D0) THEN C RATIO = 1.D0/PHI(-T) RATIO = PHI(-T) ELSE C RATIO = DIFF/(SINDIF+SIGN(1.D0,SINDIF)*1.0D-14) RATIO = (SINDIF+SIGN(1.D0,SINDIF)*1.0D-14)/DIFF END IF C RATIO = 2.D0*RATIO RATIO = 0.5D0*RATIO COSDIF = COS(DIFF) COSSUM = COS(SUM) SINSUM = SIN(SUM) C F = RATIO*(COSDIF-COSSUM) C G = RATIO*(SINDIF+SINSUM) C GT = RATIO*(SINSUM-SINDIF) C H = RATIO*(COSDIF+COSSUM) F = (COSDIF-COSSUM)/RATIO G = (SINDIF+SINSUM)/RATIO GT = (SINSUM-SINDIF)/RATIO H = (COSDIF+COSSUM)/RATIO 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 TRACE(A,IA,N) C .. Scalar Arguments .. INTEGER IA,N C .. C .. Array Arguments .. DOUBLE PRECISION A(IA,N) C .. C .. Local Scalars .. DOUBLE PRECISION SUM INTEGER J C .. SUM = 0.D0 DO 10 J = 1,N SUM = SUM + A(J,J) 10 CONTINUE TRACE = SUM 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 END