C C --------------------------------------------------------------- C User Notes for LAPACK version of 2nd order non-selfadjoint code C --------------------------------------------------------------- C C As a user you will obviously want to know where you need to go to put in your C own problems. The answer is that the routines SL2COF (coefficient functions) C and SL2BCS (boundary conditions) are the routines which you need to edit. C C At the moment, these two routines contain a list of five problems, together C with a computed GOTO statement: at the top of each routine you will find a C statement GOTO (10,20,30,40,50), and the labels 10,20, etc. are placed next C to those bits of code which define the coefficients or BCs for problems 1,2,etc. C You can add your own problems to the lists just by extending everything: C C GOTO (10,20,30,40,50,60,70) C C 60 C ... (First User problem)... C RETURN C C 70 C ... (Second User problem)... C RETURN C C Because I have not yet written proper documentatation in LaTeX, such documentation C as exists is all in the code. If you get stuck please let me know. If *you* cannot C use the code, then there is not much hope for anyone else. You will probably want C to read down through the driver program just to try to work out what is going on. C C LAPACK C ------ C You need to install LAPACK on your system. Installing LAPACK from C www.netlib.org/lapack is very easy. You just follow the instructions carefully C and it installs itself. C C C C Finally, three words of warning. C C 1. Recursion. The code uses recursion, a very standard programming trick in which C a routine calls itself. Unfortunately it is not standard FORTRAN. All currently C available FORTRAN compilers that I know about do support recursion, but the DEC C Alpha system supports it in a crazy way where it tries to put everything on stack C and runs out of memory. The cure is to compile the routine ALLEVS separately from C all the others, then link them afterwards. C C 2. Reliability. The reliability is much improved from the early days of this code. C Nevertheless it can be broken. The most common failure arises when there is an C eigenvalue close to a boundary of a box in the complex plane in which you want to C find eigenvalues. The code cannot then tell whether the eigenvalue is inside or C outside the box, and fails. Shifting the box slightly usually cures this. C C 3. Log files. The code produces log files. These are largely for my benefit. If the C code crashes or fails and you cannot work out why, then I may ask you to send me C the log files. C C Marco Marletta mm7@mcs.le.ac.uk C C PROGRAM MAIN IMPLICIT NONE C .. Parameters .. INTEGER NZMAX,NFMAX,IMESH,IMLOC,NEV PARAMETER (NZMAX=1000,NFMAX=NZMAX,IMESH=1000,IMLOC=10,NEV=100) C .. C .. Local Scalars .. DOUBLE COMPLEX ANS,ANSBC,ANSLC,ANSRC,ANSTC, & ANSWR1,ANSWR2,DMISS1,DMISS2,EBL,EBR,ELAM, & ELAM1,ELAM2,ETL,ETR,F,ZMIN,FMIN DOUBLE PRECISION A,B,CONST,CONST1,CONST2,RELAM,TOL INTEGER I,IFAIL,ILOC,IMATCH,IREVCM,NLOC,NMESH,NZACT1,NZACT2, & NZBC,NZLC,NZRC,NZTC,NEIGS,NZ CHARACTER*1 CYN C .. C .. Local Arrays .. DOUBLE COMPLEX FARRY1(0:NFMAX,1:2),FARRY2(0:NFMAX,1:2), & FB(0:NFMAX,1:2,NEV),FBC1(0:NFMAX,1:2),FL(0:NFMAX,1:2,NEV), & FLC1(0:NFMAX,1:2),FR(0:NFMAX,1:2,NEV),FRC1(0:NFMAX,1:2), & FT(0:NFMAX,1:2,NEV),FTC1(0:NFMAX,1:2),ZARRY1(0:NZMAX), & ZARRY2(0:NZMAX),ZB(0:NZMAX,NEV),ZBC1(0:NZMAX), & ZL(0:NZMAX,NEV),ZLC1(0:NZMAX),ZR(0:NZMAX,NEV), & ZRC1(0:NZMAX),ZT(0:NZMAX,NEV),ZTC1(0:NZMAX), & FBC2(0:NFMAX,1:2),FRC2(0:NFMAX,1:2),FTC2(0:NFMAX,1:2), & FLC2(0:NFMAX,1:2),ZBC2(0:NZMAX),ZLC2(0:NZMAX), & ZRC2(0:NZMAX),ZTC2(0:NZMAX) DOUBLE COMPLEX ANSL(NEV),ANSR(NEV),ANST(NEV),ANSB(NEV), & EIGVAL(NEV),ZMAX(2) INTEGER NZL(NEV),NZR(NEV),NZT(NEV),NZB(NEV),IEIGS(NEV) DOUBLE PRECISION XLOC(0:IMLOC),XMESH(0:IMESH),ARGQHM(2) C .. C .. Scalars in Common .. DOUBLE PRECISION BETA,NU INTEGER IPROB C .. C .. External Subroutines .. EXTERNAL BCHILD,BXSIDE,LCHILD,LRCHLD,MKBOX,RCHILD,RFNBOX,SL2BCS, & SL2COF,SLMISS,SLSMBC,SLSMPC,TCHILD,TBCHLD,ALLEVS,GETMSH, & NZERO C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,CMPLX,DBLE,DIMAG,NINT C .. C .. Common blocks .. COMMON BETA,NU,IPROB C .. WRITE (6,FMT=*) 'Specify problem:' READ (5,FMT=*) IPROB IF (IPROB.LE.0.OR.IPROB.GT.5) STOP IF (IPROB.EQ.1) THEN BETA = 10.D0 B = 2.D0*ATAN(1.D0) A = -B ELSE IF (IPROB.EQ.2) THEN A = 0.D0 B = 10.D0 ELSE IF (IPROB.EQ.4) THEN A = 0.D0 B = 10.D0 ELSE IF (IPROB.EQ.5) THEN WRITE (6,FMT=*) 'Specify rotation angle theta:' READ (5,FMT=*) BETA C WRITE (6,FMT=*) 'Specify parameter NU:' C READ (5,FMT=*) NU A = 0.D0 B = 40.D0 ELSE A = 0.D0 B = 1.D0 END IF WRITE (6,FMT=*) 'Specify endpoints:' READ (5,FMT=*) A,B WRITE (6,FMT=*) 'Specify meshing tolerance:' READ (5,FMT=*) TOL NMESH = IMESH IFAIL = 0 ELAM = CMPLX(1.D0,1.D0) CALL GETMSH(A,B,XMESH,NMESH,IMATCH,ELAM,SL2COF,SL2BCS,TOL,IFAIL) IF (IFAIL.NE.0) THEN WRITE (6,FMT=*) 'Meshing failed, IFAIL =',IFAIL STOP END IF WRITE (6,FMT=*) 'Here is the mesh:' DO 31 I = 0,NMESH WRITE (6,FMT=*) 'I,XMESH(I):',I,XMESH(I) 31 CONTINUE NLOC = 10 DO 30 I = 0,NLOC XLOC(I) = (DBLE(NLOC-I)*A+DBLE(I)*B)/DBLE(NLOC) 30 CONTINUE ILOC = NINT(DBLE(NLOC)/2.D0) C WRITE (6,FMT=*) 'NMESH in driver:',NMESH C WRITE (6,FMT=*) 'Suggested IMATCH:',IMATCH C WRITE (6,FMT=*) 'Set IMATCH:' C READ (5,FMT=*) IMATCH IF (IPROB.EQ.4) THEN ILOC = 0 END IF IF (IPROB.NE.4) THEN WRITE (6,FMT=*) 'Do you want to do a half-plane count (Y/N)?' READ (5,FMT=9000) CYN ELSE CYN = 'N' END IF 9000 FORMAT(A) IF (CYN.EQ.'y'.OR.CYN.EQ.'Y') THEN WRITE (6,FMT=*) 'Counting all eigenvalues with real part less', & ' than RELAM; specify RELAM:' READ (5,FMT=*) RELAM IFAIL = 0 C C NZERO finds all eigenvalues with real part less than RELAM. In order to C use this routine you must be able to find a second eigenvalue problem C for which you know (a) that there are no eigenvalues with real part C less than or equal to RELAM; (b) the **total order** of the boundary C conditions is the same as for the original problem in which you are C interestes. For example, if your problem has boundary conditions C C y'(0) - y(0) = 0, y(1) = 0, C C then the total order of the boundary conditions is 1 + 0 = 0 (the C sums of the highest orders of the derivatives appearing). In this case C an appropriate second eigenproblem might use the simpler boundary C conditions C C y'(0) = 0, y(1) = 0, C C for which it is easier to find a problem all of whose eigenvalues are C to the right of RELAM. C C This second problem must be specified by the routines SLSMPC and SLSMBC C (coefficients and boundary conditions respectively) just as the original C problem in which you are really interested is specified by SL2COF and C SL2BCS. C CALL NZERO(XMESH,NMESH,IMATCH,RELAM,SL2COF,SL2BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQHM,IFAIL) IF (IFAIL.NE.0) THEN WRITE (6,FMT=*) 'Failure in call to NZERO, IFAIL=',IFAIL STOP END IF WRITE (6,FMT=*) 'No. of eigenvalues with real part less than', & RELAM,' is',NZ END IF WRITE (6,FMT=*) 'Do you want to find all the eigenvalues in', & ' a box (Y/N)?' READ (5,FMT=9000) CYN IF (CYN.NE.'Y'.AND.CYN.NE.'y') STOP C Now decide how many eigenvalues are in a given box. WRITE (6,FMT=*) 'Specify bottom right corner of box:' READ (5,FMT=*) EBR WRITE (6,FMT=*) 'Specify top left corner of box:' READ (5,FMT=*) ETL ETR = DCMPLX(DBLE(EBR),DIMAG(ETL)) EBL = DCMPLX(DBLE(ETL),DIMAG(EBR)) WRITE (6,FMT=*) 'Top right corner:',ETR WRITE (6,FMT=*) 'Bottom left corner:',EBL IFAIL = 0 C C The routine MKBOX integrates round a box in the complex plane with corners C EBR (bottom right), ETR (top right), ETL (top left) and EBL (bottom left), C thereby counting the eigenvalues in the box. It must be called before any C attempt to find the eigenvalues in the box. CALL MKBOX(EBR,ETR,ETL,EBL,ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT, & ZB,NZMAX,NZB,FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL, & ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH,SL2COF,SL2BCS, & IFAIL) NEIGS = NINT(DIMAG(ANS)) WRITE (6,FMT=*) 'IFAIL was',IFAIL WRITE (6,FMT=*) 'Number of eigenvalues was',NEIGS WRITE (6,FMT=*) 'ANS was',ANS WRITE (6,FMT=*) 'Numbers NZR,NZT,NZL,NZB:',NZR(1),NZT(1),NZL(1), & NZB(1) IF (NEIGS.LE.NEV) THEN WRITE (6,FMT=*) 'Computing all eigenvalues:' WRITE (6,FMT=*) 'Specify tolerance:' TOL = 1.D-7 READ (5,FMT=*) TOL IFAIL = 0 DO 77 I = 1,NEIGS EIGVAL(I) = CMPLX(0.D0,0.D0) 77 CONTINUE IEIGS(1) = NEIGS C C The routine ALLEVS can be called after a call to MKBOX to find all the C eigenvalues in the box specified. The storage arrays supplied must be C big enough for the number of eigenvalues in the box, hence the C IF (NEIGS.LE.NEV) statement above. TOL controls the rootfinding C accuracy and EIGVAL(1)...EIGVAL(NEIGS) will contain the eigenvalues C on exit. C CALL ALLEVS(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB, & FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST, & ANSB,ZLC1,NZMAX,ZRC1,NZMAX,ZTC1,NZMAX,ZBC1,NZMAX, & FLC1,NFMAX,FRC1,NFMAX,FTC1,NFMAX,FBC1,NFMAX, & ZLC2,NZMAX,ZRC2,NZMAX,ZTC2,NZMAX,ZBC2,NZMAX, & FLC2,NFMAX,FRC2,NFMAX,FTC2,NFMAX,FBC2,NFMAX, & XMESH,NMESH,IMATCH,SL2COF,SL2BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) WRITE (6,FMT=*) 'Completed problem',IPROB IF (IPROB.NE.5) THEN WRITE (6,FMT=*) 'Here are the eigenvalues:' DO 17 I = 1,NEIGS WRITE (6,FMT=*) EIGVAL(I) WRITE (16,FMT=*) EIGVAL(I) 17 CONTINUE WRITE (6,FMT=*) 'IFAIL was',IFAIL ELSE WRITE (6,FMT=*) 'Here are the resonances:' DO 19 I = 1,NEIGS WRITE (6,FMT=*) EIGVAL(I)/DCMPLX(COS(2.D0*BETA), & SIN(2.D0*BETA)) WRITE (16,FMT=*) EIGVAL(I)/DCMPLX(COS(2.D0*BETA), & SIN(2.D0*BETA)) 19 CONTINUE WRITE (6,FMT=*) 'IFAIL was',IFAIL END IF END IF STOP END C ------------------------------------------------------------------------------ SUBROUTINE SL2COF(X,ELAM,P,Q1,Q0) C C This subroutine specifies the coefficients in the differential equation. The C differential equation is C C -(P(x)y')' + Q1(x)y' + Q0(x,lambda)y = 0. C C X denotes x, ELAM denotes lambda. The user must write this routine to specify C the values of the coefficients which correspond to a given X and ELAM. C C C .. Scalar Arguments .. DOUBLE COMPLEX ELAM,Q0,Q1 DOUBLE PRECISION P,X C .. C .. Local Scalars .. DOUBLE PRECISION Q,S DOUBLE COMPLEX C,Z C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,COS,SIN C .. C .. Scalars in Common .. DOUBLE PRECISION BETA,NU INTEGER IPROB C .. C .. Common blocks .. COMMON BETA,NU,IPROB C .. GOTO (10,20,30,40,50) IPROB C The first problem is the Coffey-Evans equation with a complex shifted potential. 10 P = 1.D0 Q1 = CMPLX(0.D0,0.D0) Q = -2.D0*BETA*COS(2.D0*X) + (BETA*SIN(2.D0*X))**2 Q0 = CMPLX(Q,1.D0) - ELAM RETURN C The second problem is Brian Davies rotated complex harmonic oscillator. 20 P = 1.D0 Q1 = CMPLX(0.D0,0.D0) Q0 = CMPLX(1.D0,3.D0)*X**2 - ELAM RETURN C This is just the Fourier problem -y'' = \lambda y. 30 P = 1.D0 Q1 = CMPLX(0.D0,0.D0) Q0 = CMPLX(0.D0,0.D0) - ELAM RETURN C This is -y'' + exp(-x^2)y = \lambda y 40 P = 1.D0 Q1 = DCMPLX(0.D0,0.D0) Q0 = DCMPLX(-EXP(-X*X),0.D0) - ELAM RETURN C This is a resonance problem for Malcolm Brown 50 P = 1.D0 Q1 = DCMPLX(0.D0,0.D0) C = DCMPLX(COS(BETA),SIN(BETA)) Z = C*DCMPLX(X,0.D0) C Q0 = DCMPLX(-16.D0,0.D0)*EXP(-Z/4.D0)*COS(Z) C Q0 = -EXP(-Z*Z) Q0 = 16.D0*Z*Z*EXP(-Z) C Q0 = Z*Z*EXP(-NU*Z*Z) C Q0 = 1.D0/(CMPLX(1.D0,0.D0)+Z)**2 Q0 = C*C*Q0 - ELAM RETURN END C ----------------------------------------------------------- SUBROUTINE SL2BCS(IEND,X,ELAM,Y,PDY) C .. Scalar Arguments .. DOUBLE COMPLEX ELAM,Y,PDY DOUBLE PRECISION X C .. C .. Scalars in Common .. DOUBLE PRECISION BETA,NU INTEGER IPROB C .. C .. Common blocks .. COMMON BETA,NU,IPROB C .. C .. External Functions .. DOUBLE COMPLEX MYSQRT EXTERNAL MYSQRT C .. Local Scalars .. DOUBLE PRECISION THETA,P DOUBLE COMPLEX IOMEG1,IOMEG2,E1,E2,Q1,Q GOTO (10,20,30,40,50) IPROB C Dirichlet BCs for problem 1: 10 Y = CMPLX(0.D0,0.D0) PDY = CMPLX(1.D0,0.D0) RETURN C Dirichlet BCs for problem 2: 20 Y = CMPLX(0.D0,0.D0) PDY = CMPLX(1.D0,0.D0) RETURN C Dirichlet BCs for problem 3: 30 Y = CMPLX(0.D0,0.D0) PDY = CMPLX(1.D0,0.D0) RETURN C C Neumann BCs at left hand end, weird lambda-dependent BCs at right hand C end for problem 4: 40 IF (IEND.EQ.0) THEN Y = DCMPLX(1.D0,0.D0) PDY = DCMPLX(0.D0,0.D0) ELSE Y = DCMPLX(1.D0,0.D0) PDY = DCMPLX(0.D0,1.D0)*MYSQRT(ELAM) C WRITE (10,FMT=*) 'LAMBDA was',ELAM C WRITE (10,FMT=*) 'MYSQRT(LAMBDA) was',MYSQRT(ELAM) END IF RETURN C Various boundary conditions for problem 5, the resonance problem: 50 IF (IEND.EQ.0) THEN Y = CMPLX(0.D0,0.D0) PDY = CMPLX(1.D0,0.D0) ELSE C CALL SL2COF(X,ELAM,P,Q1,Q) C Y = CMPLX(1.D0,0.D0) C PDY = -SQRT(Q) C IF (DBLE(PDY).GT.0.D0) PDY = -PDY Y = CMPLX(0.D0,0.D0) PDY = CMPLX(1.D0,0.D0) END IF RETURN END C ------------------------------------------------------------------- DOUBLE COMPLEX FUNCTION MYSQRT(Z) DOUBLE COMPLEX Z DOUBLE PRECISION X,Y,THETA C Computes the square root whose argument lies in [Pi,2.Pi). X = DBLE(Z) Y = DIMAG(Z) C Get the argument of Z: THETA = ATAN2(Y,X) IF (THETA.LT.0.D0) THEN C Add 2.Pi to the argument of Z THETA = THETA + 8.D0*ATAN(1.D0) END IF C Arg(Z) now lies in [0,2.Pi) THETA = THETA/2.D0 MYSQRT = -((X*X+Y*Y)**0.25D0)*DCMPLX(COS(THETA),SIN(THETA)) RETURN END C ------------------------------------------------------------------ SUBROUTINE SLSMPC(X,ELAM,P,Q1,Q0) C .. Scalar Arguments .. DOUBLE COMPLEX ELAM,Q0,Q1 DOUBLE PRECISION P,X C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DIMAG C .. C C The comparison problem is -y'' + (1+Re(lambda))y = \lambda y C C NB Re(lambda) is constant on a vertical line in the complex plane. C P = 1.D0 Q1 = CMPLX(0.D0,0.D0) Q0 = CMPLX(1.D0,-DIMAG(ELAM)) RETURN END C ------------------------------------------------------------------ SUBROUTINE SLSMBC(ISING,X,ELAM,Y,PDY) C .. Scalar Arguments .. DOUBLE COMPLEX ELAM,Y,PDY DOUBLE PRECISION X C .. C .. Intrinsic Functions .. INTRINSIC CMPLX C .. C C Only Dirichlet BCs are currently implemented for the comparison problem C Y = CMPLX(0.D0,0.D0) PDY = CMPLX(1.D0,0.D0) RETURN END C ========================== END OF DRIVER ========================== C ============================= START of SLNSA2.F ========================= C Routines for 2nd order non-selfadjoint problems SUBROUTINE SLMISS(XMESH,NMESH,IMATCH,ELAM,SL2COF,SL2BCS,F,CONST, & IFLAG) IMPLICIT NONE INTEGER NMESH,IMATCH,IFLAG DOUBLE PRECISION XMESH(0:NMESH),CONST DOUBLE COMPLEX ELAM,F EXTERNAL SL2COF,SL2BCS C This routine computes the miss-distance function for the second order C equation C C -(p(x)y')' + q_1(x) y' + q_0(x) y = 0 C C Here any of the coefficients is (in principle) allowed to depend on the C eigenparameter lambda = ELAM. C C The equation is cast in the quasi-Hamiltonian form C C J z' = S z C C where J = /0 -1\ and S = /-q_0(x) -q_1(x)/p(x)\ C \1 0/ \ 0 1/p(x) / C C and z = (y,py')^T. C C For this problem the Riccati-like variables are y and py', so the C Riccati-like system is just the same as the system for z. C INTEGER IEND DOUBLE COMPLEX ZL(2),ZR(2),Y,PDY DOUBLE PRECISION CONSTL, CONSTR C Get boundary conditions at the left hand endpoint: IEND = 0 CALL SL2BCS(IEND,XMESH(0),ELAM,Y,PDY) ZL(1) = Y ZL(2) = PDY C WRITE (7,FMT=*) 'Left BCS Y, PDY:',Y,PDY C Integrate from left endpoint to centre: CALL EXPINT(ZL,CONSTL,XMESH,NMESH,0,IMATCH,SL2COF,ELAM,IFLAG) IF (IFLAG.NE.0) RETURN C WRITE (7,FMT=*) 'At matchpoint, ZL, CONST:',ZL(1),ZL(2),CONSTL C Get boundary conditions at the right hand endpoint: IEND = 1 CALL SL2BCS(IEND,XMESH(NMESH),ELAM,Y,PDY) ZR(1) = Y ZR(2) = PDY C WRITE (7,FMT=*) 'Right BCS Y, PDY:',Y,PDY CALL EXPINT(ZR,CONSTR,XMESH,NMESH,NMESH,IMATCH,SL2COF,ELAM,IFLAG) IF (IFLAG.NE.0) RETURN C WRITE (7,FMT=*) 'At matchpoint, ZR, CONST:',ZR(1),ZR(2),CONSTL C Miss distance is the determinant of the matrix C C / yl yr \ C \ pyl' pyr'/ C C F = ZL(1)*ZR(2)-ZL(2)*ZR(1) C WRITE (7,FMT=*) 'ELAM,F,CONST:',ELAM,F,CONSTL+CONSTR CONST = CONSTL+CONSTR RETURN END C ------------------------------------------------------------------- SUBROUTINE EXPINT(Z,CONST,XMESH,NMESH,ISTART,IEND,SL2COF,ELAM, & IFAIL) C C This routine integrates the linear system for the 2nd order C Riccati variables across a mesh in each interval of which the C coefficients in the 2nd order ODE are approximated using the C Magnus method (and are therefore treated as 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 SL2COF: specifies the coefficients in the ODE, which is C C -(P 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 IMPLICIT NONE C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION CONST INTEGER IEND,IFAIL,ISTART,NMESH C .. C .. Array Arguments .. DOUBLE COMPLEX Z(2) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. C .. Local Scalars .. DOUBLE COMPLEX CONE,CZERO DOUBLE PRECISION CNLOC,ZLMAX INTEGER I,IDI,II C .. C .. Local Arrays .. DOUBLE COMPLEX F(2,2),ZLOC(2) C .. C .. External Subroutines .. EXTERNAL EXPSOL,ZGEMV,F06TFM C .. C .. Intrinsic Functions .. INTRINSIC ABS,DCMPLX,LOG,MAX,MIN C .. CONE = DCMPLX(1.D0,0.D0) CZERO = DCMPLX(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 C WRITE (7,FMT=*) '--- ELAM = ',ELAM DO 30 I = ISTART + IDI,IEND,IDI CNLOC = 0.D0 IFAIL = 0 CALL EXPSOL(SL2COF,ELAM,XMESH(I-IDI),XMESH(I),CNLOC,F,2,IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 2 RETURN END IF CONST = CONST + CNLOC C WRITE (7,FMT=*) 'Here is the fundamental matrix:' C WRITE (7,FMT=*) F(1,1),F(1,2) C WRITE (7,FMT=*) F(2,1),F(2,2) CALL ZGEMV('N',2,2,CONE,F,2,Z,1,CZERO,ZLOC,1) CALL F06TFM('G',2,1,ZLOC,2,Z,2) C Normalise to get the Riccati variables to a reasonable scale ZLMAX = 0.D0 DO 10 II = 1,2 ZLMAX = MAX(ZLMAX,ABS(Z(II))) 10 CONTINUE DO 20 II = 1,2 Z(II) = Z(II)/ZLMAX 20 CONTINUE CONST = CONST + LOG(ZLMAX) C WRITE (7,FMT=*) 'X,Z(1),Z(2),CONST:',XMESH(I),Z(1),Z(2),CONST 30 CONTINUE IFAIL = 0 RETURN END C ------------------------------------------------------------------- SUBROUTINE EXPSOL(SL2COF,ELAM,XSTART,XEND,CONST,F,IF,IFAIL) C C The Wronskians in the 2nd order system satisfy a linear system of C the form z' = Az, where A is a 2x2 matrix. Here we C attempt to solve this system. C C The equation being solved is C C -(P 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 IMPLICIT NONE C .. Parameters .. INTEGER IX,IXOUT,ID,IN,LDAF,IWORK PARAMETER (IX=2,IXOUT=IX,ID=IX,IN=IX,LDAF=IX,IWORK=2*IX) C .. C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION CONST,XEND,XSTART INTEGER IF,IFAIL C .. C .. Array Arguments .. DOUBLE COMPLEX F(IF,2) C .. C .. Local Scalars .. DOUBLE COMPLEX DIAG,EQR,FACC,MIN1K,OFDIAG,TERM DOUBLE PRECISION DUMMY,FAC,ERREST,RCOND INTEGER II,IJ,IQ,J,K CHARACTER*1 EQUED C .. C .. Local Arrays .. DOUBLE COMPLEX D(ID,IX),N(IN,IX),X(IX,IX),X1(IX,IX),XOUT(IXOUT,IX) DOUBLE COMPLEX AF(LDAF,IX),WORK(IWORK) DOUBLE PRECISION WREAL(IX),R(IX),FERR(IX),BERR(IX),RWORK(IX) INTEGER IPIV(IX) DOUBLE COMPLEX D2(ID,IX),N2(IN,IX),F2(2,2) C .. C .. External Functions .. DOUBLE PRECISION ZLANGE,X02MEP EXTERNAL ZLANGE,X02MEP C .. C .. External Subroutines .. EXTERNAL ZGESVX,F06TFM,F06THM,ZGEMM,GETSIG C .. C .. Intrinsic Functions .. INTRINSIC ABS,DCMPLX,DBLE,INT,LOG,MAX C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. OFDIAG = DCMPLX(0.D0,0.D0) DIAG = DCMPLX(1.D0,0.D0) CALL F06THM('G',2,2,OFDIAG,DIAG,X,IX) CALL GETSIG(SL2COF,ELAM,XSTART,XEND,XOUT,IXOUT,ERREST) C XOUT contains the matrix A (not in sparse form). ZLANGE computes its C sup norm. C ZLANGE here gets the (infinity) norm of a matrix J = MAX(0,1+INT(LOG(ZLANGE('I',2,2,XOUT,IXOUT,WREAL)/LOG(2.D0)))) FAC = 0.5D0**J FACC = DCMPLX(FAC,0.D0) C WRITE (8,FMT=*) 'J,FAC:',J,FAC IQ = 1 EQR = DCMPLX(1.D0/6.D0,0.D0) 10 EQR = EQR/DCMPLX(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.X02MEP(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 F06THM('G',2,2,OFDIAG,DIAG,X,IX) CALL F06THM('G',2,2,OFDIAG,DIAG,N,IN) CALL F06THM('G',2,2,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 ZGEMM('N','N',2,2,2,FACC,XOUT,IXOUT,X,IX,OFDIAG,X1,IX) CALL F06TFM('G',2,2,X1,IX,X,IX) MIN1K = (-1.D0)**K DO 30 IJ = 1,2 DO 20 II = 1,2 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 C WRITE (7,FMT=*) 'Input D:' C WRITE (7,FMT=*) D(1,1),D(1,2) C WRITE (7,FMT=*) D(2,1),D(2,2) C WRITE (7,FMT=*) 'Input N:' C WRITE (7,FMT=*) N(1,1),N(1,2) C WRITE (7,FMT=*) N(2,1),N(2,2) D2(1,1) = D(1,1) D2(2,1) = D(2,1) D2(1,2) = D(1,2) D2(2,2) = D(2,2) N2(1,1) = N(1,1) N2(2,1) = N(2,1) N2(1,2) = N(1,2) N2(2,2) = N(2,2) EQUED = 'N' C CALL F04ADF(D,ID,N,IN,2,2,F,IF,WREAL,IFAIL) CALL ZGESVX('N','N',2,2,D2,ID,AF,LDAF,IPIV, & EQUED,R,WREAL,N2,IN,F2,2,RCOND, & FERR,BERR,WORK,RWORK,IFAIL) C WRITE (7,FMT=*) 'Difference between NAG and LAPACK:', C & ABS(F(1,1)-F2(1,1)) + ABS(F(1,2)-F2(1,2)) C & + ABS(F(2,1)-F2(2,1)) + ABS(F(2,2)-F2(2,2)) F(1,1) = F2(1,1) F(2,1) = F2(2,1) F(1,2) = F2(1,2) F(2,2) = F2(2,2) 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,ZLANGE('M',2,2,F,IF,WREAL)) DO 60 IJ = 1,2 DO 50 II = 1,2 F(II,IJ) = F(II,IJ)/CONST 50 CONTINUE 60 CONTINUE CALL F06TFM('G',2,2,F,IF,X,IX) CONST = LOG(CONST) C WRITE (8,FMT=*) 'CONST (1):',DBLE(2**J)*CONST DO 90 K = 1,J CALL ZGEMM('N','N',2,2,2,DIAG,F,IF,X,IX,OFDIAG,XOUT,IXOUT) CALL F06TFM('G',2,2,XOUT,IXOUT,F,IF) CALL F06TFM('G',2,2,XOUT,IXOUT,X,IX) CONST = CONST + CONST DUMMY = ZLANGE('M',2,2,F,IF,WREAL) C WRITE (8,FMT=*) 'DUMMY=',DUMMY DO 80 IJ = 1,2 DO 70 II = 1,2 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 GETSIG(SL2COF,ELAM,XSTART,XEND,SIGMA,ISIGMA,ERREST) IMPLICIT NONE C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION XEND,XSTART,ERREST INTEGER ISIGMA C .. C .. Array Arguments .. DOUBLE COMPLEX SIGMA(ISIGMA,2) C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. C .. Local Scalars .. DOUBLE PRECISION H,T1,T2 INTEGER I,J C .. C .. Local Arrays .. DOUBLE COMPLEX A1(2,2),A2(2,2),ADIFF(2,2),C1(2,2),C2(2,2) 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(SL2COF,ELAM,T1,A1,2) CALL GETA(SL2COF,ELAM,T2,A2,2) CALL LIEB(A2,2,A1,2,C1,2,2) DO 20 J = 1,2 DO 10 I = 1,2 ADIFF(I,J) = A2(I,J) - A1(I,J) 10 CONTINUE 20 CONTINUE CALL LIEB(ADIFF,2,C1,2,C2,2,2) ERREST = 0.D0 DO 40 J = 1,2 DO 30 I = 1,2 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 GETMSH(A,B,XMESH,NMESH,IMATCH,ELAM,SL2COF,SL2BCS, & 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(2,2) EXTERNAL SL2COF,SL2BCS 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(SL2COF,ELAM,XMESH(I),XMESH(I+1),SIGMA,2,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(SL2COF,ELAM,XMESH(I),XMESH(I+1),SIGMA,2,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 LIEB(A,IA,B,IB,C,IC,N) IMPLICIT NONE 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 ZGEMM C .. CALL ZGEMM('N','N',N,N,N,CONE,A,IA,B,IB,CZERO,C,IC) CALL ZGEMM('N','N',N,N,N,CMONE,B,IB,A,IA,CONE,C,IC) RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETA(SL2COF,ELAM,X,A,IA) IMPLICIT NONE C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION X INTEGER IA C .. C .. Array Arguments .. DOUBLE COMPLEX A(IA,2) C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. C .. Local Scalars .. DOUBLE COMPLEX Q0,Q1 DOUBLE PRECISION P INTEGER I,J C .. CALL SL2COF(X,ELAM,P,Q1,Q0) A(1,1) = DCMPLX(0.D0,0.D0) A(2,1) = Q0 A(1,2) = DCMPLX(1.D0/P,0.D0) A(2,2) = Q1/P RETURN END C ============================= END OF SLNSA2.F =========================== C ============================= START of SLNSAGEN.F ======================= SUBROUTINE STRSCH(XL,XR,NUMEIG,TOL,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB, & IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ZLC1,ZRC1,ZTC1,ZBC1,ZLC2, & ZRC2,ZBC2,ZTC2,FLC1,FRC1,FTC1,FBC1,FLC2, & FRC2,FBC2,FTC2,NEV,EIGS,IEIGS,IFAIL) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C This routine STRSCH (Strip Search) computes all eigenvalues in a C vertical strip specified by the x-coordinates XL, XR. C C Input Variables: XL,XR,NUMEIG C C NUMEIG is the number of eigenvalues in the strip. This is C calculated using NZERO, and is input into STRSCH. C C Output Variables: EIGS, IFAIL. C C EIGS is an array containing the eigenvalues in the strip. C C EBR,ETR,ETL, and EBL are the corners of boxes used in STRSCH. C These are passed to MKBOX. NUMLOC is the number of eigenvalues C found in the box by MKBOX. This number is passed to ALLEVS. C NUMFND is the cumulative number of eigenvalues found in the C strip. STRSCH exits when either NUMFND = NUMEIG, or the next C box would be higher than HMAX. C C All other variables in the parameter list are passed to and from C MKBOX and ALLEVS, and are used in those subroutines. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C .. Local Scalars .. INTEGER ICOUNT,IDIRECT,IFAIL,IMATCH,IFB,IFL,IFR,IFT, & IZB,IZL,IZR,IZT,NMESH,NUMEIG,NUMLOC,NUMFND DOUBLE PRECISION XL,XR,YB,YBOLD,YT,HEIGHT,WIDTH,HMAX,TOL DOUBLE COMPLEX EBL,EBR,ETL,ETR,ANS C .. Local Arrays .. DOUBLE COMPLEX FB(0:IFB,1:2,NEV),FBC1(0:IFB,1:2), & FL(0:IFL,1:2,NEV),FLC1(0:IFL,1:2), & FR(0:IFR,1:2,NEV),FRC1(0:IFR,1:2), & FT(0:IFT,1:2,NEV),FTC1(0:IFT,1:2), & ZB(0:IZB,NEV),ZBC1(0:IZB), & ZL(0:IZL,NEV),ZLC1(0:IZL),ZR(0:IZR,NEV), & ZRC1(0:IZR),ZT(0:IZT,NEV),ZTC1(0:IZT), & FBC2(0:IFB,1:2),FRC2(0:IFR,1:2),FTC2(0:IFT,1:2), & FLC2(0:IFL,1:2),ZBC2(0:IZB),ZLC2(0:IZL), & ZRC2(0:IZR),ZTC2(0:IZT),EIGS(NEV) DOUBLE COMPLEX ANSL(NEV),ANSR(NEV),ANST(NEV),ANSB(NEV) INTEGER NZL(NEV),NZR(NEV),NZT(NEV),NZB(NEV),IEIGS(NEV) DOUBLE PRECISION XMESH(0:NMESH) C C C23456789112345678921234567893123456789412345678951234567896123456789712 C C EXTERNAL ALLEVS,MKBOX,SL4BCS,SL4COF INTRINSIC CMPLX,DIMAG,NINT IF (NUMEIG.EQ.0) THEN WRITE(6,FMT=*) 'NUMEIG =' ,0 STOP ENDIF HMAX = 1.D8 WIDTH = XR - XL HEIGHT = DMAX1(2.D0*WIDTH,1.D0) YB = 0.D0 YT = HEIGHT J = 1 IDIRECT = 1 ICOUNT = 0 ANS = CMPLX(0.D0,0.D0) NUMFND = 0 30 EBR = CMPLX(XR,YB) ETR = CMPLX(XR,YT) ETL = CMPLX(XL,YT) EBL = CMPLX(XL,YB) CALL MKBOX(EBR,ETR,ETL,EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT, & ZB,IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH,SL4COF,SL4BCS, & IFAIL) IF (IFAIL.NE.0) THEN WRITE(6,FMT=*) 'For MKBOX' WRITE(6,FMT=*) 'IFAIL =' ,IFAIL STOP ENDIF ICOUNT = ICOUNT + 1 IF (IDIRECT.EQ.1) THEN YBOLD = YB YB = -YT YT = -YBOLD IDIRECT = 2 ELSE YBOLD = YB YB = -YB YT = YB + 2.D0*(YT - YBOLD) IDIRECT = 1 ENDIF NUMLOC = NINT(DIMAG(ANS)) WRITE(6,FMT=*) 'NUMLOC = ' ,NUMLOC IF (NINT(DIMAG(ANS)).EQ.0) GOTO 30 IFAIL = 0 IEIGS(1) = NUMLOC DO 35 I = 1, NUMLOC EIGS(NUMFND+I) = CMPLX(0.D0,0.D0) 35 CONTINUE CALL ALLEVS(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC1,IZL,ZRC1,IZR,ZTC1,IZT,ZBC1,IZB, & FLC1,IFL,FRC1,IFR,FTC1,IFT,FBC1,IFB, & ZLC2,IZL,ZRC2,IZR,ZTC2,IZT,ZBC2,IZB, & FLC2,IFL,FRC2,IFR,FTC2,IFT,FBC2,IFB, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGS(NUMFND+1),NEV,IEIGS,IFAIL) IF (IFAIL.NE.0) THEN WRITE(6,FMT=*) 'For ALLEVS' WRITE(6,FMT=*) 'IFAIL =' ,IFAIL STOP ENDIF WRITE(6,FMT=*) 'Eigenvalues found in current box:' DO 40 I = 1, NUMLOC WRITE(6,FMT=*) EIGS(NUMFND+I) C WRITE(100,FMT=*) EIGS(NUMFND+I) 40 CONTINUE NUMFND = NUMFND + NUMLOC IF (NUMFND.EQ.NUMEIG) THEN RETURN ELSEIF (YT.GE.HMAX) THEN WRITE(6,FMT=*) 'Height exceeds' , HMAX WRITE(6,FMT=*) 'Number of eigenvalues in the strip is' ,NUMEIG WRITE(6,FMT=*) 'Number of eigenvalues found is' ,NUMFND WRITE(6,FMT=*) 'Here are the eigenvalues found' DO 60 I = 1, NUMFND WRITE(6,FMT=*) EIGS(I) 60 CONTINUE STOP ELSE GOTO 30 ENDIF END C ------------------------------------------------------------------------- SUBROUTINE ALLEVS(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1,ZBC1,IZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) C This routine computes all the eigenvalues in a box specified by the C arrays ZL, ZR, ZT, ZB, FL, FR, FT, FB. C C The arguments are as follows. C C ZL(0:IZL,1:NEIGS): double complex array such that the elements ZL(0:NZL(1),1) C specify the points on the left hand side of the box, starting C from the top. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZL(0:IZL,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZL is liable to have been C overwritten. C C IZL: integer, the first dimension of ZL as declared in the calling (sub)program. C Unchanged on exit. C C NZL(1:NEIGS): integer array. The nodes on the left hand side of the box are in C ZL(0:NZL(1),1). C C NZL is overwritten on exit. C C ZR(0:IZR,1:NEIGS): double complex array such that the elements ZR(0:NZR(1),1) C specify the points on the right hand side of the box, starting C from the bottom. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZR(0:IZR,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZR is liable to have been C overwritten. C C IZR: integer, the first dimension of ZR as declared in the calling (sub)program. C Unchanged on exit. C C NZR(1:NEIGS): integer array. The nodes on the right hand side of the box are in C ZR(0:NZR(1),1). C C NZR is overwritten on exit. C C ZT(0:IZT,1:NEIGS): double complex array such that the elements ZT(0:NZT(1),1) C specify the points on the top side of the box, starting C from the right. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZT(0:IZT,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZT is liable to have been C overwritten. C C IZT: integer, the first dimension of ZT as declared in the calling (sub)program. C Unchanged on exit. C C NZT(1:NEIGS): integer array. The nodes on the top side of the box are in C ZT(0:NZT(1),1). C C NZT is overwritten on exit. C C C ZB(0:IZB,1:NEIGS): double complex array such that the elements ZB(0:NZB(1),1) C specify the points on the bottom side of the box, starting C from the left. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZB(0:IZB,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZB is liable to have been C overwritten. C C IZB: integer, the first dimension of ZB as declared in the calling (sub)program. C Unchanged on exit. C C NZB(1:NEIGS): integer array. The nodes on the left hand side of the box are in C ZB(0:NZB(1),1). C C NZB is overwritten on exit. C C FL(0:IFL,1:2,1:NEIGS): double complex array such that the elements FL(0:NZL(1),1:2,1) C specify the function values at points on the left hand side of the box, starting C from the top. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FL(0:IFL,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FL is liable to have been C overwritten. C C IFL: integer, the first dimension of FL as declared in the calling (sub)program. C Unchanged on exit. C C C FR(0:IFR,1:2,1:NEIGS): double complex array such that the elements FR(0:NZR(1),1:2,1) C specify the function values at points on the right hand side of the box, starting C from the bottom. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FR(0:IFR,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FR is liable to have been C overwritten. C C IFR: integer, the first dimension of FR as declared in the calling (sub)program. C Unchanged on exit. C C C FT(0:IFT,1:2,1:NEIGS): double complex array such that the elements FT(0:NZT(1),1:2,1) C specify the function values at points on the top side of the box, starting C from the right. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FT(0:IFT,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FT is liable to have been C overwritten. C C IFT: integer, the first dimension of FT as declared in the calling (sub)program. C Unchanged on exit. C C C FB(0:IFB,1:2,1:NEIGS): double complex array such that the elements FB(0:NZB(1),1:2,1) C specify the function values at points on the bottom side of the box, starting C from the left. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FT(0:IFB,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FB is liable to have been C overwritten. C C IFB: integer, the first dimension of FB as declared in the calling (sub)program. C Unchanged on exit. C C ANSL(NEIGS): double complex array, such that ANSL(1) has been assigned to C the appropriate value corresponding to ZL and FL by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANSL may be overwritten by ALLEVS. C C ANSR(NEIGS): double complex array, such that ANSL(1) has been assigned to C the appropriate value corresponding to ZR and FR61q by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANSR may be overwritten by ALLEVS. C C ANST(NEIGS): double complex array, such that ANST(1) has been assigned to C the appropriate value corresponding to ZT and FT by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANST may be overwritten by ALLEVS. C C ANSB(NEIGS): double complex array, such that ANSB(1) has been assigned to C the appropriate value corresponding to ZB and FB by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANSB may be overwritten by ALLEVS. C C ZLC1(0:IZLC1), FLC1(0:IFLC1,1:2), ZRC1(0:IZRC1), FRC1(0:IFRC1,1:2), C ZTC1(0:IZTC1), FTC1(0:IFTC1,1:2), ZBC1(0:IZBC1), FBC1(0:IFBC1,1:2), C ZLC2(0:IZLC2), FLC2(0:IFLC2,1:2), ZRC2(0:IZRC2), FRC2(0:IFRC2,1:2), C ZTC2(0:IZTC2), FTC2(0:IFTC2,1:2), ZBC2(0:IZBC2), FBC2(0:IFBC2,1:2): C double complex workspace arrays of the dimensions shown. These C need not be assigned on entry. C C IZLC1, IZRC1, IZTC1, IZBC1, IFLC1, IFRC1, IFBC1, IFTC1, C IZLC2, IZRC2, IZTC2, IZBC2, IFLC2, IFRC2, IFBC2, IFTC2: C integers specifying the array dimensions. In general these should C be large enough to ensure that C C IZLC1 > = NZL(1), IZLC2 > = NZL(1), C IZRC1 > = NZR(1), IZRC2 > = NZR(1), C IZTC1 > = NZT(1), IZTC2 > = NZT(1), C IZBC1 > = NZB(1), IZBC2 > = NZB(1), C IFLC1 > = NZL(1), IFLC2 > = NZL(1), C IFRC1 > = NZR(1), IFRC2 > = NZR(1), C IFTC1 > = NZT(1), IFTC2 > = NZT(1), C IFBC1 > = NZB(1), IFBC2 > = NZB(1). C C XMESH(0:NMESH): double precision real; the mesh used for the X-integration. C Unchanged on exit. C C NMESH: integer, number of mesh intervals for X-integration. C C IMATCH: integer, index of shooting matchpoint XMESH(IMATCH). C C SL4COF, SL4BCS: the subroutines which define the coefficients C and boundary conditions for the problem. C C TOL: double precision real. TOL > 0 must be set on input and C speficies the accuracy for the rootfinding. C C C EIGVAL(NEIGS): double complex array, *output*. On exit this array C contains the eigenvalues found. WARNING: if a box contains several C eigenvalues and the routine does not manage to find all of them, C then several elements of EIGVAL may be undefined on exit. C C NEIGS: integer. On input, NEIGS should specify the last dimensions C of various arrays (see above). NEIGS should be greater than IEIGS(1) C (see below). NEIGS is unchanged on exit. C C IEIGS(NEIGS): integer array. On entry, IEIGS(1) must specify the number C of eigenvalues in the box, while the other elements need not be C defined. On exit, the whole of the array IEIGS is liable to have C been overwritten. C C IFAIL: integer, error flag. To be set to 0 on entry. A non-zero C value on exit denotes an error. C C IMPLICIT NONE C .. Scalar Arguments .. INTEGER IZL,IZR,IZT,IZB,IFL,IFR,IFT,IFB,NMESH,IMATCH, & NEIGS,IFAIL,IZLC1,IZRC1,IZTC1,IZBC1,IZLC2,IZRC2,IZTC2,IZBC2, & IFLC1,IFLC2,IFRC1,IFRC2,IFTC1,IFTC2,IFBC1,IFBC2 DOUBLE PRECISION TOL C .. C .. Array Arguments .. INTEGER NZL(NEIGS),NZR(NEIGS),NZT(NEIGS),NZB(NEIGS),IEIGS(*) DOUBLE COMPLEX ANSL(NEIGS),ANSR(NEIGS),ANST(NEIGS),ANSB(NEIGS), & EIGVAL(NEIGS),ZL(0:IZL,NEIGS),ZR(0:IZR,NEIGS),ZT(0:IZT,NEIGS), & ZB(0:IZB,NEIGS),FL(0:IFL,2,NEIGS),FR(0:IFR,2,NEIGS), & FT(0:IFT,2,NEIGS),FB(0:IFB,2,NEIGS) DOUBLE COMPLEX ZLC1(0:IZLC1),ZRC1(0:IZRC1),ZTC1(0:IZTC1), & ZBC1(0:IZBC1),ZLC2(0:IZLC2),ZRC2(0:IZRC2),ZTC2(0:IZTC2), & ZBC2(0:IZBC2) DOUBLE COMPLEX FLC1(0:IFLC1,2),FRC1(0:IFRC1,2),FTC1(0:IFTC1,2), & FBC1(0:IFBC1,2),FLC2(0:IFLC2,2),FRC2(0:IFRC2,2),FTC2(0:IFTC2,2), & FBC2(0:IFBC2,2) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4COF,SL4BCS C .. C .. Local Scalars .. DOUBLE COMPLEX ELOC,ANSLC1,ANSLC2,ANSRC1, & ANSRC2,ANSTC1,ANSTC2,ANSBC1,ANSBC2 INTEGER I,NZLC1,NZLC2,NZRC1,NZRC2,NZTC1,NZTC2,NZBC1, & NZBC2,IFLAG1,IFLAG2,IEIGS1,IEIGS2,INDEX DOUBLE PRECISION RLENT1,RLENT2 C .. WRITE (9,FMT=*) '****** NMESH,IMATCH:',NMESH,IMATCH IF (IMATCH.LT.0.OR.IMATCH.GT.NMESH) THEN IMATCH = NMESH/2 END IF WRITE (9,FMT=*) '****** NMESH,IMATCH:',NMESH,IMATCH IF (IEIGS(1).LE.0) THEN IFAIL = MAX(1,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF ELOC = (ZL(0,1)+ZR(0,1))/2.D0 IF (ABS(ZL(0,1)-ZR(0,1)).LT.TOL*ABS(ELOC)) THEN C The diagonal of the box has length less than TOL times the C midpoint. In this case we assign a cluster of IEIGS(1) eigenvalues C all to be equal to the midpoint. DO 10 I = 1,IEIGS(1) EIGVAL(I) = ELOC 10 CONTINUE WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF C If we are here then the box has diameter greater than we C want and the box contains at least one eigenvalue. IF (IEIGS(1).EQ.1) THEN C We have found a box containing just one eigenvalue, which C we now proceed to locate. IFLAG1 = 0 CALL RFNBOX(ZL,IZL,NZL(1),ZR,IZR,NZR(1),ZT,IZT,NZT(1),ZB,IZB, & NZB(1),FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL(1), & ANSR(1),ANST(1),ANSB(1),ZLC1,IZLC1,NZLC1,ZRC1, & IZRC1,NZRC1,ZTC1,IZTC1,NZTC1,ZBC1,IZBC1,NZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ANSLC1,ANSRC1,ANSTC1,ANSBC1,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,TOL,'C',EIGVAL(1),IFLAG1) IFAIL = MAX(IFAIL,IFLAG1) IF (IFLAG1.NE.0) THEN C If the error is not too dire we return the eigenvalue approximations C anyway. IF (ABS(ZL(0,1)-ZR(0,1)).LT.SQRT(TOL)*ABS(ELOC)) THEN EIGVAL(1) = 0.5D0*(ZL(0,1)+ZR(0,1)) END IF END IF WRITE (9,FMT=*) 'Eigenvalue located,IFLAG1, IEIGS(1)=',IFLAG1, & IEIGS(1) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN C ELSE C We form appropriate child boxes and use recursion RLENT1 = MIN(ABS(DBLE(ZL(NZL(1),1))),ABS(DBLE(ZR(0,1)))) RLENT1 = MAX(1.D0,SQRT(RLENT1)) RLENT2 = MIN(ABS(DIMAG(ZT(0,1))),ABS(DIMAG(ZB(NZB(1),1)))) RLENT2 = MAX(1.D0,SQRT(RLENT2)) C IF (ABS(ZL(NZL(1),1)-ZR(0,1))/RLENT1.GT. C & ABS(ZT(0,1)-ZB(NZB(1),1))/RLENT2) THEN IF (ABS(ZL(NZL(1),1)-ZR(0,1)).GT.ABS(ZT(0,1)-ZB(NZB(1),1))) THEN C We have a wide, short box. Form the left and right children. WRITE (9,FMT=*) 'IMATCH,NMESH,IEIGS(1) before LCHILD:',IMATCH, & NMESH,IEIGS(1) IFLAG1 = 0 CALL LCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,NZLC1,ZRC1,IZRC1,NZRC1, & ZTC1,IZTC1,NZTC1,ZBC1,IZBC1,NZBC1,FLC1,IFLC1,FRC1, & IFRC1,FTC1,IFTC1,FBC1,IFBC1,ANSLC1,ANSRC1,ANSTC1, & ANSBC1,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF WRITE (9,FMT=*) 'IMATCH,NMESH,IEIGS(1) after LCHILD:',IMATCH, & NMESH,IEIGS(1) IEIGS1 = NINT(DIMAG(ANSLC1+ANSRC1+ANSTC1+ANSBC1)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'LCHILD FOUND',IEIGS1,' EIGENVALUES' IF (IEIGS1.EQ.IEIGS(1)) THEN C All the eigenvalues are in the left box. Copy the child box C arrays onto the parent box arrays and continue. IF (NZLC1.GT.IZL.OR.NZRC1.GT.IZR.OR.NZTC1.GT.IZT. & OR.NZBC1.GT.IZB) THEN IFAIL = MAX(25,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF DO 23 I = 0,NZLC1 ZL(I,1) = ZLC1(I) FL(I,1,1) = FLC1(I,1) FL(I,2,1) = FLC1(I,2) 23 CONTINUE ANSL(1) = ANSLC1 NZL(1) = NZLC1 DO 33 I = 0,NZRC1 ZR(I,1) = ZRC1(I) FR(I,1,1) = FRC1(I,1) FR(I,2,1) = FRC1(I,2) 33 CONTINUE ANSR(1) = ANSRC1 NZR(1) = NZRC1 DO 43 I = 0,NZTC1 ZT(I,1) = ZTC1(I) FT(I,1,1) = FTC1(I,1) FT(I,2,1) = FTC1(I,2) 43 CONTINUE ANST(1) = ANSTC1 NZT(1) = NZTC1 DO 53 I = 0,NZBC1 ZB(I,1) = ZBC1(I) FB(I,1,1) = FBC1(I,1) FB(I,2,1) = FBC1(I,2) 53 CONTINUE ANSB(1) = ANSBC1 NZB(1) = NZBC1 WRITE (9,FMT=*) 'Calling ALLEVS with IEIGS(1)=',IEIGS(1) CALL ALLEVS(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1, & ZBC1,IZBC1,FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1, & IFBC1,ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS with IEIGS(1)=',IEIGS(1) RETURN END IF C Form right child: WRITE (9,FMT=*) 'IEIGS(1) before RCHILD:',IEIGS(1) IFLAG1 = 0 CALL RCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC2,IZLC2,NZLC2,ZRC2,IZRC2,NZRC2, & ZTC2,IZTC2,NZTC2,ZBC2,IZBC2,NZBC2,FLC2,IFLC2,FRC2, & IFRC2,FTC2,IFTC2,FBC2,IFBC2,ANSLC2,ANSRC2,ANSTC2, & ANSBC2,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) WRITE (9,FMT=*) 'IMATCH,NMESH,IEIGS(1),IFLAG1 after RCHILD:', & IMATCH,NMESH,IEIGS(1),IFLAG1 IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IEIGS2 = NINT(DIMAG(ANSLC2+ANSRC2+ANSTC2+ANSBC2)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'RCHILD FOUND',IEIGS2,' EIGENVALUES' WRITE (9,FMT=*) '-* IEIGS1,IEIGS2,NEIGS:',IEIGS1,IEIGS2,IEIGS(1) IF (IEIGS1+IEIGS2.NE.IEIGS(1)) THEN C The eigenvalue counts do not match up. Return error flag. C Eventually we shall have to work out what to do in this C case. WRITE (9,FMT=*) '* IEIGS1,IEIGS2,NEIGS:',IEIGS1,IEIGS2, & IEIGS(1) IFAIL = MAX(24,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF ELSE C We have a long, thin box. Form the top and bottom children. C Form top child: WRITE (9,FMT=*) 'IEIGS(1) before TCHILD:',IEIGS(1) IFLAG1 = 0 CALL TCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,NZLC1,ZRC1,IZRC1,NZRC1, & ZTC1,IZTC1,NZTC1,ZBC1,IZBC1,NZBC1,FLC1,IFLC1,FRC1, & IFRC1,FTC1,IFTC1,FBC1,IFBC1,ANSLC1,ANSRC1,ANSTC1, & ANSBC1,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) WRITE (9,FMT=*) 'NMESH,IMATCH,IEIGS(1) after TCHILD:',NMESH, & IMATCH,IEIGS(1) IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IEIGS1 = NINT(DIMAG(ANSLC1+ANSRC1+ANSTC1+ANSBC1)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'TCHILD FOUND',IEIGS1,' EIGENVALUES' IF (IEIGS1.EQ.IEIGS(1)) THEN C All the eigenvalues are in the top box. Copy the child box C arrays onto the parent box arrays and continue. IF (NZLC1.GT.IZL.OR.NZRC1.GT.IZR.OR.NZTC1.GT.IZT. & OR.NZBC1.GT.IZB) THEN IFAIL = MAX(25,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF DO 20 I = 0,NZLC1 ZL(I,1) = ZLC1(I) FL(I,1,1) = FLC1(I,1) FL(I,2,1) = FLC1(I,2) 20 CONTINUE ANSL(1) = ANSLC1 NZL(1) = NZLC1 DO 30 I = 0,NZRC1 ZR(I,1) = ZRC1(I) FR(I,1,1) = FRC1(I,1) FR(I,2,1) = FRC1(I,2) 30 CONTINUE ANSR(1) = ANSRC1 NZR(1) = NZRC1 DO 40 I = 0,NZTC1 ZT(I,1) = ZTC1(I) FT(I,1,1) = FTC1(I,1) FT(I,2,1) = FTC1(I,2) 40 CONTINUE ANST(1) = ANSTC1 NZT(1) = NZTC1 DO 50 I = 0,NZBC1 ZB(I,1) = ZBC1(I) FB(I,1,1) = FBC1(I,1) FB(I,2,1) = FBC1(I,2) 50 CONTINUE ANSB(1) = ANSBC1 NZB(1) = NZBC1 WRITE (9,FMT=*) 'Current box:' WRITE (9,FMT=*) ZT(0,1),ZL(0,1),ZB(0,1),ZR(0,1) WRITE (9,FMT=*) 'Mesh sizes (T,L,B,R):',NZT(1), & NZL(1),NZB(1),NZR(1) WRITE (9,FMT=*) '* Calling ALLEVS with IEIGS(1)=',IEIGS(1) CALL ALLEVS(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1, & ZBC1,IZBC1,FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1, & IFBC1,ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS with IEIGS(1)=',IEIGS(1) RETURN END IF C Form bottom child: WRITE (9,FMT=*) 'Before BCHILD, IEIGS(1)=',IEIGS(1) IFLAG1 = 0 CALL BCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC2,IZLC2,NZLC2,ZRC2,IZRC2,NZRC2, & ZTC2,IZTC2,NZTC2,ZBC2,IZBC2,NZBC2,FLC2,IFLC2,FRC2, & IFRC2,FTC2,IFTC2,FBC2,IFBC2,ANSLC2,ANSRC2,ANSTC2, & ANSBC2,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) WRITE (9,FMT=*) 'IMATCH,IEIGS(1) after BCHILD:',IMATCH,IEIGS(1) IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF WRITE (9,FMT=*) 'IMATCH after BCHILD:',IMATCH IEIGS2 = NINT(DIMAG(ANSLC2+ANSRC2+ANSTC2+ANSBC2)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'BCHILD FOUND',IEIGS2,' EIGENVALUES' IF (IEIGS1+IEIGS2.NE.IEIGS(1)) THEN C The eigenvalue counts do not match up. Return error flag. C Eventually we shall have to work out what to do in this C case. WRITE (9,FMT=*) 'IEIGS1,IEIGS2,NEIGS:',IEIGS1,IEIGS2, & IEIGS(1) IFAIL = MAX(23,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF END IF C Have now formed the child boxes. C Copy the child arrays into appropriate parts of the parent C arrays. IF (NZLC1.GT.IZL.OR.NZRC1.GT.IZR.OR.NZTC1.GT.IZT. & OR.NZBC1.GT.IZB) THEN IFAIL = MAX(25,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IF (NZLC2.GT.IZL.OR.NZRC2.GT.IZR.OR.NZTC2.GT.IZT. & OR.NZBC2.GT.IZB) THEN IFAIL = MAX(25,IFAIL) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IF (IEIGS1.NE.0) THEN DO 21 I = 0,NZLC1 ZL(I,1) = ZLC1(I) FL(I,1,1) = FLC1(I,1) FL(I,2,1) = FLC1(I,2) 21 CONTINUE ANSL(1) = ANSLC1 NZL(1) = NZLC1 DO 31 I = 0,NZRC1 ZR(I,1) = ZRC1(I) FR(I,1,1) = FRC1(I,1) FR(I,2,1) = FRC1(I,2) 31 CONTINUE ANSR(1) = ANSRC1 NZR(1) = NZRC1 DO 41 I = 0,NZTC1 ZT(I,1) = ZTC1(I) FT(I,1,1) = FTC1(I,1) FT(I,2,1) = FTC1(I,2) 41 CONTINUE ANST(1) = ANSTC1 NZT(1) = NZTC1 DO 51 I = 0,NZBC1 ZB(I,1) = ZBC1(I) FB(I,1,1) = FBC1(I,1) FB(I,2,1) = FBC1(I,2) 51 CONTINUE ANSB(1) = ANSBC1 NZB(1) = NZBC1 END IF INDEX = IEIGS1 + 1 DO 22 I = 0,NZLC2 ZL(I,INDEX) = ZLC2(I) FL(I,1,INDEX) = FLC2(I,1) FL(I,2,INDEX) = FLC2(I,2) 22 CONTINUE ANSL(INDEX) = ANSLC2 NZL(INDEX) = NZLC2 DO 32 I = 0,NZRC2 ZR(I,INDEX) = ZRC2(I) FR(I,1,INDEX) = FRC2(I,1) FR(I,2,INDEX) = FRC2(I,2) 32 CONTINUE ANSR(INDEX) = ANSRC2 NZR(INDEX) = NZRC2 DO 42 I = 0,NZTC2 ZT(I,INDEX) = ZTC2(I) FT(I,1,INDEX) = FTC2(I,1) FT(I,2,INDEX) = FTC2(I,2) 42 CONTINUE ANST(INDEX) = ANSTC2 NZT(INDEX) = NZTC2 DO 52 I = 0,NZBC2 ZB(I,INDEX) = ZBC2(I) FB(I,1,INDEX) = FBC2(I,1) FB(I,2,INDEX) = FBC2(I,2) 52 CONTINUE ANSB(INDEX) = ANSBC2 NZB(INDEX) = NZBC2 C Very important that the next statement be executed before C any recursive calls to ALLEVS IEIGS(INDEX) = IEIGS2 IF (IEIGS1.NE.0) THEN IEIGS(1) = IEIGS1 WRITE (13,FMT=*) 'INDEX before ALLEVS:',INDEX C IFLAG2 = 0 WRITE (9,FMT=*) '*** Calling ALLEVS with IEIGS(INDEX)=', & IEIGS(INDEX) CALL ALLEVS(ZL(0,INDEX),IZL,NZL(INDEX),ZR(0,INDEX),IZR, & NZR(INDEX),ZT(0,INDEX),IZT,NZT(INDEX),ZB(0,INDEX), & IZB,NZB(INDEX),FL(0,1,INDEX),IFL,FR(0,1,INDEX),IFR, & FT(0,1,INDEX),IFT,FB(0,1,INDEX),IFB,ANSL(INDEX), & ANSR(INDEX),ANST(INDEX),ANSB(INDEX),ZLC1,IZLC1, & ZRC1,IZRC1,ZTC1,IZTC1,ZBC1,IZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL(INDEX),NEIGS,IEIGS(INDEX),IFAIL) C IF (IFLAG2.NE.0) THEN C WRITE (9,FMT=*) 'IFLAG2 was ',IFLAG2 C END IF C IFLAG1 = 0 WRITE (9,FMT=*) '** Calling ALLEVS with IEIGS(1)=',IEIGS(1) WRITE (9,FMT=*) 'Check IMATCH before in:',IMATCH CALL ALLEVS(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1, & ZBC1,IZBC1,FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1, & IFBC1,ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) C IF (IFLAG1.NE.0) THEN C WRITE (9,FMT=*) 'IFLAG1 was ',IFLAG1 C END IF C IFAIL = MAX(IFLAG1,IFLAG2) WRITE (9,FMT=*) 'Return from ALLEVS' RETURN ELSE WRITE (13,FMT=*) 'INDEX before ALLEVS:',INDEX C IFLAG2 = 0 WRITE (9,FMT=*) '*** Calling ALLEVS with IEIGS(INDEX)=', & IEIGS(INDEX) CALL ALLEVS(ZL(0,INDEX),IZL,NZL(INDEX),ZR(0,INDEX),IZR, & NZR(INDEX),ZT(0,INDEX),IZT,NZT(INDEX),ZB(0,INDEX), & IZB,NZB(INDEX),FL(0,1,INDEX),IFL,FR(0,1,INDEX),IFR, & FT(0,1,INDEX),IFT,FB(0,1,INDEX),IFB,ANSL(INDEX), & ANSR(INDEX),ANST(INDEX),ANSB(INDEX),ZLC1,IZLC1, & ZRC1,IZRC1,ZTC1,IZTC1,ZBC1,IZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL(INDEX),NEIGS,IEIGS(INDEX),IFAIL) C IF (IFLAG2.NE.0) THEN C WRITE (9,FMT=*) 'IFLAG2 was ',IFLAG2 C END IF C IFAIL = MAX(IFAIL,IFLAG2) RETURN END IF END IF END C ----------------------------------------------------------------------- SUBROUTINE STRIP(XMESH,NMESH,IMATCH,SL4COF,SL4BCS,SLSMPC,SLSMBC, & RELAM,EPS,K,TOL,RLOW,RUP,DL,DU,IFAIL) C Subroutine to get a strip containing the `K-th'eigenvalue of a C non-selfadjoint SL problem C C XMESH(0:NMESH): input. This array specifies a mesh covering C the interval on which the problem is defined. C C NMESH: the number of intervals in the mesh. C C IMATCH: the matchpoint index, between 0 and NMESH. C C SL4COF: external SUBROUTINE defining the coefficients in C the differential equation. C C SL4BCS: external SUBROUTINE defining the boundary conditions C for the problem. C C SLSMPC: external SUBROUTINE defining the CONSTANT coefficients in C a simple differential equation used for asymptotics / comparison C purposes. C C SLSMBC: external SUBROUTINE defining the boundary conditions C for the simple comparison problem. C C RELAM: initial estimate of the real part of the eigenvalue(s) C sought. Need not be at all accurate. C C EPS: estimate of error in RELAM. Thus, the routine will look C initially for eigenvalues in the strip C C RELAM - EPS < Re(Lambda) < RELAM + EPS C C but this strip will be expanded or shifted if necessary C to locate eigenvalues. The final strip in which the C eigenvalue(s) is(are) located is C C RLOW < Re(Lambda) < RUP C C (see below). C C K: the routine will attempt to ensure that there are precisely C K eigenvalues with Re(Lambda) < RLOW and precisely K+1 C with Re(Lambda) < RUP. However for some values of K this C may not be possible because of the presence of eigenvalues C of multiplicity greater than 1. In such cases, the routine C will try to narrow the strip until RUP-RLOW < TOL (see below C for description of TOL). C C TOL: In case the eigenvalue sought turns out to have real part C very close to the real part of another eigenvalue, TOL C specifies the width of the strip required around the C eigenvalue. Must be set on entry; unchanged on exit. C C RLOW, RUP: If IFAIL = 0 or IFAIL > 9 on exit, then RLOW and C RUP specify a strip RLOW < Re(Lambda) < RUP containing C at least one eigenvalue of the problem. C C DL, DU: IF IFAIL = 0 or IFAIL > 9 on exit, then C C DL = -0.5 + no. of eigenvalues with real part < RLOW C DU = -0.5 + no. of eigenvalues with real part < RUP C C and hence DU - DL is the number of eigenvalues in the C strip. C C IFAIL: the error flag. Should be set to 0 on entry. On exit: C C IFAIL = 0: successful exit. RLOW and RUP specify a strip containing C an eigenvalue and DU - DL specifies the number of eigenvalues in C this strip. C IFAIL between 1 and 7: failure in the routine NZERO. See NZERO C for the meaning of this value if IFAIL. C IFAIL = 8: the routine was able to determine RLOW but not RUP. C IFAIL = 9: the routine was able to determine RUP but not RLOW. C IFAIL = 10: the strip RLOW < Re(Lambda) < RUP contains more C than one eigenvalue but the code has not managed to narrow C its width to be less than TOL, due to an internal restriction C on the number of calls to NZERO. C IFAIL = 13: the routine located RLOW and RUP with RUP-RLOW < TOL C but the strip appears to contain several eigenvalues. C C .. Parameters .. INTEGER N,LDFJAC,LR PARAMETER (N=2,LDFJAC=N,LR=N* (N+1)/2+1) C .. C .. Scalar Arguments .. DOUBLE PRECISION DL,DU,EPS,RELAM,RLOW,RUP,TOL INTEGER IFAIL,IMATCH,K,NMESH C .. C .. Local Scalars .. DOUBLE COMPLEX FMIN,FUN,ZMIN DOUBLE PRECISION CONST,EPSFCN,FACTOR,RELOC,W1,W2,XTOL INTEGER IFLAG,ILOW,IREVCM,ISECT,IUP,ML,MODE,MU,NZ,IFEVAL LOGICAL DONELO,DONEUP,NONVAL C .. C .. Local Arrays .. DOUBLE COMPLEX ZL(2),ZMAX(2),ZU(2) DOUBLE PRECISION ARGQH(2),ARGQL(2),ARGQU(2),DIAG(N), & FJAC(LDFJAC,N),FVEC(N),QTF(N),R(LR),W(N,4),X(N) C .. C .. External Functions .. DOUBLE PRECISION X02MEP EXTERNAL X02MEP C .. C .. External Subroutines .. EXTERNAL C05NDF,NZERO,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX,DBLE,DIMAG,MAX,REAL,SQRT C .. C .. Array Arguments .. DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF,SLSMBC,SLSMPC C .. C First use bisection to find a strip a < Re(lambda) < b containing C just one eigenvalue, if this is possible. DONEUP = .FALSE. DONELO = .FALSE. IUP = 0 ILOW = 0 RELOC = RELAM IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.LT.0.D0) THEN C RELOC is less than Re(lambda) RLOW = RELOC DL = DBLE(NZ-K) - 0.5D0 ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) DONELO = .TRUE. ELSE RUP = RELOC DU = DBLE(NZ-K) - 0.5D0 ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) DONEUP = .TRUE. END IF IF (.NOT.DONEUP) THEN C Look for upper bound on eigenvalue RELOC = RELAM + EPS 10 IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.LT.0.D0) THEN C Upper bound not found; tighten lower bound RLOW = RELOC DL = DBLE(NZ-K) - 0.5D0 ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) EPS = 2.D0*EPS RELOC = RELOC + EPS IUP = IUP + 1 IF (IUP.GT.40) THEN IFAIL = 8 RETURN END IF GO TO 10 END IF C If we are here then we have found an upper bound DU = DBLE(NZ-K) - 0.5D0 RUP = RELOC ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) DONEUP = .TRUE. ELSE C IF (.NOT.DONELO) THEN C Look for lower bound on eigenvalue RELOC = RELAM - EPS 20 IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.GT.0.D0) THEN C Lower bound not found; tighten upper bound RUP = RELOC DU = DBLE(NZ-K) - 0.5D0 ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) EPS = 2.D0*EPS RELOC = RELOC - EPS ILOW = ILOW + 1 IF (ILOW.GT.40) THEN IFAIL = 9 RETURN END IF GO TO 20 END IF C If we are here then we have found a lower bound DL = DBLE(NZ-K) - 0.5D0 RLOW = RELOC ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) DONELO = .TRUE. END IF C We now have both upper and lower bounds on the eigenvalue. ISECT = 0 C NONVAL = ARGQL(1) .EQ. 0.D0 .OR. ARGQU(2) .EQ. 0.D0 C 30 IF ((NONVAL.OR.ABS(DU-DL).GT.1.5D0) .AND. C & ABS(RUP-RLOW).GT.TOL) THEN 30 IF (ABS(DU-DL).GT.1.5D0 .AND. & ABS(RUP-RLOW).GT.TOL) THEN C We want to narrow down the strip containing the eigenvalue RELOC = (ABS(DL)*RUP+ABS(DU)*RLOW)/ (ABS(DL)+ABS(DU)) IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.GT.0.D0) THEN RUP = RELOC DU = DBLE(NZ-K) - 0.5D0 ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) ELSE RLOW = RELOC DL = DBLE(NZ-K) - 0.5D0 ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) END IF C NONVAL = ARGQL(1) .EQ. 0.D0 .OR. ARGQU(2) .EQ. 0.D0 ISECT = ISECT + 1 IF (ISECT.GT.40) THEN IFAIL = 10 RETURN END IF GO TO 30 END IF C If we are here then we have either found a strip of width C less than TOL containing the eigenvalue, plus maybe some other C eigenvalues as well, or else we have found a strip containing C just one eigenvalue. IF (ABS(DU-DL).GT.1.5D0) THEN IFAIL = 13 ELSE IFAIL = 0 END IF RETURN END C ----------------------------------------------------------------------- SUBROUTINE FNDBOX(XL,XR,ICOUNT,IDIRCT,EBR,ETR,ETL, & EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB, & IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) INTEGER ICOUNT,IDIRCT,IFAIL,IMATCH,NMESH,IFB,IFL,IFR,IFT, & IZB,IZL,IZR,IZT,NZB,NZL,NZR,NZT DOUBLE PRECISION XL,XR,YB,YBOLD,YT DOUBLE PRECISION XMESH(0:NMESH) DOUBLE COMPLEX ANSB,ANSL,ANSR,ANST,ANS,EBL,EBR,ETL,ETR DOUBLE COMPLEX FB(0:IFB,1:2),FL(0:IFL,1:2),FR(0:IFR,1:2), & FT(0:IFT,1:2),ZB(0:IZB),ZL(0:IZL),ZR(0:IZR), & ZT(0:IZT) EXTERNAL MKBOX,SL4BCS,SL4COF,X01MPI DOUBLE PRECISION X01MPI INTRINSIC CMPLX,DIMAG,NINT YB = 0.D0 YT = SQRT(X01MPI(0.D0)/2.D0) IDIRCT = 1 ICOUNT = 0 ANS = CMPLX(0.D0,0.D0) 30 EBR = CMPLX(XR,YB) ETR = CMPLX(XR,YT) ETL = CMPLX(XL,YT) EBL = CMPLX(XL,YB) ICOUNT = ICOUNT + 1 CALL MKBOX(EBR,ETR,ETL,EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT, & ZB,IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH,SL4COF,SL4BCS, & IFAIL) WRITE(6,FMT=*) 'ICOUNT =' ,ICOUNT WRITE(6,FMT=*) 'ANS =' ,ANS WRITE(6,FMT=*) 'EBR =' ,EBR WRITE(6,FMT=*) 'ETR =' ,ETR WRITE(6,FMT=*) 'ETL =' ,ETL WRITE(6,FMT=*) 'EBL =' ,EBL IF (NINT(DIMAG(ANS)).EQ.0) THEN IF (IDIRCT.EQ.1) THEN YBOLD = YB YB = -YT YT = -YBOLD IDIRCT = 2 GOTO 30 ELSE YBOLD = YB YB = -YB YT = YB + 2.D0*(YT - YBOLD) IDIRCT = 1 GOTO 30 ENDIF ELSE WRITE(6,FMT=*) 'The current box has corners:' WRITE(6,FMT=*) 'EBR = ' ,EBR WRITE(6,FMT=*) 'ETR = ' ,ETR WRITE(6,FMT=*) 'ETL = ' ,ETL WRITE(6,FMT=*) 'EBL = ' ,EBL WRITE(6,FMT=*) 'The number of eigenvalues in the box is' & ,NINT(DIMAG(ANS)) ENDIF RETURN END C --------------------------------------------------------------------- SUBROUTINE RFNBOX(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,TOL,MODE,EIGVAL,IFAIL) C C This routine refines a box until one of the following is C true: C C IF (MODE.EQ.'B') THEN the code bisects the box until both sides C are of length less than TOL. C C IF (MODE.EQ.'C') THEN the code bisects the box until it appears C to be safe to use a linear approximation to the miss-distance C function to find the eigenvalue. The eigenvalue is located C to a tolerance TOL. C C C Get `length' of box C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC,EIGVAL DOUBLE PRECISION TOL INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC CHARACTER MODE C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX ELAM,F,FD1,FD2,ZNEW,ZNEW1,ZNEW2 DOUBLE PRECISION CMAX,CONST,ERROR,LENGTH INTEGER I,IEVAL LOGICAL BISECT,NONLIN,OUT C .. C .. External Subroutines .. EXTERNAL NXTBOX,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,DIMAG,EXP,MAX,MIN C .. IEVAL = 0 10 LENGTH=MAX( & ABS(ZL(0)-ZB(0))/ & MAX(1.D0,MIN(ABS(DIMAG(ZL(0))),ABS(DIMAG(ZB(0))))), & ABS(ZB(0)-ZR(0))/ & MAX(1.D0,MIN(ABS(DBLE(ZR(0))),ABS(DBLE(ZB(0)))))) C WRITE (9,FMT=*) 'Bottom right corner:',ZR(0) WRITE (9,FMT=*) 'Top right corner:',ZT(0) WRITE (9,FMT=*) 'Top left corner:',ZL(0) WRITE (9,FMT=*) 'Bottom left corner:',ZB(0) C Decide whether or not miss-distance is well-approximated C by a linear function. CALL SLMISS(XMESH,NMESH,IMATCH,ZL(0),SL4COF,SL4BCS,F,CONST,IFAIL) WRITE (7,FMT=*) 'F,FL(0,1),CONST,FL(0,2):',F,FL(0,1),CONST,FL(0,2) CALL SLMISS(XMESH,NMESH,IMATCH,ZR(0),SL4COF,SL4BCS,F,CONST,IFAIL) WRITE (7,FMT=*) 'F,FR(0,1),CONST,FR(0,2):',F,FR(0,1),CONST,FR(0,2) CALL SLMISS(XMESH,NMESH,IMATCH,ZT(0),SL4COF,SL4BCS,F,CONST,IFAIL) WRITE (7,FMT=*) 'F,FT(0,1),CONST,FT(0,2):',F,FT(0,1),CONST,FT(0,2) CALL SLMISS(XMESH,NMESH,IMATCH,ZB(0),SL4COF,SL4BCS,F,CONST,IFAIL) WRITE (7,FMT=*) 'F,FB(0,1),CONST,FB(0,2):',F,FB(0,1),CONST,FB(0,2) CMAX = MAX(DBLE(FL(0,2)),DBLE(FR(0,2))) CMAX = MAX(CMAX,DBLE(FT(0,2))) CMAX = MAX(CMAX,DBLE(FB(0,2))) FD1 = (FL(0,1)*EXP(DBLE(FL(0,2))-CMAX)- & FR(0,1)*EXP(DBLE(FR(0,2))-CMAX))/ (ZL(0)-ZR(0)) FD2 = (FT(0,1)*EXP(DBLE(FT(0,2))-CMAX)- & FB(0,1)*EXP(DBLE(FB(0,2))-CMAX))/ (ZT(0)-ZB(0)) ZNEW1 = ZL(0) - FL(0,1)*EXP(DBLE(FL(0,2))-CMAX)/FD1 ZNEW2 = ZT(0) - FT(0,1)*EXP(DBLE(FT(0,2))-CMAX)/FD2 NONLIN = ABS(FD1-FD2) .GT. 0.1D0*MIN(ABS(FD1),ABS(FD2)) WRITE (9,FMT=*) 'Initial NONLIN:',NONLIN OUT = MAX(DBLE(ZNEW1),DBLE(ZNEW2)) .GT. DBLE(ZR(0)) + TOL OUT = OUT .OR. MIN(DBLE(ZNEW1),DBLE(ZNEW2)) .LT. DBLE(ZL(0)) - TOL OUT = OUT .OR. MAX(DIMAG(ZNEW1),DIMAG(ZNEW2)) .GT. & DIMAG(ZT(0)) + TOL OUT = OUT .OR. MIN(DIMAG(ZNEW1),DIMAG(ZNEW2)) .LT. & DIMAG(ZB(0)) - TOL WRITE (9,FMT=*) 'OUT:',OUT WRITE (7,FMT=*) 'ZR(0),ZT(0),ZL(0),ZB(0):',ZR(0),ZT(0),ZL(0),ZB(0) WRITE (7,FMT=*) 'ZR(NZR),ZT(NZT),ZL(NZL),ZB(NZB):',ZR(NZR), & ZT(NZT),ZL(NZL),ZB(NZB) NONLIN = NONLIN .OR. OUT WRITE (9,FMT=*) 'FD1,FD2:',FD1,FD2 WRITE (9,FMT=*) 'ABS(FD1-FD2),MIN(ABS(FD1),ABS(FD2)):', & ABS(FD1-FD2),MIN(ABS(FD1),ABS(FD2)) WRITE (9,FMT=*) 'ZNEW1,ZNEW2:',ZNEW1,ZNEW2 C NONLIN = .TRUE. BISECT = (LENGTH.GT.TOL .AND. MODE.EQ.'B') .OR. & (NONLIN .AND. MODE.EQ.'C' .AND. LENGTH.GT.TOL) WRITE (9,FMT=*) 'MODE,NONLIN,BISECT:',MODE,NONLIN,BISECT IF (.NOT.BISECT) THEN GO TO 70 END IF C We do the bisection of the box. 20 CALL NXTBOX(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL,IFL,FR, & IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB,ZLC,IZLC,NZLC, & ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC,IZBC,NZBC,FLC,IFLC, & FRC,IFRC,FTC,IFTC,FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN C Copy new box onto old box NZL = NZLC NZR = NZRC NZT = NZTC NZB = NZBC ANSL = ANSLC ANSR = ANSRC ANST = ANSTC ANSB = ANSBC DO 30 I = 0,NZL ZL(I) = ZLC(I) FL(I,1) = FLC(I,1) FL(I,2) = FLC(I,2) 30 CONTINUE DO 40 I = 0,NZR ZR(I) = ZRC(I) FR(I,1) = FRC(I,1) FR(I,2) = FRC(I,2) 40 CONTINUE DO 50 I = 0,NZT ZT(I) = ZTC(I) FT(I,1) = FTC(I,1) FT(I,2) = FTC(I,2) 50 CONTINUE DO 60 I = 0,NZB ZB(I) = ZBC(I) FB(I,1) = FBC(I,1) FB(I,2) = FBC(I,2) 60 CONTINUE C This trap deals with the case where bisection has C been re-entered after a failed attempt to use Newton's C method. IF (IEVAL.GT.5) THEN BISECT = .TRUE. IEVAL = IEVAL - 1 GO TO 20 END IF IEVAL = 0 GO TO 10 70 CONTINUE C If we are here then we have finished bisection. IF (MODE.EQ.'B' .OR. LENGTH.LT.TOL) THEN EIGVAL = (ZL(0)+ZR(0)+ZT(0)+ZB(0))/4.D0 IFAIL = 0 IF (MODE.EQ.'C' .AND. LENGTH.LT.TOL) THEN WRITE (9,FMT=*) 'Bisected to convergence' END IF RETURN END IF C If we are here then we should find the eigenvalue by Newton's C method. IEVAL = 0 ELAM = 0.5D0* (ZNEW1+ZNEW2) 80 CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST,IFAIL) IF (IFAIL.NE.0) RETURN ZNEW = ELAM - 2.D0*F*EXP(CONST-CMAX)/ (FD1+FD2) WRITE (9,FMT=*) 'Best approx so far:',ZNEW WRITE (9,FMT=*) 'IEVAL=',IEVAL IF (ABS(ZNEW-ELAM).GT.TOL .AND. IEVAL.LT.20) THEN IEVAL = IEVAL + 1 ELAM = ZNEW IF (DBLE(ELAM).GT.DBLE(ZR(0))+TOL .OR. & DBLE(ELAM).LT.DBLE(ZL(0))-TOL .OR. & DIMAG(ELAM).GT.DIMAG(ZT(0))+TOL .OR. & DIMAG(ELAM).LT.DIMAG(ZB(0))-TOL) THEN BISECT = .TRUE. GO TO 20 END IF GO TO 80 END IF ERROR = ABS(DBLE(ZNEW-ELAM))/ & MAX(1.D0,MIN(ABS(DBLE(ELAM)),ABS(DBLE(ZNEW)))) & + ABS(DIMAG(ZNEW-ELAM))/ & MAX(1.D0,MIN(ABS(DIMAG(ELAM)),ABS(DIMAG(ZNEW)))) IF (ERROR.GT.TOL .AND. IEVAL.GE.5) THEN C EIGVAL = ZNEW C IFAIL = 22 BISECT = .TRUE. IEVAL = IEVAL - 1 WRITE (9,FMT=*) 'Return to bisection' WRITE (8,FMT=*) 'Return to bisection' WRITE (8,FMT=*) 'NZL,NZR,NZT,NZB:',NZL,NZR,NZT,NZB WRITE (8,FMT=*) 'ZL(NZL),ZB(0):',ZL(NZL),ZB(0) WRITE (8,FMT=*) 'ZB(NZB),ZR(0):',ZB(NZB),ZR(0) WRITE (8,FMT=*) 'ZR(NZR),ZT(0):',ZR(NZR),ZT(0) WRITE (8,FMT=*) 'ZT(NZT),ZL(0):',ZT(NZT),ZL(0) GO TO 20 C RETURN END IF C If we are here then we have found an eigenvalue EIGVAL = ZNEW IFAIL = 0 RETURN END C --------------------------------------------------------------------- SUBROUTINE NXTBOX(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box containing one eigenvalue, this routine bisects the box C in one of two directions, and forms a new box containing one C eigenvalue. C C First decide which way to bisect C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. INTEGER IEIGS DOUBLE PRECISION RLENT1,RLENT2 C .. C .. External Subroutines .. EXTERNAL BCHILD,LCHILD,LRCHLD,RCHILD,TCHILD,TBCHLD C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,DIMAG,NINT C .. RLENT1 = MIN(ABS(DBLE(ZL(NZL))),ABS(DBLE(ZR(0)))) RLENT1 = MAX(1.D0,SQRT(RLENT1)) RLENT2 = MIN(ABS(DIMAG(ZT(0))),ABS(DIMAG(ZB(NZB)))) RLENT2 = MAX(1.D0,SQRT(RLENT2)) IF (ABS(ZL(NZL)-ZR(0))/RLENT1.GT. & ABS(ZT(0)-ZB(NZB))/RLENT2) THEN C We have a wide, short box. Form the left and right children. IFAIL = 0 CALL LCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB,ZLC, & IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC,IZBC, & NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC,ANSLC, & ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(DIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.GT.0) THEN C The new box is the left child. GO TO 10 ELSE C The new box is the right child IFAIL = 0 CALL RCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(DIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN C Try doing an LRCHLD: IFAIL = 0 CALL LRCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(DIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN IFAIL = 20 RETURN END IF END IF END IF ELSE C We have a long, thin box. Form the top and bottom children. IFAIL = 0 CALL TCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB,ZLC, & IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC,IZBC, & NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC,ANSLC, & ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(DIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.GT.0) THEN C The new box is the left child. GO TO 10 ELSE C The new box is the right child IFAIL = 0 CALL BCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(DIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN C Try doing a TBCHLD: IFAIL = 0 CALL TBCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(DIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN IFAIL = 21 RETURN END IF END IF END IF END IF 10 IFAIL = 0 RETURN END C --------------------------------------------------------------------- SUBROUTINE LCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the left half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the left-half child, the left-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBR,ELAM,ETR,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DBLE,LOG,MIN C .. IFAIL = 0 IF (MIN(IZLC,IFLC).LT.NZL) THEN IFAIL = 18 RETURN END IF NZLC = NZL DO 10 I = 0,NZL ZLC(I) = ZL(I) FLC(I,1) = FL(I,1) FLC(I,2) = FL(I,2) ANSLC = ANSL 10 CONTINUE C That does the left. Now do the right, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get top and bottom endpoints of the right side of the child box: C these are, respectively, the midpoints of the top and bottom C sides of the parent box. ETR = 0.5D0* (ZT(0)+ZT(NZT)) EBR = 0.5D0* (ZB(0)+ZB(NZB)) C Do the integral from EBR to ETR: 20 CALL BXSIDE(EBR,ETR,ELAM,F,CONST,IREVCM,ZRC,IZRC,NZRC,FRC,IFRC, & ANSRC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done the sides. We must do the top and bottom. C Do bottom. I = 0 30 IF (DBLE(ZB(I)).LT.DBLE(EBR)) THEN ZBC(I) = ZB(I) FBC(I,1) = FB(I,1) FBC(I,2) = FB(I,2) I = I + 1 NZBC = I GO TO 30 END IF ZBC(NZBC) = ZRC(0) FBC(NZBC,1) = FRC(0,1) FBC(NZBC,2) = FRC(0,2) ANSBC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZBC ANSBC = ANSBC + LOG(FBC(I,1)/FBC(I-1,1)) 40 CONTINUE C Do top. I = 0 50 IF (DBLE(ZT(I)).GE.DBLE(ETR)) THEN I = I + 1 GO TO 50 END IF ZTC(0) = ZRC(NZRC) FTC(0,1) = FRC(NZRC,1) FTC(0,2) = FRC(NZRC,2) NZTC = NZT - I + 1 DO 60 J = 1,NZTC IND = J + I - 1 ZTC(J) = ZT(IND) FTC(J,1) = FT(IND,1) FTC(J,2) = FT(IND,2) 60 CONTINUE ANSTC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZTC ANSTC = ANSTC + LOG(FTC(I,1)/FTC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE RCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the right half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the right-half child, the right-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBL,ELAM,ETL,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DBLE,LOG,MIN C .. IFAIL = 0 IF (MIN(IZRC,IFRC).LT.NZR) THEN IFAIL = 18 RETURN END IF NZRC = NZR DO 10 I = 0,NZR ZRC(I) = ZR(I) FRC(I,1) = FR(I,1) FRC(I,2) = FR(I,2) ANSRC = ANSR 10 CONTINUE C That does the right. Now do the left, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get top and bottom endpoints of the left side of the child box: C these are, respectively, the midpoints of the top and bottom C sides of the parent box. ETL = 0.5D0* (ZT(0)+ZT(NZT)) EBL = 0.5D0* (ZB(0)+ZB(NZB)) C Do the integral from ETL to EBL: 20 CALL BXSIDE(ETL,EBL,ELAM,F,CONST,IREVCM,ZLC,IZLC,NZLC,FLC,IFLC, & ANSLC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done the sides. We must do the top and bottom. C Do top. I = 0 30 IF (DBLE(ZT(I)).GT.DBLE(ETL)) THEN ZTC(I) = ZT(I) FTC(I,1) = FT(I,1) FTC(I,2) = FT(I,2) I = I + 1 NZTC = I GO TO 30 END IF ZTC(NZTC) = ZLC(0) FTC(NZTC,1) = FLC(0,1) FTC(NZTC,2) = FLC(0,2) ANSTC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZTC ANSTC = ANSTC + LOG(FTC(I,1)/FTC(I-1,1)) 40 CONTINUE C Do bottom. I = 0 50 IF (DBLE(ZB(I)).LE.DBLE(EBL)) THEN I = I + 1 GO TO 50 END IF ZBC(0) = ZLC(NZLC) FBC(0,1) = FLC(NZLC,1) FBC(0,2) = FLC(NZLC,2) NZBC = NZB - I + 1 DO 60 J = 1,NZBC IND = J + I - 1 ZBC(J) = ZB(IND) FBC(J,1) = FB(IND,1) FBC(J,2) = FB(IND,2) 60 CONTINUE ANSBC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZBC ANSBC = ANSBC + LOG(FBC(I,1)/FBC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE LRCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the middle vertical half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBR,ELAM,ETR,F,EBL,ETL DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DBLE,LOG,MIN C .. IFAIL = 0 IREVCM = 1 IFLAG = 0 C C Get top and bottom endpoints of the right side of the child box: C ETR = 0.75D0*ZT(0)+0.25D0*ZT(NZT) EBR = 0.25D0*ZB(0)+0.75D0*ZB(NZB) C Do the integral from EBR to ETR: 20 CALL BXSIDE(EBR,ETR,ELAM,F,CONST,IREVCM,ZRC,IZRC,NZRC,FRC,IFRC, & ANSRC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done the right side. Now do the left side. ETL = 0.25D0*ZT(0)+0.75D0*ZT(NZT) EBL = 0.75D0*ZB(0)+0.25D0*ZB(NZB) C Do the integral from ETL to EBL: 21 CALL BXSIDE(ETL,EBL,ELAM,F,CONST,IREVCM,ZLC,IZLC,NZLC,FLC,IFLC, & ANSLC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 21 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done the sides. We must do the top and bottom. C Do bottom. I = 0 ZBC(I) = ZLC(NZLC) FBC(I,1) = FLC(NZLC,1) FBC(I,2) = FLC(NZLC,2) I = 1 30 IF (DBLE(ZB(I)).LE.DBLE(ZBC(0))) THEN I = I + 1 GOTO 30 END IF J = 1 31 IF (DBLE(ZB(I)).LT.DBLE(ZRC(0))) THEN ZBC(J) = ZB(I) FBC(J,1) = FB(I,1) FBC(J,2) = FB(I,2) I = I + 1 J = J + 1 NZBC = J GO TO 31 END IF ZBC(NZBC) = ZRC(0) FBC(NZBC,1) = FRC(0,1) FBC(NZBC,2) = FRC(0,2) ANSBC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZBC ANSBC = ANSBC + LOG(FBC(I,1)/FBC(I-1,1)) 40 CONTINUE C Do top. I = 0 ZTC(I) = ZRC(NZRC) FTC(I,1) = FRC(NZRC,1) FTC(I,2) = FRC(NZRC,2) I = 1 50 IF (DBLE(ZT(I)).GE.DBLE(ZRC(NZRC))) THEN I = I + 1 GO TO 50 END IF J = 1 51 IF (DBLE(ZT(I)).GT.DBLE(ZLC(0))) THEN ZTC(J) = ZT(I) FTC(J,1) = FT(I,1) FTC(J,2) = FT(I,2) I = I + 1 J = J + 1 NZTC = J GOTO 51 END IF ZTC(NZTC) = ZLC(0) FTC(NZTC,1) = FLC(0,1) FTC(NZTC,2) = FLC(0,2) ANSTC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZTC ANSTC = ANSTC + LOG(FTC(I,1)/FTC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE BCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the bottom half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the bottom-half child, the bottom-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX ELAM,ETL,ETR,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DIMAG,LOG,MIN C .. IFAIL = 0 IF (MIN(IZBC,IFBC).LT.NZB) THEN IFAIL = 18 RETURN END IF NZBC = NZB DO 10 I = 0,NZB ZBC(I) = ZB(I) FBC(I,1) = FB(I,1) FBC(I,2) = FB(I,2) ANSBC = ANSB 10 CONTINUE C That does the bottom. Now do the top, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the top of the child box: C these are, respectively, the midpoints of the left and right C sides of the parent box. ETL = 0.5D0* (ZL(0)+ZL(NZL)) ETR = 0.5D0* (ZR(0)+ZR(NZR)) C Do the integral from ETR to ETL: 20 CALL BXSIDE(ETR,ETL,ELAM,F,CONST,IREVCM,ZTC,IZTC,NZTC,FTC,IFTC, & ANSTC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done top and bottom. We must do the sides. C Do right side. I = 0 30 IF (DIMAG(ZR(I)).LT.DIMAG(ETR)) THEN ZRC(I) = ZR(I) FRC(I,1) = FR(I,1) FRC(I,2) = FR(I,2) I = I + 1 NZRC = I GO TO 30 END IF ZRC(NZRC) = ZTC(0) FRC(NZRC,1) = FTC(0,1) FRC(NZRC,2) = FTC(0,2) ANSRC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZRC ANSRC = ANSRC + LOG(FRC(I,1)/FRC(I-1,1)) 40 CONTINUE C Do left side. I = 0 50 IF (DIMAG(ZL(I)).GE.DIMAG(ETL)) THEN I = I + 1 GO TO 50 END IF ZLC(0) = ZTC(NZTC) FLC(0,1) = FTC(NZTC,1) FLC(0,2) = FTC(NZTC,2) NZLC = NZL - I + 1 DO 60 J = 1,NZLC IND = J + I - 1 ZLC(J) = ZL(IND) FLC(J,1) = FL(IND,1) FLC(J,2) = FL(IND,2) 60 CONTINUE ANSLC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZLC ANSLC = ANSLC + LOG(FLC(I,1)/FLC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE TCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the top half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the top-half child, the top-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBL,EBR,ELAM,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DIMAG,LOG,MIN C .. IFAIL = 0 IF (MIN(IZTC,IFTC).LT.NZT) THEN IFAIL = 18 RETURN END IF NZTC = NZT DO 10 I = 0,NZT ZTC(I) = ZT(I) FTC(I,1) = FT(I,1) FTC(I,2) = FT(I,2) ANSTC = ANST 10 CONTINUE C That does the top. Now do the bottom, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the bottom of the child box: C these are, respectively, the midpoints of the left and right C sides of the parent box. EBL = 0.5D0* (ZL(0)+ZL(NZL)) EBR = 0.5D0* (ZR(0)+ZR(NZR)) C Do the integral from EBL to EBR: 20 CALL BXSIDE(EBL,EBR,ELAM,F,CONST,IREVCM,ZBC,IZBC,NZBC,FBC,IFBC, & ANSBC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done top and bottom. We must do the sides. C Do left side. I = 0 30 IF (DIMAG(ZL(I)).GT.DIMAG(EBL)) THEN ZLC(I) = ZL(I) FLC(I,1) = FL(I,1) FLC(I,2) = FL(I,2) I = I + 1 NZLC = I GO TO 30 END IF ZLC(NZLC) = ZBC(0) FLC(NZLC,1) = FBC(0,1) FLC(NZLC,2) = FBC(0,2) ANSLC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZLC ANSLC = ANSLC + LOG(FLC(I,1)/FLC(I-1,1)) 40 CONTINUE C Do right side. I = 0 50 IF (DIMAG(ZR(I)).LE.DIMAG(EBR)) THEN I = I + 1 GO TO 50 END IF ZRC(0) = ZBC(NZBC) FRC(0,1) = FBC(NZBC,1) FRC(0,2) = FBC(NZBC,2) NZRC = NZR - I + 1 DO 60 J = 1,NZRC IND = J + I - 1 ZRC(J) = ZR(IND) FRC(J,1) = FR(IND,1) FRC(J,2) = FR(IND,2) 60 CONTINUE ANSRC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZRC ANSRC = ANSRC + LOG(FRC(I,1)/FRC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C ----------------------------------------------------------- SUBROUTINE TBCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the middle horizontal half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBL,EBR,ELAM,F,ETL,ETR DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DIMAG,LOG,MIN C .. IFAIL = 0 C C Do the bottom: C IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the bottom of the child box: EBL = 0.25D0*ZL(0)+0.75D0*ZL(NZL) EBR = 0.75D0*ZR(0)+0.25D0*ZR(NZR) C Do the integral from EBL to EBR: 20 CALL BXSIDE(EBL,EBR,ELAM,F,CONST,IREVCM,ZBC,IZBC,NZBC,FBC,IFBC, & ANSBC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C C Do the top: C IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the top of the child box: ETL = 0.75D0*ZL(0)+0.25D0*ZL(NZL) ETR = 0.25D0*ZR(0)+0.75D0*ZR(NZR) C Do the integral from ETR to ETL: 21 CALL BXSIDE(ETR,ETL,ELAM,F,CONST,IREVCM,ZTC,IZTC,NZTC,FTC,IFTC, & ANSTC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 21 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done top and bottom. We must do the sides. C Do left side. ZLC(0) = ZTC(NZTC) FLC(0,1) = FTC(NZTC,1) FLC(0,2) = FTC(NZTC,2) I = 1 30 IF (DIMAG(ZL(I)).GE.DIMAG(ZTC(NZTC))) THEN I = I + 1 GOTO 30 END IF J = 1 31 IF (DIMAG(ZL(I)).GT.DIMAG(ZBC(0))) THEN ZLC(J) = ZL(I) FLC(J,1) = FL(I,1) FLC(J,2) = FL(I,2) I = I + 1 J = J + 1 GO TO 31 END IF NZLC = J ZLC(NZLC) = ZBC(0) FLC(NZLC,1) = FBC(0,1) FLC(NZLC,2) = FBC(0,2) ANSLC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZLC ANSLC = ANSLC + LOG(FLC(I,1)/FLC(I-1,1)) 40 CONTINUE C Do right side. ZRC(0) = ZBC(NZBC) FRC(0,1) = FBC(NZBC,1) FRC(0,2) = FBC(NZBC,2) I = 1 50 IF (DIMAG(ZR(I)).LE.DIMAG(ZBC(NZBC))) THEN I = I + 1 GO TO 50 END IF J = 1 51 IF (DIMAG(ZR(I)).LT.DIMAG(ZTC(0))) THEN ZRC(J) = ZR(I) FRC(J,1) = FR(I,1) FRC(J,2) = FR(I,2) I = I + 1 J = J + 1 GO TO 51 END IF NZRC = J ZRC(NZRC) = ZTC(0) FRC(NZRC,1) = FTC(0,1) FRC(NZRC,2) = FTC(0,2) ANSRC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZRC ANSRC = ANSRC + LOG(FRC(I,1)/FRC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C ----------------------------------------------------------- SUBROUTINE MKBOX(EBR,ETR,ETL,EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT, & ZB,IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR, & ANST,ANSB,ANSWER,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,IFAIL) C This routine creates a box in the complex plane with corners at C EBR (bottom right), ETR (top right), ETL (top left) and EBL (bottom C left). It also computes the integrals of f'/f round the sides of the C box (right, top, left and bottom sides) taken in anti-clockwise sense C round the boundary. The integrals are in ANSR (right side), ANST (top), C ANSL (left side) and ANSB (bottom). The meshpoints used for these C integrals are (respectively) in C C ZR(I), I = 0,..,NZR; ZT(I), I = 0,..,NZT; C ZL(I), I = 0,..,NZL; ZB(I), I = 0,..,NZB. C C Note that ZR(0) = EBR, ZR(NZR) = ETR = ZT(0), C ZT(NZT) = ETL = ZL(0), ZL(NZL) = EBL = ZB(0), C ZB(NZB) = EBR = ZR(0). C C The function values are stored in C C FR(I,1), I = 0,..,NZR; FT(I,1), I = 0,..,NZT; C FL(I,1), I = 0,..,NZL; FB(I,1), I = 0,..,NZB. C C and the scaling factors are stored in C C FR(I,2), I = 0,..,NZR; FT(I,2), I = 0,..,NZT; C FL(I,2), I = 0,..,NZL; FB(I,2), I = 0,..,NZB. C C (the value of f(z(i)) is F(i,1)*exp(F(i,2))). C C C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSL,ANSR,ANST,ANSWER,EBL,EBR,ETL,ETR INTEGER IFAIL,IFB,IFL,IFR,IFT,IMATCH,IZB,IZL,IZR,IZT,NMESH,NZB, & NZL,NZR,NZT C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FL(0:IFL,1:2),FR(0:IFR,1:2), & FT(0:IFT,1:2),ZB(0:IZB),ZL(0:IZL),ZR(0:IZR), & ZT(0:IZT) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX ELAM,F DOUBLE PRECISION CONST INTEGER IFLAG,IREVCM C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC ATAN C .. C WRITE (6,FMT=*) 'Successfully called MKBOX:' C Do right side: IREVCM = 1 IFLAG = 0 C WRITE (6,FMT=*) 'Doing right side:' 10 CALL BXSIDE(EBR,ETR,ELAM,F,CONST,IREVCM,ZR,IZR,NZR,FR,IFR,ANSR, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) C WRITE (6,FMT=*) 'Returned from call to SLMISS' IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 10 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done right side. C Do top side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing top:' 20 CALL BXSIDE(ETR,ETL,ELAM,F,CONST,IREVCM,ZT,IZT,NZT,FT,IFT,ANST, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done top side. C Do left side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing left side:' 30 CALL BXSIDE(ETL,EBL,ELAM,F,CONST,IREVCM,ZL,IZL,NZL,FL,IFL,ANSL, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 30 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done left side C Do bottom side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing bottom:' 40 CALL BXSIDE(EBL,EBR,ELAM,F,CONST,IREVCM,ZB,IZB,NZB,FB,IFB,ANSB, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 40 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Have now done all sides ANSWER = (ANSL+ANSR+ANST+ANSB)/ (8.D0*ATAN(1.D0)) IFAIL = 0 RETURN END C ----------------------------------------------------------- SUBROUTINE BXSIDE(ELAM1,ELAM2,ELAM,F,CONST,IREVCM,ZARRAY,NZMAX, & NZACT,FARRAY,NFMAX,ANSWER,IFAIL) C C This routine performs the integral of f'(lambda)/f(lambda) C over the line segment from ELAM1 to ELAM2. Generally this C line segment will be parallel to either the X-axis or the C Y-axis in the complex lambda-plane. Function values are C obtained by reverse communication, for flexibility. C C Documentation: C C ELAM1, ELAM2: The endpoints of the line segment. Must be C set on input. Unchanged on output. C C ELAM: output. On intermediate return with IREVCM.NE.0, ELAM C specifies the point on the line segment at which a value C of f is required. This value of f must be assigned to the C variable F before calling BXSIDE again with all other C parameters UNCHANGED. C C F: input. On first call with IREVCM.EQ.1 F need not be assigned. C On subsequent calls with IREVCM greater than 1, F must be C assigned the value f(ELAM)*exp(-CONST). The value of CONST C (see below) should be chosen by the user to avoid overflow C in the evaluation of F. C C CONST: input. On first call with IREVCM.EQ.1 CONST need not be C assigned. On subsequent calls with IREVCM greater than 1, C CONST should be such that f(ELAM) = F*exp(CONST). C C IREVCM: input/output. On first call, IREVCM must be assigned the C value 1. Thereafter, the user must not change the value of C IREVCM. On an intermediate exit with IREVCM > 1 the user C must assign F and CONST as described above and call BXSIDE C again. On exit with IREVCM = 0, the routine has finished the C integration. C C ZARRAY(0:NZMAX): output. On final exit with IREVCM.EQ.0 this array C stores all the points at which function evaluations were made, C with ZARRAY(0) = ELAM1 and ZARRAY(NZACT) = ELAM2. C C NMAX: input. The dimension of ZARRAY as declared in the calling C (sub)program. C C NZACT: output. On final exit, NZACT specifies the number of mesh C intervals used to perform the integration. C C FARRAY(0:NFMAX,1:2): output. On final exit with IREVCM.EQ.0 the C elements of FARRAY are assigned as follows: C C FARRAY(I,1) = f(ZARRAY(I))*exp(-CONST(ZARRAY(I))) C FARRAY(I,2) = CONST(ZARRAY(I)) C C for I = 0,1,..,NZACT. C C NFMAX: input. The first dimension of FARRAY as declared in the C calling (sub)program. C C ANSWER: output. On final exit with IFAIL.EQ.0, ANSWER specifies C the value of the integral. C C IFAIL: input/output. On entry IFAIL should be assigned the value C 0. On output IFAIL.EQ.0 unless a failure has occurred. The C three nonzero values of IFAIL which may be returned are: C C IFAIL = 15: cannot manage to reduce the stepsize to a small C enough size in a reasonable number of attempts. May C indicate that f is very close to zero, or that there C is an error in some of the user-supplied values of f. C C IFAIL = 16: the stepsize has become so small that reasonable C progress cannot be made. Possible causes: same as for C IFAIL = 15. C C IFAIL = 17: one of NZMAX, NFMAX is too small to contain all C the points and function values required to do the C integral. C C .. Scalar Arguments .. DOUBLE COMPLEX ANSWER,ELAM,ELAM1,ELAM2,F DOUBLE PRECISION CONST INTEGER IFAIL,IREVCM,NFMAX,NZACT,NZMAX C .. C .. Array Arguments .. DOUBLE COMPLEX FARRAY(0:NFMAX,1:2),ZARRAY(0:NZMAX) C .. C .. Local Scalars .. DOUBLE COMPLEX FOLD,HMAX,HSTEP,LGQUOT,QUOT,Z,ZOLD DOUBLE PRECISION ARGQ,ARGQH,COLD,EPS INTEGER I,ISTEP,NQUAD,NMIN LOGICAL GTFAR1,GTFAR2,GTFAR3,GTFAR4 C .. C .. External Functions .. DOUBLE PRECISION X02MEP EXTERNAL X02MEP C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,CMPLX,DBLE,DIMAG,LOG,MAX,MIN C .. C .. Save statement .. SAVE I,NQUAD,HMAX,HSTEP,ZOLD,Z,FOLD,COLD,QUOT,LGQUOT,ISTEP,EPS C .. C WRITE (6,FMT=*) 'Successfully called BXSIDE' IF (IREVCM.EQ.0) RETURN GO TO (10,20,50) IREVCM 10 ZOLD = ELAM1 ZARRAY(0) = ELAM1 I = 0 ISTEP = 0 NQUAD = 100 NMIN = 10 EPS = 1000.D0*X02MEP(0.D0) HMAX = (ELAM2-ELAM1)/NMIN C WRITE (6,FMT=*) 'ELAM2,ELAM1,HSTEP:',ELAM2,ELAM1,HSTEP ELAM = ZOLD IREVCM = 2 RETURN 20 FARRAY(0,1) = F FARRAY(0,2) = CONST FOLD = F COLD = CONST ANSWER = CMPLX(0.D0,0.D0) C Get initial stepsize HSTEP = (ELAM2-ELAM1)/NQUAD 30 IF (ABS(HSTEP).GT.1.D-2) THEN HSTEP = HSTEP/10.D0 C WRITE (6,FMT=*) 'HSTEP (provisional):',HSTEP GO TO 30 END IF 40 Z = ZOLD + HSTEP C WRITE (6,FMT=*) 'HSTEP:',HSTEP C WRITE (6,FMT=*) 'ZOLD + HSTEP:',ZOLD + HSTEP C Trap possible overshooting of the interval of integration GTFAR1 = DBLE(ELAM1) .LT. DBLE(ELAM2) .AND. & DBLE(Z) .GT. (DBLE(ELAM2)+EPS) GTFAR2 = DBLE(ELAM1) .GT. DBLE(ELAM2) .AND. & DBLE(Z) .LT. (DBLE(ELAM2)-EPS) GTFAR3 = DIMAG(ELAM1) .LT. DIMAG(ELAM2) .AND. & DIMAG(Z) .GT. (DIMAG(ELAM2)+EPS) GTFAR4 = DIMAG(ELAM1) .GT. DIMAG(ELAM2) .AND. & DIMAG(Z) .LT. (DIMAG(ELAM2)-EPS) C WRITE (7,FMT=*) 'GTFAR1,GTFAR2,GTFAR3,GTFAR4:',GTFAR1,GTFAR2, C & GTFAR3,GTFAR4 IF (GTFAR1 .OR. GTFAR2 .OR. GTFAR3 .OR. GTFAR4) THEN C We have stepped too far: reset Z to the end of the integration C range. Z = ELAM2 END IF ELAM = Z IREVCM = 3 RETURN 50 QUOT = F/FOLD LGQUOT = LOG(QUOT) ARGQ = ATAN2(DIMAG(QUOT),DBLE(QUOT)) ARGQ = ABS(ARGQ) WRITE (7,FMT=*) 'ZOLD:',ZOLD WRITE (7,FMT=*) 'FOLD,COLD:',FOLD,COLD WRITE (7,FMT=*) 'Z:',Z WRITE (7,FMT=*) 'F,CONST:',F,CONST WRITE (7,FMT=*) 'QUOT,ARGQ:',QUOT,ARGQ WRITE (7,FMT=*) ARGQH = ARGQ/ABS(HSTEP) IF (ARGQ.GT.1.0D0) THEN C Step was too large WRITE (7,FMT=*) 'Step too large' HSTEP = HSTEP/ (1.4D0*ARGQ) ISTEP = ISTEP + 1 IF (ISTEP.LT.20) GO TO 40 IFAIL = 15 IREVCM = 0 RETURN END IF IF (ABS(HSTEP).LT.1.D-14) THEN C The stepsize is too small to make decent progress C IFAIL = 16 IREVCM = 0 RETURN END IF C If we are here then the step was small enough C Can we increase the stepsize? IF (ISTEP.EQ.0 .AND. ARGQ.LT.0.6D0) THEN C We may consider increasing the stepsize for the next step C WRITE (6,FMT=*) 'HSTEP,HMAX:',HSTEP,HMAX HSTEP = HSTEP/MAX(0.2D0,1.2D0*ABS(ARGQ)) IF (ABS(HSTEP).GT.ABS(HMAX)) & HSTEP = HSTEP*ABS(HMAX)/ABS(HSTEP) END IF ANSWER = ANSWER + LGQUOT I = I + 1 ZARRAY(I) = Z FARRAY(I,1) = F FARRAY(I,2) = CMPLX(COLD,0.D0) WRITE (7,FMT=*) 'ELAM,F,CONST:',Z,F,REAL(COLD) NZACT = I ZOLD = Z FOLD = F COLD = CONST IF (ISTEP.GT.0) ISTEP = ISTEP - 1 IF (ABS(Z-ELAM2).GT.EPS) THEN C We have not yet reached the end of the range of integration IF (I+1.LT.MIN(NZMAX,NFMAX)) GO TO 40 C Out of storage! IFAIL = 17 IREVCM = 0 RETURN END IF C If we are here then we have reached the end of the range C of integration. IREVCM = 0 IFAIL = 0 RETURN END C ------------------------------------------------------------------ SUBROUTINE NZERO(XMESH,NMESH,IMATCH,RELAM,SL4COF,SL4BCS, & SLSMPC,SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQHM, & IFAIL) C C This routine computes the number of eigenvalues of our C SLP whose real part is less than RELAM C C .. Parameters .. INTEGER ILOC PARAMETER (ILOC=10) C .. C .. Scalar Arguments .. DOUBLE PRECISION ARGQHM(2),RELAM INTEGER IFAIL,IMATCH,NMESH,NZ DOUBLE COMPLEX ZMIN,FMIN,ZMAX(2) C .. C .. Array Arguments .. DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF,SLSMPC,SLSMBC C .. C .. Local Scalars .. DOUBLE COMPLEX CSUM,CSUM2,ELAM,FUN,FUNSMP,ZEVAL, & ZMAX1(2),ZMAX2(2),FUN1 DOUBLE PRECISION ARGQH1(2),ARGQH2(2),CNSMP,CONST, & RHIGH,RLOW,AFMIN INTEGER I,IMLOC,IREVCM,IZC,NEVAL2,NEVALS C .. C .. Local Arrays .. DOUBLE PRECISION XLOC(0:10) C .. C .. External Subroutines .. EXTERNAL SLMISS,ZCOUNT C .. C .. External Functions .. DOUBLE PRECISION X01MPI EXTERNAL X01MPI C .. C .. Intrinsic Functions .. INTRINSIC DBLE,DIMAG,EXP,NINT C .. C Remark: On exit, ARGQHM(1) = most negative value of ARG(Q)/H, C ARGQHM(2) = most positive value of ARG(Q)/H, and ZMAX(1) = value C of lambda where ARGQHM(1) occurs, ZMAX(2) = value of lambda C where ARGQHM(2) occurs. C ZMIN = CMPLX(RELAM,0.D0) AFMIN = 1.D40 FMIN = CMPLX(AFMIN,AFMIN) XLOC(0) = XMESH(0) XLOC(ILOC) = XMESH(NMESH) DO 10 I = 1,ILOC XLOC(I) = (DBLE(ILOC-I)*XLOC(0)+DBLE(I)*XLOC(ILOC))/DBLE(ILOC) 10 CONTINUE IMLOC = NINT(DBLE(ILOC)/2.D0) RLOW = 0.D0 C Set least acceptable approximation to +infinity: RHIGH = 1.D8 IREVCM = 1 IFAIL = 0 20 CALL ZCOUNT(RLOW,RHIGH,RELAM,ZEVAL,FUN,IREVCM,CSUM,IZC,NEVALS, & ZMAX1,ARGQH1,IFAIL) IF (IFAIL.NE.0) THEN RETURN END IF IF (IREVCM.NE.0) THEN ELAM = ZEVAL CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,FUN,CONST, & IFAIL) IF (IFAIL.NE.0) RETURN C The reference problem has constant coefficients so C we can get away with using a really coarse mesh to C find its miss-distance function. CALL SLMISS(XLOC,ILOC,IMLOC,ELAM,SLSMPC,SLSMBC,FUNSMP,CNSMP, & IFAIL) IF (IFAIL.NE.0) RETURN FUN = (FUN/FUNSMP) C WRITE (7,FMT=*) ABS(FMIN),ABS(FUN),ELAM,ZMIN IF (ABS(FUN).LT.AFMIN) THEN ZMIN = ELAM FMIN = FUN AFMIN = ABS(FMIN) END IF C WRITE (9,FMT=*) 'Z,F(Z):',REAL(ELAM),REAL(DIMAG(ELAM)), C & REAL(FUN),REAL(DIMAG(FUN)) C WRITE (7,FMT=*) GO TO 20 END IF RLOW = 0.D0 C Set greatest acceptable approximation to -infinity: RHIGH = -1.D8 IREVCM = 1 IFAIL = 0 30 CALL ZCOUNT(RLOW,RHIGH,RELAM,ZEVAL,FUN1,IREVCM,CSUM2,IZC,NEVAL2, & ZMAX2,ARGQH2,IFAIL) IF (IFAIL.NE.0) THEN RETURN END IF IF (IREVCM.NE.0) THEN ELAM = ZEVAL C WRITE (9,FMT=*) 'Calling SLMISS (1)' CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,FUN1,CONST, & IFAIL) IF (IFAIL.NE.0) RETURN C The reference problem has constant coefficients so C we can get away with using a really coarse mesh to C find its miss-distance function. C WRITE (9,FMT=*) 'Call SLMISS (2)' CALL SLMISS(XLOC,ILOC,IMLOC,ELAM,SLSMPC,SLSMBC,FUNSMP,CNSMP, & IFAIL) IF (IFAIL.NE.0) RETURN FUN1 = FUN1/FUNSMP C WRITE (9,FMT=*) 'Z,F(Z):',REAL(ELAM),REAL(DIMAG(ELAM)), C & REAL(FUN1),REAL(DIMAG(FUN1)) C WRITE (7,FMT=*) ABS(FMIN),ABS(FUN1),ELAM,ZMIN IF (ABS(FUN1).LT.AFMIN) THEN C WRITE (7,FMT=*) 'Updating ZMIN' ZMIN = ELAM FMIN = FUN1 AFMIN = ABS(FMIN) END IF GO TO 30 END IF C NZ = NINT(DIMAG(CSUM-CSUM2)) NZ = NINT(DIMAG(CSUM-CSUM2+0.5D0*LOG(FUN1/FUN)/X01MPI(0.D0))) WRITE (9,FMT=*) 'CSUM-CSUM2 + LOG(FUN1/FUN)/2Pi:', & CSUM-CSUM2 + 0.5D0*LOG(FUN1/FUN)/X01MPI(0.D0) WRITE (9,FMT=*) 'No of fn evals:',NEVALS + NEVAL2 C WRITE (9,FMT=*) '0.5D0*LOG(FUN1/FUN)/X01MPI(0.D0):', C & 0.5D0*LOG(FUN1/FUN)/X01MPI(0.D0) IF (ARGQH1(1).LT.ARGQH2(1)) THEN ZMAX(1) = ZMAX1(1) ARGQHM(1) = ARGQH1(1) ELSE ZMAX(1) = ZMAX2(1) ARGQHM(1) = ARGQH2(1) END IF IF (ARGQH1(2).GT.ARGQH2(2)) THEN ZMAX(2) = ZMAX1(2) ARGQHM(2) = ARGQH1(2) ELSE ZMAX(2) = ZMAX2(2) ARGQHM(2) = ARGQH2(2) END IF WRITE (9,FMT=*) WRITE (9,FMT=*) 'Re(Lambda):',RELAM WRITE (9,FMT=*) 'NZ:',NZ IFAIL = 0 RETURN END C ------------------------------------------------------------------ SUBROUTINE ZCOUNT(RLOW,RHIGH,RELAM,ZEVAL,FUN,IREVCM,CSUM,IZC, & NEVALS,ZMAX,ARGQHM,IFAIL) C C Program to do zero count based on integration along C the imaginary axis C C .. Scalar Arguments .. DOUBLE COMPLEX CSUM,FUN,ZEVAL DOUBLE PRECISION RELAM,RHIGH,RLOW INTEGER IFAIL,IREVCM,IZC,NEVALS C .. C .. Local Scalars .. DOUBLE COMPLEX CHIGH,CLOW,F,FOLD,LGQUOT,QUOT,Z,ZOLD DOUBLE PRECISION ARGQ,ARGQH,HSTEP,IMBIT,IMBITN,REBIT INTEGER ISTEP,NQUAD LOGICAL GFENUF C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,ATAN2,CMPLX,DBLE,DIMAG,LOG,MAX,NINT,REAL C .. C .. Save statement .. SAVE NQUAD,REBIT,TRNDIF,CLOW,CHIGH,CTERM,FOLD,F,ZOLD,Z,LGQUOT,ZTR, & ZTRN,ZN,QUOT,IMBIT,IMBITN,HSTEP,ARGQH,ISTEP C .. C .. Array Arguments .. DOUBLE COMPLEX ZMAX(2) DOUBLE PRECISION ARGQHM(2) C .. IF (IREVCM.EQ.0) RETURN C WRITE (7,FMT=*) 'Z,ZOLD',Z,ZOLD C WRITE (7,FMT=*) ' *** REBIT = ***',REBIT GO TO (10,20,40) IREVCM 10 REBIT = RELAM CLOW = CMPLX(REBIT,RLOW) CHIGH = CMPLX(REBIT,RHIGH) ARGQHM(1) = 0.D0 ARGQHM(2) = 0.D0 ZMAX(1) = CLOW ZMAX(2) = CLOW ZOLD = CLOW IMBIT = RLOW C IMBIT = RLOW/ (1.D0+ABS(RLOW)) C IMBIT = 1.D2*RLOW/ (1.D0+1.D2*ABS(RLOW)) IREVCM = 2 ZEVAL = ZOLD RETURN 20 FOLD = FUN C FOLD = FUN(ZOLD) CSUM = CMPLX(0.D0,0.D0) NQUAD = 100 NEVALS = 0 HSTEP = 1.D0/NQUAD IF (RLOW.GT.RHIGH) HSTEP = -HSTEP ISTEP = 0 30 IMBITN = IMBIT + HSTEP C* C Avoid overshooting ends: C* IF ((HSTEP.GT.0.D0.AND.IMBITN.GT.RHIGH) .OR. C* & (HSTEP.LT.0.D0.AND.IMBITN.LT.RHIGH)) THEN C* Z = CHIGH C* ELSE Z = CMPLX(REBIT,IMBITN) C* C Z = CMPLX(REBIT,IMBITN/ (1.D0-ABS(IMBITN))) C* C Z = CMPLX(REBIT,1.D-2*IMBITN/ (1.D0-ABS(IMBITN))) C* C* END IF C* IF ((DIMAG(Z).GT.RHIGH.AND.RLOW.LT.RHIGH) .OR. C* & (DIMAG(Z).LT.RHIGH.AND.RLOW.GT.RHIGH)) THEN C* Z = CMPLX(REBIT,RHIGH) C* END IF ZEVAL = Z IREVCM = 3 RETURN 40 F = FUN C F = FUN(Z) NEVALS = NEVALS + 1 QUOT = F/FOLD LGQUOT = LOG(QUOT) C WRITE (7,FMT=*) 'QUOT=',QUOT C WRITE (7,FMT=*) 'REAL(QUOT):',REAL(QUOT) C WRITE (7,FMT=*) 'IM(QUOT):',DIMAG(QUOT) C WRITE (7,FMT=*) 'ATAN2(0.00001,1.0000):', C & ATAN2(1.D-5,1.D0) ARGQ = ATAN2(DIMAG(QUOT),DBLE(REAL(QUOT))) ARGQH = ARGQ/HSTEP IF (ARGQH.LT.ARGQHM(1)) THEN ZMAX(1) = CMPLX(REBIT,DIMAG(Z+ZOLD)*0.5D0) ARGQHM(1) = ARGQH END IF IF (ARGQH.GT.ARGQHM(2)) THEN ZMAX(2) = CMPLX(REBIT,DIMAG(Z+ZOLD)*0.5D0) ARGQHM(2) = ARGQH END IF C WRITE (7,FMT=*) 'ZMAX(1),ZMAX(2):',ZMAX(1),ZMAX(2) C WRITE (10,FMT=*) 'Z,ARG(Q):',Z,ARGQ C WRITE (7,FMT=*) 'Z,ARG(Q):',Z, C & ARGQ ARGQ = ABS(ARGQ) C WRITE (7,FMT=*) REAL(DIMAG(Z)),REAL(F),REAL(DIMAG(F)) C WRITE (7,FMT=*) 'Z,ZOLD,F,FOLD:',Z,ZOLD,F,FOLD C WRITE (7,FMT=*) 'HSTEP=',HSTEP IF (ARGQ.GT.1.0D0) THEN C WRITE (7,FMT=*) 'Stepsize was too large' C WRITE (7,FMT=*) 'New stepsize is',HSTEP/(1.4D0*ARGQ) C Step was too large HSTEP = HSTEP/ (1.4D0*ARGQ) ISTEP = ISTEP + 1 IF (ISTEP.LT.10) GO TO 30 IFAIL = 5 IREVCM = 0 RETURN END IF IF (ABS(HSTEP).LT.1.D-14) THEN C The stepsize is too small to make decent progress C IFAIL = 6 IREVCM = 0 RETURN END IF C If we are here then the step was small enough C Can we increase the stepsize? C WRITE (7,FMT=*) 'Stepsize is',HSTEP C WRITE (7,FMT=*) 'ISTEP=',ISTEP IF (ISTEP.EQ.0 .AND. ABS(ARGQ).LT.0.6D0) THEN C We may consider increasing the stepsize for the next step C WRITE (7,FMT=*) 'Increasing stepsize' HSTEP = HSTEP/MAX(0.2D0,1.2D0*ABS(ARGQ)) C WRITE (7,FMT=*) 'New stepsize is',HSTEP C IF (ABS(HSTEP).GT.0.1D0) HSTEP = SIGN(1.D0,HSTEP)*0.1D0 END IF CSUM = CSUM + LGQUOT ZOLD = Z FOLD = F IMBIT = IMBITN IF (ISTEP.GT.0) ISTEP = ISTEP - 1 C Check to see if we have gone far enough: GFENUF = (RLOW.LT.RHIGH.AND.DIMAG(Z).GT.RHIGH).OR. & (RLOW.GT.RHIGH.AND.DIMAG(Z).LT.RHIGH) GFENUF = GFENUF.AND.(ABS(ARGQ).LT.0.1D0) IF (.NOT.GFENUF) GO TO 30 C** IF (ABS(RHIGH-DIMAG(Z)).GT.1.D-6) GO TO 30 C We have now finished the quadrature C Divide by 2*Pi CSUM = CSUM/ (8.D0*ATAN(1.D0)) IREVCM = 0 IZC = NINT(DIMAG(CSUM)) RETURN END C **** REPLACEMENTS FOR NAG ROUTINES FOR THOSE WHO DO NOT HAVE THE NAG **** C **** LIBRARY. DOUBLE PRECISION FUNCTION X02MEP(X) IMPLICIT NONE DOUBLE PRECISION X LOGICAL FRSTME DATA FRSTME/.TRUE./ DOUBLE PRECISION PREC SAVE PREC IF (FRSTME) THEN PREC = 1.D0/2.D0**53 X02MEP = PREC FRSTME = .FALSE. ELSE X02MEP = PREC END IF C WRITE (7,FMT=*) 'Value of machine precision was',X02MEP RETURN END DOUBLE PRECISION FUNCTION X01MPI(X) IMPLICIT NONE DOUBLE PRECISION X LOGICAL FRSTME DATA FRSTME/.TRUE./ DOUBLE PRECISION PI SAVE PI IF (FRSTME) THEN PI = 4.D0*ATAN(1.D0) X01MPI = PI FRSTME = .FALSE. ELSE X01MPI = PI END IF WRITE (7,FMT=*) 'Value of PI was',X01MPI RETURN END SUBROUTINE F06TFM(G,N,M,X,IX,XOUT,IXOUT) C WARNING: THIS ROUTINE DOES NOT HAVE THE FULL FUNCTIONALITY OF THE C NAG VERSION F06TFF IMPLICIT NONE CHARACTER*1 G INTEGER N,M,IXOUT,IX DOUBLE COMPLEX XOUT(IXOUT,*),X(IX,*) INTEGER I,J C NOTES: C C G: CHARACTER*1. NOT USED. INCLUDED FOR ONE-WAY COMPATIBILITY WITH C NAG. C C N,M: INTEGERS. The routine does the copy C C XOUT(I,J) = X(I,J) FOR I = 1,..,N, J = 1,..,M C C DO 10 J = 1,M DO 20 I = 1,N XOUT(I,J) = X(I,J) 20 CONTINUE 10 CONTINUE RETURN END SUBROUTINE F06THM(G,N,M,OFDIAG,DIAG,X,IX) C WARNING: THIS ROUTINE DOES NOT HAVE THE FULL FUNCTIONALITY OF THE C NAG VERSION F06THF IMPLICIT NONE CHARACTER*1 G INTEGER N,M,IX DOUBLE COMPLEX OFDIAG,DIAG,X(IX,*) INTEGER I,J C NOTES: C C G: CHARACTER*1. NOT USED. INCLUDED FOR ONE-WAY COMPATIBILITY WITH C NAG. C C N,M: INTEGERS. The routine does the assignment C C XOUT(I,J) = OFDIAG FOR I = 1,N, J = 1,M, I not equal to J; C XOUT(J,J) = DIAG FOR J = 1,M. C C DO 10 J = 1,M DO 20 I = 1,N X(I,J) = OFDIAG 20 CONTINUE X(J,J) = DIAG 10 CONTINUE RETURN END