SUBROUTINE SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,DMISS, & CONST,IFAIL) C .. Parameters .. INTEGER IU,IV PARAMETER (IU=2,IV=2) C .. C .. Scalar Arguments .. DOUBLE COMPLEX DMISS,ELAM DOUBLE PRECISION CONST INTEGER IFAIL,IMATCH,NMESH C .. C .. Array Arguments .. DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE PRECISION CONSTL,CONSTR INTEGER IEND C .. C .. Local Arrays .. DOUBLE COMPLEX U(IU,IU),V(IV,IV),ZL(6),ZR(6) C .. C .. External Subroutines .. EXTERNAL EXPINT C .. IF (IMATCH.LT.0 .OR. IMATCH.GT.NMESH) THEN WRITE (6,FMT=*) 'IMATCH,NMESH:',IMATCH,NMESH IFAIL = 4 RETURN END IF C Get boundary conditions at the left hand endpoint: IEND = 0 CALL SL4BCS(IEND,XMESH(0),ELAM,U,IU,V,IV) C Convert to the Riccati form: ZL(1) = U(1,1)*U(2,2) - U(2,1)*U(1,2) ZL(2) = V(1,1)*V(2,2) - V(2,1)*V(1,2) ZL(3) = U(1,1)*V(2,2) - U(1,2)*V(2,1) ZL(4) = -U(1,1)*V(1,2) + U(1,2)*V(1,1) ZL(5) = U(2,1)*V(2,2) - U(2,2)*V(2,1) ZL(6) = -U(2,1)*V(1,2) + U(2,2)*V(1,1) CONSTL = 0.D0 C Integrate to the matchpoint: CALL EXPINT(ZL,CONSTL,XMESH,NMESH,0,IMATCH,SL4COF,ELAM,IFAIL) IF (IFAIL.NE.0) RETURN WRITE (12,FMT=*) 'ELAM,X,DETU,DETV:',ELAM,XMESH(IMATCH), & ZL(1)*EXP(CONSTL),ZL(2)*EXP(CONSTL) C Get boundary conditions at the right hand endpoint: IEND = 1 CALL SL4BCS(IEND,XMESH(NMESH),ELAM,U,IU,V,IV) C Convert to the Riccati form: ZR(1) = U(1,1)*U(2,2) - U(2,1)*U(1,2) ZR(2) = V(1,1)*V(2,2) - V(2,1)*V(1,2) ZR(3) = U(1,1)*V(2,2) - U(1,2)*V(2,1) ZR(4) = -U(1,1)*V(1,2) + U(1,2)*V(1,1) ZR(5) = U(2,1)*V(2,2) - U(2,2)*V(2,1) ZR(6) = -U(2,1)*V(1,2) + U(2,2)*V(1,1) CONSTR = 0.D0 C Integrate to the matchpoint: CALL EXPINT(ZR,CONSTR,XMESH,NMESH,NMESH,IMATCH,SL4COF,ELAM,IFAIL) IF (IFAIL.NE.0) RETURN C Now form the miss distance function: this is the C determinant of the matrix /U_L U_R\ C \V_L V_R/ DMISS = ZL(1)*ZR(2) + ZL(2)*ZR(1) - ZL(3)*ZR(6) + ZL(4)*ZR(5) - & ZL(6)*ZR(3) + ZL(5)*ZR(4) CONST = CONSTL + CONSTR WRITE (8,FMT=*) 'ELAM,F,CONST:',ELAM,DMISS,CONST IFAIL = 0 RETURN END C ------------------------------------------------------------------- SUBROUTINE EXPINT(Z,CONST,XMESH,NMESH,ISTART,IEND,SL4COF,ELAM, & IFAIL) C C This routine integrates the linear system for the 4th order C Riccati variables across a mesh in each interval of which the C coefficients in the 4th order ODE are constant. C C Z: on entry, the array of Riccati variables at XMESH(ISTART) C on exit, the Riccati variables at XMESH(IEND) scaled by C exp(-CONST) to avoid overflow. C C CONST: on entry, CONST = 0. On exit, CONST specifies the number C such that exp(CONST).Z recovers the Riccati variables. N.B. C computation of exp(CONST) may cause overflow. C C XMESH: input. The mesh, on each interval of which the coefficients C in the ODE are assumed to be constant. C C SL4COF: specifies the coefficients in the ODE, which is C C (P y'')'' + (Q3 y'')' + (Q2 y')' + Q1 y' + Q0 y = 0. C C The dependence on the eigenparameter ELAM is hidden in the C coefficients. C C IFAIL: usual NAG-type error flag. C C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION CONST INTEGER IEND,IFAIL,ISTART,NMESH C .. C .. Array Arguments .. DOUBLE COMPLEX Z(6) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX CONE,CZERO DOUBLE PRECISION CNLOC,ZLMAX INTEGER I,IDI,II C .. C .. Local Arrays .. DOUBLE COMPLEX F(6,6),ZLOC(6) C .. C .. External Subroutines .. EXTERNAL EXPSOL,F06SAF,F06TFF C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX,LOG,MAX,MIN C .. CONE = CMPLX(1.D0,0.D0) CZERO = CMPLX(0.D0,0.D0) IF (MIN(ISTART,IEND).LT.0 .OR. MAX(ISTART,IEND).GT.NMESH) THEN C Cannot integrate over this range IFAIL = 1 RETURN END IF CONST = 0.D0 IDI = 1 IF (IEND.LT.ISTART) IDI = -1 IF (ISTART.EQ.IEND) RETURN DO 30 I = ISTART + IDI,IEND,IDI CNLOC = 0.D0 IFAIL = 0 CALL EXPSOL(SL4COF,ELAM,XMESH(I-IDI),XMESH(I),CNLOC,F,6,IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 2 RETURN END IF CONST = CONST + CNLOC CALL F06SAF('N',6,6,CONE,F,6,Z,1,CZERO,ZLOC,1) CALL F06TFF('G',6,1,ZLOC,6,Z,6) C Normalise to get the Riccati variables to a reasonable scale ZLMAX = 0.D0 DO 10 II = 1,6 ZLMAX = MAX(ZLMAX,ABS(Z(II))) 10 CONTINUE DO 20 II = 1,6 Z(II) = Z(II)/ZLMAX 20 CONTINUE CONST = CONST + LOG(ZLMAX) 30 CONTINUE IFAIL = 0 RETURN END C ------------------------------------------------------------------- SUBROUTINE EXPSOL(SL4COF,ELAM,XSTART,XEND,CONST,F,IF,IFAIL) C C The Wronskians in the 4th order system satisfy a linear system of C the form z' = Az, where A is a sparse 6x6 matrix. Here we C attempt to solve this system. C C The equation being solved is C C (P y'')'' + (Q3 y'')' + (Q2 y')' + Q1 y' + Q0 y = 0 C C Routine to compute exp(A*(xend-xstart)) using routine from Golub and Van Loan C The output is exp(const).F; const is a real constant, a scaling used C to avoid overflow. C C .. Parameters .. INTEGER IX,IXOUT,ID,IN PARAMETER (IX=6,IXOUT=IX,ID=IX,IN=IX) C .. C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION CONST,XEND,XSTART INTEGER IF,IFAIL C .. C .. Array Arguments .. DOUBLE COMPLEX F(IF,6) C .. C .. Local Scalars .. DOUBLE COMPLEX DIAG,EQR,FACC,MIN1K,OFDIAG,TERM DOUBLE PRECISION DUMMY,FAC,ERREST INTEGER II,IJ,IQ,J,K C .. C .. Local Arrays .. DOUBLE COMPLEX D(ID,IX),N(IN,IX),X(IX,IX),X1(IX,IX),XOUT(IXOUT,IX) DOUBLE PRECISION WREAL(IX) C .. C .. External Functions .. DOUBLE PRECISION F06UAF,X02AJF EXTERNAL F06UAF,X02AJF C .. C .. External Subroutines .. EXTERNAL F04ADF,F06TFF,F06THF,F06ZAF,GETSIG C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX,DBLE,INT,LOG,MAX C .. C .. Subroutine Arguments .. EXTERNAL SL4COF C .. OFDIAG = CMPLX(0.D0,0.D0) DIAG = CMPLX(1.D0,0.D0) CALL F06THF('G',6,6,OFDIAG,DIAG,X,IX) CALL GETSIG(SL4COF,ELAM,XSTART,XEND,XOUT,IXOUT,ERREST) C XOUT contains the matrix A (not in sparse form). F06UAF computes its C sup norm. C F06UAF here gets the (infinity) norm of a matrix J = MAX(0,1+INT(LOG(F06UAF('I',6,6,XOUT,IXOUT,WREAL)/LOG(2.D0)))) FAC = 0.5D0**J FACC = CMPLX(FAC,0.D0) C WRITE (8,FMT=*) 'J,FAC:',J,FAC IQ = 1 EQR = CMPLX(1.D0/6.D0,0.D0) 10 EQR = EQR/CMPLX(DBLE(16*IQ*IQ* (2*IQ+1)* (2*IQ+3)),0.D0) C WRITE (8,FMT=*) 'E_q:',EQR IQ = IQ + 1 IF (ABS(EQR).GT.X02AJF(1.D0)) GO TO 10 C WRITE (6,FMT=*) 'Q is',IQ C Initialise matrices D, N and X all to be equal to the identity CALL F06THF('G',6,6,OFDIAG,DIAG,X,IX) CALL F06THF('G',6,6,OFDIAG,DIAG,N,IN) CALL F06THF('G',6,6,OFDIAG,DIAG,D,ID) CONST = 1.D0 DO 40 K = 1,IQ CONST = CONST* (DBLE(IQ-K+1)/DBLE(IQ+IQ-K+1))/DBLE(K) CALL F06ZAF('N','N',6,6,6,FACC,XOUT,IXOUT,X,IX,OFDIAG,X1,IX) CALL F06TFF('G',6,6,X1,IX,X,IX) MIN1K = (-1.D0)**K DO 30 IJ = 1,6 DO 20 II = 1,6 TERM = CONST*X(II,IJ) N(II,IJ) = N(II,IJ) + TERM D(II,IJ) = D(II,IJ) + MIN1K*TERM 20 CONTINUE 30 CONTINUE 40 CONTINUE C Solve DF = N for F using a suitable routine IFAIL = 0 C CALL F04ADF(D,ID,N,IN,6,6,F,IF,WREAL,IFAIL) C For k = 1 to j do F = F*F. N.B. F is almost certainly C full. C Normalise F: CONST = MAX(1.D0,F06UAF('M',6,6,F,IF,WREAL)) DO 60 IJ = 1,6 DO 50 II = 1,6 F(II,IJ) = F(II,IJ)/CONST 50 CONTINUE 60 CONTINUE CALL F06TFF('G',6,6,F,IF,X,IX) CONST = LOG(CONST) C WRITE (8,FMT=*) 'CONST (1):',DBLE(2**J)*CONST DO 90 K = 1,J CALL F06ZAF('N','N',6,6,6,DIAG,F,IF,X,IX,OFDIAG,XOUT,IXOUT) CALL F06TFF('G',6,6,XOUT,IXOUT,F,IF) CALL F06TFF('G',6,6,XOUT,IXOUT,X,IX) CONST = CONST + CONST DUMMY = F06UAF('M',6,6,F,IF,WREAL) C WRITE (8,FMT=*) 'DUMMY=',DUMMY DO 80 IJ = 1,6 DO 70 II = 1,6 F(II,IJ) = F(II,IJ)/DUMMY X(II,IJ) = X(II,IJ)/DUMMY 70 CONTINUE 80 CONTINUE CONST = CONST + LOG(DUMMY) 90 CONTINUE C F now contains an approximation to exp(A(xend-xstart))/exp(CONST). 9000 FORMAT (A) RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETMSH(A,B,XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS, & TOL,IFAIL) IMPLICIT NONE C C This subroutine uses the natural error estimate from the Magnus formula C as a means of automatic mesh generation. C INTEGER NMESH,IMATCH,IFAIL,I DOUBLE PRECISION A,B,XMESH(0:NMESH),TOL,H,ERREST,ERRMIN DOUBLE COMPLEX ELAM,SIGMA(6,6) EXTERNAL SL4COF,SL4BCS C The direction of integration is not important for setting up the mesh, so C we do it in one big sweep from x=A to x=B. C Set the initial stepsize: IF (B.LE.A) THEN IFAIL = -1 RETURN END IF H = (B-A)/DBLE(MAX(NMESH,100)) I = 0 XMESH(I) = A 100 IF ((I+1).GT.NMESH) THEN IFAIL = 2 RETURN END IF XMESH(I+1) = XMESH(I) + H IF (XMESH(I+1).GT.B) XMESH(I+1) = B C Get error estimate CALL GETSIG(SL4COF,ELAM,XMESH(I),XMESH(I+1),SIGMA,6,ERREST) IF (ABS(ERREST).LE.TOL) THEN C The step was successful. We can go to the next step, but first we C may want to increase the stepsize. IF (2.D0*ABS(ERREST).LT.TOL) THEN H = H/MAX(0.2D0,(2.D0*ABS(ERREST)/TOL)**0.2D0) END IF I = I + 1 IF (XMESH(I).NE.B) GOTO 100 C If we are here it means we have finished meshing; re-set the value of C NMESH. NMESH = I ELSE C The step was unsuccessful. We must reduce the stepsize and repeat C the step. H = H*MIN(0.8D0,(TOL/ABS(ERREST))**0.2D0) IF (ABS(H).LT.1.D-8) THEN C The stepsize is too small to make further progress IFAIL = 1 RETURN END IF GOTO 100 END IF C We have now finished meshing. We have re-set the value of NMESH. We must C also choose the best matching point. Choose the point where there is `the C least going on'. IMATCH = 0 ERRMIN = 1.0D38 DO 10 I = 1,NMESH-1 CALL GETSIG(SL4COF,ELAM,XMESH(I),XMESH(I+1),SIGMA,6,ERREST) H = XMESH(I+1)-XMESH(I) H = H/MAX(1.D0,MAX(XMESH(I)**2,XMESH(I+1)**2)) ERREST = ERREST/(H**3) IF (ERREST.LT.ERRMIN) IMATCH = I ERRMIN = MIN(ERRMIN,ERREST) 10 CONTINUE C The matching point has now been reset. IFAIL = 0 RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETSIG(SL4COF,ELAM,XSTART,XEND,SIGMA,ISIGMA,ERREST) C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION XEND,XSTART,ERREST INTEGER ISIGMA C .. C .. Array Arguments .. DOUBLE COMPLEX SIGMA(ISIGMA,6) C .. C .. Subroutine Arguments .. EXTERNAL SL4COF C .. C .. Local Scalars .. DOUBLE PRECISION H,T1,T2 INTEGER I,J C .. C .. Local Arrays .. DOUBLE COMPLEX A1(6,6),A2(6,6),ADIFF(6,6),C1(6,6),C2(6,6) C .. C .. Parameters .. DOUBLE PRECISION Q1,Q2,W1,W2 PARAMETER (Q1=0.21132486540519D0,Q2=0.78867513459481D0, & W1=0.14433756729741D0,W2=1.D0/80.D0) C .. C .. External Subroutines .. EXTERNAL GETA,LIEB C .. T1 = XSTART + Q1* (XEND-XSTART) T2 = XSTART + Q2* (XEND-XSTART) H = XEND - XSTART CALL GETA(SL4COF,ELAM,T1,A1,6) CALL GETA(SL4COF,ELAM,T2,A2,6) CALL LIEB(A2,6,A1,6,C1,6,6) DO 20 J = 1,6 DO 10 I = 1,6 ADIFF(I,J) = A2(I,J) - A1(I,J) 10 CONTINUE 20 CONTINUE CALL LIEB(ADIFF,6,C1,6,C2,6,6) ERREST = 0.D0 DO 40 J = 1,6 DO 30 I = 1,6 SIGMA(I,J) = H* ((A1(I,J)+A2(I,J))/2.D0+ & H* (W1*C1(I,J)+H*W2*C2(I,J))) ERREST = MAX(ERREST,ABS(W2*C2(I,J))) 30 CONTINUE 40 CONTINUE H = H/MAX(1.D0,MAX(XSTART**2,XEND**2)) ERREST = ABS(ERREST*H*H*H) RETURN END C --------------------------------------------------------------------------- SUBROUTINE LIEB(A,IA,B,IB,C,IC,N) C Compute AB-BA and put the answer in C C .. Parameters .. DOUBLE COMPLEX CONE,CMONE,CZERO PARAMETER (CONE= (1.D0,0.D0),CMONE= (-1.D0,0.D0), & CZERO= (0.D0,0.D0)) C .. C .. Scalar Arguments .. INTEGER IA,IB,IC,N C .. C .. Array Arguments .. DOUBLE COMPLEX A(IA,N),B(IB,N),C(IC,N) C .. C .. External Subroutines .. EXTERNAL F06ZAF C .. CALL F06ZAF('N','N',N,N,N,CONE,A,IA,B,IB,CZERO,C,IC) CALL F06ZAF('N','N',N,N,N,CMONE,B,IB,A,IA,CONE,C,IC) RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETA(SL4COF,ELAM,X,A,IA) C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION X INTEGER IA C .. C .. Array Arguments .. DOUBLE COMPLEX A(IA,6) C .. C .. Subroutine Arguments .. EXTERNAL SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX Q0,Q1,Q2,Q3 DOUBLE PRECISION P INTEGER I,J C .. CALL SL4COF(X,ELAM,P,Q3,Q2,Q1,Q0) DO 20 I = 1,6 DO 10 J = 1,6 A(J,I) = 0.D0 10 CONTINUE 20 CONTINUE A(1,3) = 1.D0/P A(2,2) = -Q3/P A(2,3) = Q0 A(2,5) = Q1 A(2,6) = -Q2 A(3,1) = -Q2 A(3,3) = -Q3/P A(3,4) = 1.D0 A(3,5) = 1.D0 A(4,1) = -Q1 A(4,6) = 1.D0 A(5,5) = -Q3/P A(5,6) = 1.D0 A(6,1) = Q0 A(6,2) = 1.D0/P RETURN END C --------------------------------------------------------------------------- C Routine to do multiplication A times v, where v is a vector in R^{6}. SUBROUTINE AVMULT(P,Q3,Q2,Q1,Q0,SCAL,VIN,IVIN,VOUT,IVOUT) C ODE is C C (P y'')'' + (Q3 y'')' + (Q2 y')' + Q1 y' + Q0 y = 0 C C C .. Scalar Arguments .. DOUBLE COMPLEX Q0,Q1,Q2,Q3 DOUBLE PRECISION P,SCAL INTEGER IVIN,IVOUT C .. C .. Array Arguments .. DOUBLE COMPLEX VIN(IVIN),VOUT(IVOUT) C .. C The contents of this routine specify the linear system satisfied by the C Riccati-like variables for the 4th order ODE above. C The Riccati-like variables are as follows: C C z(1) = det(U), z(2) = det(V), UV^A = /z(3) z(4)\ C \z(5) z(6)/ C C The quasiderivatives are given by (y,y',-(P y'')' + Q3 y'' + Q2 y', Py''). C VOUT(1) = SCAL* (VIN(3)/P) VOUT(2) = SCAL* (- (Q3/P)*VIN(2)+Q0*VIN(3)+Q1*VIN(5)-Q2*VIN(6)) VOUT(3) = SCAL* (-Q2*VIN(1)- (Q3/P)*VIN(3)+VIN(4)+VIN(5)) VOUT(4) = SCAL* (-Q1*VIN(1)+VIN(6)) VOUT(5) = SCAL* (- (Q3/P)*VIN(5)+VIN(6)) VOUT(6) = SCAL* (Q0*VIN(1)+VIN(2)/P) RETURN END C Routine to do multiplication A times RM, where RM is a 6x6 matrix SUBROUTINE AMMULT(P,Q3,Q2,Q1,Q0,SCAL,RMIN,IRMIN,RMOUT,IRMOUT) C ODE is C C (P y'')'' + (Q3 y'')' + (Q2 y')' + Q1 y' + Q0 y = 0 C C C Do the multiplication column-by-column: C .. Scalar Arguments .. DOUBLE COMPLEX Q0,Q1,Q2,Q3 DOUBLE PRECISION P,SCAL INTEGER IRMIN,IRMOUT C .. C .. Array Arguments .. DOUBLE COMPLEX RMIN(IRMIN,6),RMOUT(IRMOUT,6) C .. C .. Local Scalars .. INTEGER I C .. C .. External Subroutines .. EXTERNAL AVMULT C .. DO 10 I = 1,6 CALL AVMULT(P,Q3,Q2,Q1,Q0,SCAL,RMIN(1,I),IRMIN,RMOUT(1,I), & IRMOUT) 10 CONTINUE RETURN END C ------------------------------------------------------------------- SUBROUTINE NEST(MNEST,IMNEST,N,SCALAR,FACTOR,IFACT,STORE,ISTOR) C This routine performs the operation C C mnest = Identity + scalar.factor.mnest, C C where scalar is a real number, factor and mnest are nxn matrices, C imnest is the first dimension of mnest as declared in the C calling program and ifact is the first dimension of factor C as declared in the calling program. C C .. Scalar Arguments .. DOUBLE COMPLEX SCALAR INTEGER IFACT,IMNEST,ISTOR,N C .. C .. Array Arguments .. DOUBLE COMPLEX FACTOR(IFACT,N),MNEST(IMNEST,N),STORE(ISTOR,N) C .. C .. External Subroutines .. EXTERNAL F06TFF,F06THF,F06ZAF C .. C .. Local Scalars .. DOUBLE COMPLEX CONE,CZERO C .. C .. Intrinsic Functions .. INTRINSIC CMPLX C .. CZERO = CMPLX(0.D0,0.D0) CONE = CMPLX(1.D0,0.D0) CALL F06TFF('G',N,N,MNEST,IMNEST,STORE,ISTOR) CALL F06THF('G',N,N,CZERO,CONE,MNEST,IMNEST) CALL F06ZAF('N','N',N,N,N,SCALAR,FACTOR,IFACT,STORE,ISTOR,CONE, & MNEST,IMNEST) RETURN END