*********************************************************************** * * * Routines for SLEDGE * * (Sturm-Liouville Estimates Determined by Global Errors) * * * *********************************************************************** C C Release 2.2 04/12/93 C C Steven Pruess, Colorado School of Mines C spruess@mines.colorado.edu C Charles Fulton, Florida Institute of Technology C fulton@zach.fit.edu C *********************************************************************** * These routines estimate eigenvalues, eigenfunctions and/or * * spectral density functions for Sturm-Liouville problems. * * The differential equation has the form: * * * * -(p(x)u')' + q(x)u = EV*r(x)u for x in [A,B] * * * * with boundary conditions (at regular points) * * * * A1*u - A2*(pu') = EV*(A1'*u - A2'*(pu')) at A * * B1*u + B2*(pu') = 0 at B . * * * * The functions p(x) and r(x) are assumed to be positive in * * the open interval (A,B). * *********************************************************************** C C Possible outputs are: C a set of eigenvalues; C a set of eigenvalues and tables of values for their eigen- C functions; C a table of the spectral density function (for cases with C continuous spectrum). C a classification of the problem (regular or singular; if C singular then limit point or limit circle, oscillatory C or nonoscillatory). C C The code can find eigenvalues and eigenfunctions for problems C in spectral category 1 (both endpoints NONOSC), spectral C category 2 (one endpoint NONOSC and the other O-NO), and C those discrete eigenvalues below the essential spectrum in C spectral category 10 (both endpoints O-NO). Here OSC at an C endpoint means the Sturm-Liouville equation is oscillatory for C all real values of EV at that endpoint, NONOSC at an endpoint C means the equation is nonoscillatory for all real values of EV C at that endpoint, and O-NO means there is a `cutoff' value EV' C such that the equation is nonoscillatory for real values of C EV < EV' and oscillatory for real values of EV > EV'. For C problems in other spectral categories an error return will C be generated. The manner in which SLEDGE classifies singular C endpoints of Sturm-Liouville problems as LP/LC (Limit Point/ C Limit Circle), OSC/NONOSC/O-NO, and uses this information to C determine the spectral category is explained in detail in C reference [2]. C C There is one subroutine called SLEDGE of direct interest to the C user; additionally, a secondary routine INTERV is available which C determines the indices of eigenvalues located in a specified C subinterval of the real line. C C The names of other routines in this package are AITKEN, ASYMEV, C ASYMR, BRCKET, CLASS, CLSEND, DENSEF, DSCRIP, EXTRAP, GETEF, C GETRN, MESH, POWER, PQRINT, REGULR, SHOOT, START, STEP, and C ZZERO. C C There are 4 blocks of labeled COMMON with the names SLREAL, C SLINT, SLLOG, and SLCLSS. C C This is the double precision version of the code; all floating C point variables should be declared DOUBLE PRECISION in the C calling program. In these subprograms all such local C variables and constants have been explicitly declared; also, C FORTRAN77 generic intrinsic functions have been used, so C conversion to single precision should be straightforward, if C desired. C C ACKNOWLEDGMENT: This work was partially supported by the C National Science Foundation under grants DMS-8813113 and DMS- C 8905202 to Florida Institute of Technology and DMS-8800839 and C DMS-8905232 to the Colorado School of Mines. C C References C C The following papers are available from the authors on request: C C [1]. Pruess & Fulton, Mathematical software for Sturm-Liouville C problems, ACM Trans. on Math. Software, to appear, 1993. C C [2]. Fulton, Pruess & Xie, The automatic classification of Sturm- C Liouville problems, submitted, 1992. C C [3]. Pruess, Fulton & Xie, An asymptotic numerical method for a C class of singular Sturm-Liouville problems, submitted, 1992. C C [4] Fulton and Pruess, Eigenvalue and eigenfunction asymptotics C for regular Sturm-Liouville problems, Jour. Math. Anal. and C Appls., to appear. C C [5] Fulton and Pruess, Numerical Approximation of singular C spectral functions arising from the Fourier-Jacobi problem C on a half line with continuous spectra, Sixth International C Workshop in Analysis and its Applications, June, 1992. C [6]. Pruess, Fulton & Xie, Performance of the Sturm-Liouville C software package SLEDGE, Colo. School of Mines, Dept. of C Math. and Comp. Sci., MCS-91-19, 1991. Revision 12/92. C----------------------------------------------------------------------- C Brief overview of algorithms: C C The code constructs (or takes from input) an initial mesh, C called the level 0 mesh. Subsequent meshes (for level 1,2,...) C are unions of the previous level's mesh with its midpoints. A C sequence of estimates for desired eigenvalues and eigenfunctions C is constructed, one set for each level. These estimates (the C eigenvalue is called EvHat) are exact solutions (up to the C requested tolerance) of a Sturm-Liouville problem which is an C approximation to the original one; this approximation results C from replacing the given coefficient functions with step function C approximations relative to the current level's mesh. The eigen- C functions of the resulting ODE's are piecewise trigonometric C (circular or hyperbolic) functions. C If estimates for the spectral density function are reqested, C these are computed as limits of a sequence of spectral density C functions of approximating regular problems. For these regular C problems the spectral density function is a step function, and C is computed directly from the definition making use of computed C eigenvalues and the norm reciprocals of the corresponding eigen- C functions. If verbose output is rquested by the user, there C will be displayed iterations (corresponding to the sequence of C approximating regular intervals which the code automatically C selects) and within each iteration there will be levels C (corresponding to increasingly finer meshes as described above). C A step spectral density function will be printed at each level C of each iteration. The spectral density function displayed at C the end of each iteration is the result of an h-squared extra- C polation over the regular step functions generated at each level C of this iteration. The condition for stopping at a given iter- C ation is a straightforward comparison of the spectral function C data for the current iteration with the previous iteration. There C is no extrapolation over the sequence of regular approximating C intervals as no extrapolation theory for the approximation of C the singular spectral function by regular step spectral functions C is known. (To achieve closer approximation of the regular step C spectral functions to the singular spectral function, it is C actually the piecewise linear function obtained by joining the C midpoints of successive steps by a straight line which is used as C the `regular' spectral function for the purpose of generating the C actual data used for the h-squared extrapolation.) C The classification is determined by applying standard theory C to an approximating problem, each of whose coefficient functions, C in a small neighborhood of each endpoint, consists of the leading C term in a power-like asymptotic development. For this reason C there are many problems, particularly those with oscillatory C coefficient functions, for which the code's output for the C classification information is labelled `uncertain'. For further C information on the theory used by the code to generate endpoint C classifications and spectral category information see [2] above. C----------------------------------------------------------------------- C Usage (simple explanation) - C The subroutine SLEDGE is called in the following manner: C C SUBROUTINE SLEDGE(JOB,CONS,ENDFIN,INVEC,TOL,TYPE,EV,NUMX,XEF,EF, C PDEF,T,RHO,IFLAG,STORE) C C If k eigenvalues (no eigenvectors) are sought then set C (a) the logical 5-vector JOB to (True, False, False, False, True); C (b) the real 8-vector CONS to the values of A1,A1',A2,A2',B1,B2, C A,B for the boundary condition information. It does not C matter what values are used for infinite endpoints, nor for C the boundary constants at a singular endpoint; the code C automatically selects the Friedrichs' boundary condition at C NONOSC singular endpoints, overriding user input for the C boundary condition constants; for infinite endpoints the code C also automatically selects these constants. C (c) the logical 2-vector ENDFIN to (True, True) if both endpoints C are finite, (True, False) if A is finite but B infinite, etc.; C (d) the integer vector INVEC should have C INVEC(1) = 0 (no internal printing) C INVEC(2) = 0 C INVEC(3) = k, the number of eigenvalues sought C INVEC(3+i) = index of ith eigenvalue sought,i = 1,...,k; C (e) the real 6-vector TOL should have C TOL(1) = absolute error tolerance desired, C TOL(2) = relative error tolerance desired, C the remaining 4 entries of TOL are ignored; C (f) the output estimate for the ith eigenvalue is returned C in EV(i), i = 1,...,k; C (g) the output integer k-vector IFLAG(*) should have all entries C zero; nonzero values indicate warnings or error returns and C are explained in the detailed usage section below; C (h) the auxiliary vector STORE(*) should be dimensioned at least C 155 in the calling program; C (i) the logical 4 by 2 vector TYPE, the real vectors XEF(1), C EF(1), PDEF(1), T(1), and RHO(1) can be ignored except that C they need to be declared in the calling program. The integer C scalar NUMX can also be ignored. C C If k eigenfunctions are also desired, then follow the above C pattern except make JOB(1) False and JOB(2) True. The values of C TOL(3) and TOL(4) control the absolute and relative errors in C each u(x); TOL(5) and TOL(6) control the absolute and relative C errors in each (pu')(x). It is usually appropriate to set TOL(5) C = TOL(3) = TOL(1) and TOL(6) = TOL(4) = TOL(2), but the user has C the option of entering all six tolerance parameters as desired. C The output eigenfunction information is returned in the three C real vectors X(*) for the independent variable x , EF(*) for u(x), C PDEF(*) for (pu')(x). The code automatically chooses the x C values; the number of values is returned in NUMX. If you prefer C another choice of output points, see the detailed explanation C below on usage of the code. The values for the first requested C u(x) are returned in the first NUMX locations of EF(*), those for C the second are in the next NUMX locations, etc. PDEF(*) is part- C itioned similarly; X(*) must be dimensioned at least 31 in the C calling program while EF(*) and PDEF(*) must be dimensioned at C least 31*k. The auxiliary vector STORE(*) should be dimensioned C at least 420. C C For other possibilities, see the detailed description which C follows. C----------------------------------------------------------------------- C Usage (detailed explanation) - C C SUBROUTINE SLEDGE(JOB,CONS,ENDFIN,INVEC,TOL,TYPE,EV,NUMX,XEF,EF,PDEF, C T,RHO,IFLAG,STORE) C C Input parameters; C JOB(*) = logical 5-vector, C JOB(1) = .True. iff a set of eigenvalues are to be C computed but not their eigenfunctions. C JOB(2) = .True. iff a set of eigenvalue and eigenfunc- C tion pairs are to be calculated. C JOB(3) = .True. iff the spectral function is to be C computed over some subinterval of the C essential spectrum. C JOB(4) = .True. iff the normal call to the routines for C classification (regular/singular, etc.) C is OVERRIDDEN. If JOB(4) is True then C TYPE(*,*) discussed below must be C INPUT correctly! Most users will not C want to override the classification C routines, but it would, of course, be C appropriate for users experimenting with C problems for which the coefficient C functions do not have power-like C behavior near the singular endpoints. C Note: the code may perform poorly if C the classification information is C incorrect; since the cost is usually C negligible, it is strongly recommended C that JOB(4) be False. The classifica- C tion is deemed sufficiently important C for spectral density function calcul- C ations that JOB(4) is ignored when C the input JOB(3) is True. C JOB(5) = .True. iff mesh distribution is to be chosen C by SLEDGE. If JOB(5) is True and NUMX C is zero, then the number of mesh C points is also chosen by SLEDGE; if C NUMX > 0 then NUMX mesh points will be C used. If JOB(5) is False, then the C number (NUMX) and distribution C (XEF(*)) must be input by the user. C If JOB(3) is True and JOB(5) False C then the user must set BOTH the number C NUMX and distribution. In this case, C NO global error estimates are made. C CONS(*) = real vector of length 8, values are the boundary C condition constants A1, A1', A2, A2', B1, B2, A, B. C In the case of a NONOSC singular endpoint, the class- C ification routine uses the Friedrichs' boundary C condition constants. The code cannot automatically C choose a non-Friedrichs' boundary condition; however, C interval truncation in the user's calling program can C be used, together with many calls to SLEDGE, to C compute singular eigenvalues associated with a non- C Friedrichs' boundary condition at a NONOSC endpoint C (see remark 12 below). C ENDFIN(*) = logical 2-vector, values are C ENDFIN(1) = .True. iff endpoint A is finite. C ENDFIN(2) = .True. iff endpoint B is finite. C INVEC(*) = integer vector of length 3+(number of eigenvalues C desired). This vector contains a variety of input C information. C INVEC(1) controls the amount of internal printing: values are C from 0 (no printing) to 5 (much printing). C For INVEC(1) > 0 much of the output will be to a file C attached to unit #21 which should be named in the C user's calling program via an OPEN statement. C Output for the various cases is, when INVEC(1) = C 0 no printing. C When JOB(1) or JOB(2) is True C 1 initial mesh (the first 51 or fewer points), C eigenvalue estimate at each level, C 4 the above, C at each level C matching point for eigenfunction shooting, C X(*), EF(*), PDE(*) values, C 5 all the above, C at each level C brackets for the eigenvalue search, C intermediate shooting info for the eigen- C function and eigenfunction norm. C When JOB(3) is True C 1 the actual (a,b) used at each iteration, C the total number of eigenvalues computed, C 2 the above, C switchover points to the asymptotic formulas, C some intermediate Rho(t) approximations, C 3 all the above, C initial meshes for each iteration, C index of the largest EV which may be computed, C various Ev and RsubN values, C 4 all of the above, C RhoHat values at each level, C 5 all of the above, C all Ev and RsubN values below switchover point. C When JOB(4) is False C 2 output a description of the spectrum, C 3 the above plus the constants for the C Friedrichs' boundary condition(s), C 5 all the above plus intermediate details of C classification calculation. C Some of the output may go to the default output device C (screen or printer), but all information requested is C also directed to the file attached to unit #21. C INVEC(2) gives the number (positive) of output values desired C for the array RHO(*) (not referenced if JOB(3) is C False). C INVEC(3) is the total number of eigenvalues to be output in C EV(*). C INVEC(J) for J = 4, 5, ..., 3+INVEC(3) contains the indices for C the eigenvalues sought. If JOB(1) and JOB(2) are C False, this part of INVEC(*) is not referenced. C TOL(*) = real vector of from 2 to 6 tolerances. C If JOB(1) or JOB(2) is True then C TOL(1) is the absolute error tolerance for e-values, C TOL(2) is the relative error tolerance for e-values, C TOL(3) is the abs. error tolerance for e-functions, C TOL(4) is the rel. error tolerance for e-functions, C TOL(5) is the abs. error tolerance for eigenfunction C derivatives, C TOL(6) is the rel. error tolerance for eigenfunction C derivatives. C Eigenfunction tolerances need not be set if JOB(2) C is False. C If JOB(3) is True then C TOL(1) is the absolute error tolerance, C TOL(2) is the relative error tolerance; C the output RHO values are NOT required to satisfy C these tolerances when JOB(5) is False. C All absolute error tolerances must be positive; all C relative error tolerances must be at least 100 times C the unit roundoff. C NUMX = integer whose value is C the number of output points where each eigen- C function is to be evaluated (the number of entries C in XEF(*)) when JOB(2) is True, C or C the number of points in the initial mesh used when C JOB(5) is False and NUMX>0. C If JOB(5) is False, the points in XEF(*) should be C chosen to have a reasonable distribution. Since the C endpoints A and B must be part of any mesh, NUMX C cannot be 1 in this case. If JOB(5) is FALSE and C JOB(3) is True, then NUMX must be positive. C XEF(*) = real vector of points where C eigenfunction estimates are desired (JOB(2) True) C or C where user's initial mesh is entered (JOB(5) False C and NUMX>0). C The values must satisfy C A = XEF(1) < XEF(2) < ... < XEF(NUMX) = B . C When JOB(2) is True the initial mesh corresponds to C the set of points where eigenfunction output is C desired. If JOB(2) is False and NUMX = 0, then this C vector is not referenced. When A and/or B are C infinite (as indicated through ENDFIN(*)), the C entries XEF(1) and/or XEF(NUMX) are ignored; however, C it is required that XEF(2) be negative when ENDFIN(1) C is False, and XEF(NUMX-1) be positive when ENDFIN(2) C is False (otherwise, IFLAG = -39 will result). C T(*) = real vector of INVEC(2) values where the spectral C function RHO(*) is desired (the existence and location C of continuous spectrum can be found by first calling C SLEDGE with JOB(J) False, J=1,...,4 and INVEC(1) = 1). C Vector T(*) is not referenced if JOB(3) is False. Its C entries must be in increasing order. C C Output parameters: C TYPE(*,*) = 4 by 2 logical array; column 1 carries information C about endpoint A while column 2 refers to B. C TYPE(1,*) = True iff the endpoint is regular, C TYPE(2,*) = True iff it is limit circle, C TYPE(3,*) = True iff it is nonoscillatory for all EV, C TYPE(4,*) = True iff it is oscillatory for all EV, C Important note: all of these must be correctly INPUT C if JOB(4) is True! C EV(*) = real vector containing the computed approximations to C the eigenvalues whose indices are specified in C INVEC(*); if JOB(1) and JOB(2) are False, then the C output has no meaning. C NUMX = the number of output points for eigenfunctions when C input NUMX = 0, and JOB(2) or JOB(5) is True. C XEF(*) = input values (if any) are changed only if JOB(2) and C JOB(5) are True; in this case, the output values C are chosen by the code. If JOB(2) is False then this C vector is not referenced; if JOB(2) is True and NUMX>0 C on input then XEF(*) should be dimensioned at least C NUMX+16 in the calling program. If JOB(2) is True and C NUMX=0 on input (so that the code chooses NUMX), then C dimension XEF(*) at least 31 in the calling program. C EF(*) = real vector of eigenfunction values: EF((k-1)*NUMX+i) C is the estimate of u(XEF(i)) corresponding to the C eigenvalue in EV(k). If JOB(2) is False then this C vector is not referenced. Otherwise, if JOB(2) is C True and NUMX>0 on input then EF(*) should be C dimensioned at least NUMX*INVEC(3) in the calling C program. If JOB(2) is True and NUMX=0 on input (so C that the code chooses NUMX), then dimension XEF(*) C at least 31*INVEC(3) in the calling program. C PDEF(*) = real vector of eigenfunction derivative values: C PDEF((k-1)*NUMX+i) is the estimate of (pu')(XEF(i)) C corresponding to the eigenvalue in EV(k). If JOB(2) C is False then this vector is not referenced; otherwise, C it must be dimensioned as is EF(*). C RHO(*) = real vector of values for the spectral density C function rho(t), RHO(I) = rho(T(I)). RHO(*) must be C dimensioned at least INVEC(2); this vector is not C referenced if JOB(3) is False. C IFLAG(*) = integer vector carrying information about the output. C Declared length must be at least max(1,INVEC(3)). For the Kth C requested eigenvalue (when JOB(1) or JOB(2) is true; otherwise, C only IFLAG(1) is used): C IFLAG(K) = 0, normal return, output should be reliable. C < 0, fatal error, calculations ceased: if C = -1, too many levels needed for the eigenvalue C calculation; problem seems too difficult for C this algorithm at this tolerance. Are the C coefficient functions nonsmooth? C = -2, too many levels needed for the eigenfunction C calculation; problem seems too difficult for C this algorithm at this tolerance. Are the C eigenfunctions ill-conditioned? C = -3, too many levels needed for the spectral density C calculation; problem seems too difficult for C this algorithm at this tolerance. C = -4, the user has requested the spectral density C function for a problem which has no continuous C spectrum. C = -5, the user has requested the spectral density C function for a problem with both endpoints C generating essential spectrum, i.e., both C endpoints being either OSC or O-NO. The spectral C density function calculation has not been C implemented for such cases. For spectral C category 10 (both endpoints O-NO) the spectral C multiplicity is generally two, proper normal- C izations for the solutions against which the C spectral functions will be normalized will C depend on how the user wants to express the C eigenfunction expansion. Users having problems C in spectral category 10 are encouraged to supply C them to the authors, and if possible, recommend C normalizations of the two solutions to be used in C writing the associated eigenfunction expansion. C = -6, the user has requested the spectral density C function for a problem in spectral category 2 for C which a proper normalization of solution at the C NONOSC endpoint is not known; for example, C problems with an irregular singular point or C infinite endpoint at one end and continuous C spectrum generated at the other. Users with C problems of this type are encouraged to supply C them to the authors, and if possible, recommend a C normalization of solution at the NONOSC endpoint C which they would like to see implemented. As a C rule it is best to pick a normalization which C ensures that the solution is uniquely fixed and C entire in the eigenvalue parameter EV for all x C in the Sturm-Liouville interval; for further C mathematical information on NONOSC endpoints we C refer to paper [2] above. C = -7, problems encountered in obtaining a bracket. C = -8, too small a step used in the integration; C TOL(*) values may be too small for this problem. C = -9, too small a step used in a spectral density C function calculation for which the continuous C spectrum is generated by a finite endpoint. Try C transforming to Liouville (or some other) form. C = -10, an argument to the circular trig functions is C too large. Try rerunning with a finer initial C mesh, or, on singular problems, use interval C truncation (see remark (12)). C = -15, p(x) and r(x) not positive in (A,B). C = -20, eigenvalues/functions were requested for a C problem with an OSC singular endpoint. C Interval truncation (see remark (12)) must be C used on such problems. C = -3?, illegal input, viz. C -30, NUMX = 1 when JOB(5) is True, C or NUMX = 0 when JOB(3) is True and JOB(5) is C False, C -31, B1 = B2 = 0 (at a regular endpoint), C -32, A1'*A2-A1*A2' .le. 0 when A1' or A2' nonzero, C -33, A1 = A2 = A1'= A2'= 0 (at a regular endpoint), C -34, A .ge. B (when both are finite), C -35, TOL(odd) .le. 0 , C -36, TOL(even) < 100*unit roundoff, C -37, INVEC(k) < 0 for some k>3 when INVEC(3)>0, C -38, INVEC(2) .le. 0 when JOB(3) is True , C -39, XEF(*) entries out of order or not in [A,B]. C or XEF(2), XEF(NUMX-1) have the wrong sign in C infinite interval cases, C or T(*) entries are out of order. C > 0, indicates some kind of warning, in this case the C value may contain ANY of the following digits: C = 1, failure in routine BRCKET probably due to a C cluster of eigenvalues which the code cannot C separate. Calculations have continued as best C as possible, but any eigenfunction results are C suspect. Try rerunning with tighter input C tolerances to separate the cluster. C = 2, there is uncertainty in the classification for C this problem. Because of the limitations of the C floating point arithmetic on the computer used, C and the nature of the finite sampling, the C routine is cannot be decisive about the C classification information at the requested C tolerance. C = 3, there may be some eigenvalues imbedded in the C essential spectrum; using IPRINT greater than C zero will result in additional output giving C the location of the approximating eigenvalues C for the step function problem. These could be C extrapolated to estimate the actual eigenvalue C embedded in the essential spectrum. C = 5, a change of variables was made to avoid poten- C tial slow convergence; however, the global C error estimates may not be as reliable. Some C experimentation using different tolerances is C recommended. C = 6, there were problems with eigenfunction conver- C gence in a spectral density calculation; the C output Rho(t) may not be accurate. C C Auxiliary storage: C STORE(*) = real vector of auxiliary storage, must be dimensioned C at least C max(155,NUMX+16) in general; C 26*(NUMX+16) for any eigenfunction calculation; C 2400+13*INVEC(2) for any spectral density calculation. C----------------------------------------------------------------------- C SUBROUTINE INTERV(FIRST,ALPHA,BETA,CONS,ENDFIN,NFIRST,NTOTAL, C IFLAG,STORE) C C Input parameters: C FIRST = logical; value is True if various internal variables C have not yet been set. If a prior call has been made C to INTERV with FIRST True, then a little time can C be saved by letting FIRST be False. C IMPORTANT NOTE: setting FIRST = True will clobber any C initial mesh the user has input (when NUMX > 0 or C JOB(5) is False); also, INTERV will classify the C problem irregardless of what JOB(4) is set to C for SLEDGE. C ALPHA = real value of left end point of search interval. C BETA = real value of right end point of search interval. C CONS(* ) = real vector of 8 input constants: A1, A1', A2, A2', C B1, B2, A, B. C ENDFIN(*) = logical 2-vector, same meaning as in SLEDGE. C STORE(*) = real vector holding initial mesh. C C Output parameters: C NFIRST = index of first eigenvalue > ALPHA. C NTOTAL = total number of eigenvalues in the interval. C IFLAG = integer status indicator. C IFLAG = 0 , normal return, output should be reliable, C = 11 , there are no eigenvalues in [alpha, beta], C = 12 , low confidence in NFIRST or NTOTAL or both, C = 13 , BETA and/or ALPHA exceed the cutoff for C the continuous spectrum. If only BETA C is too big then NFIRST may be OK, but C NTOTAL is meaningless. C = -11 , ALPHA .ge. BETA, C = -25 , oscillatory endpoint, output meaningless, C = -3? , illegal CONS(*) values (see above comments C on SLEDGE for an explanation). C-------------------------------------------------------------------------- C In addition, a subroutine subprogram must be provided for the C coefficient functions p(x), q(x), and r(x); the form of this C routine is C C SUBROUTINE COEFF(X,PX,QX,RX) C DOUBLE PRECISION X,PX,QX,RX C ... C PX = ... C QX = ... C RX = ... C RETURN C END C C The subroutine name MUST be COEFF, though of course the names of C arguments only need follow the usual FORTRAN77 rules. X is the C independent variable; PX, QX, and RX are the output values of the C respective coefficient functions p(x), q(x), and r(x) at X. C----------------------------------------------------------------------- C This is a simple sample driver for SLEDGE. CC CC Declare all variables: CC C INTEGER IFLAG(1),INVEC(4),NUMX, I,J,K C LOGICAL JOB(5),TYPE(4,2),ENDFIN(2) C DOUBLE PRECISION CONS(8),TOL(6),EV(1),T(3),RHO(3),STORE(2450), C & XEF(5),EF(5),PDEF(5) CC CC Load the boundary condition information into CONS(*). CC This example has a Neumann condition at A = 1, and a CC singular point at B = +infinity. CC C DATA CONS/0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0/ C DATA ENDFIN/.TRUE., .FALSE./ CC CC The eigenfunctions will be estimated at 5 points. CC C DATA NUMX,XEF/5, 1.0, 1.5D0, 2.0, 4.0, 100.0/ CC CC Initialize the vector INVEC(*): CC little printing, CC 3 output points for the density function Rho(t), CC estimates for the first (index 0) eigenvalue/function. CC C DATA INVEC/1, 3, 1, 0/ CC CC Set the JOB(*) vector: CC estimate both eigenvalues and eigenvectors, CC estimate the spectral density function, CC classify, CC force the initial mesh to be the output points. CC C DATA JOB/.FALSE.,.TRUE.,.TRUE.,.FALSE.,.FALSE./ CC CC Set the tolerances: CC C DATA TOL/1.D-5,1.D-4, 1.D-5,1.D-4, 1.D-5,1.D-4/ CC CC Initialize the 3 output points for the density function. CC C DATA T/0.0, 0.5, 2.0/ CC CC Open file for output. CC C OPEN(21,FILE = 'sample.out') C CALL SLEDGE(JOB,CONS,ENDFIN,INVEC,TOL,TYPE,EV,NUMX,XEF,EF,PDEF, C & T,RHO,IFLAG,STORE) CC CC Print results: CC C DO 30 I = 1,INVEC(3) C WRITE (*,10) INVEC(3+I),EV(I),IFLAG(I) C WRITE (21,10) INVEC(3+I),EV(I),IFLAG(I) C 10 FORMAT(' Nev =',I6,'; Ev =',D25.15,'; Flag = ',I3) C IF (IFLAG(I) .GT. -10) THEN C WRITE (*,15) C WRITE (21,15) C 15 FORMAT(13X,'x',23X,'u(x)',18X,'(pu`)(x)') C K = NUMX*(I-1) C DO 25 J = 1,NUMX C WRITE (21,20) XEF(J),EF(J+K),PDEF(J+K) C 20 FORMAT(3D25.15) C 25 CONTINUE C ENDIF C 30 CONTINUE C WRITE (*,35) C WRITE (21,35) C 35 FORMAT(/,8X,'t',21X,'Rho(t)') C DO 45 I = 1,INVEC(2) C WRITE (*,40) T(I),RHO(I) C WRITE (21,40) T(I),RHO(I) C 40 FORMAT(F11.3,D32.15) C 45 CONTINUE C CLOSE(21) C STOP C END CC C SUBROUTINE COEFF(X,PX,QX,RX) CC CC Define the coefficient functions; here a Yukawa potential. CC C DOUBLE PRECISION X,PX,QX,RX, T CC CC Be careful with potential over/underflows; here we assume the CC IEEE double precision exponent range. CC C IF (X .LT. 650.0) THEN C T = EXP(-X) C ELSE C T = 0.0 C ENDIF C PX = 1.0 C QX = -T/X C RX = PX C RETURN C END CC CC End of sample driver for SLEDGE. C----------------------------------------------------------------------- C General remarks: C (1) Two machine dependent constants must be set in a DATA C statement in routine START (in part 4 of the package): C URN - an estimate of the unit roundoff; infinite output C values are assigned the value 1/URN. C UFLOW - a number somewhat smaller than -ln(underflow level). C Values of certain variables z for which C ln(abs(z)) < -under C will be set to zero. C (2) A value of IFLAG = -1, -2, or -3 may be the result of a C lack of smoothness in the coefficient functions. In such C cases a user input mesh may perform better (see (4) below). C (3) The heuristics for generating the initial mesh distribution C work reasonably well over a wide range of examples, but C occasionally they are far from optimal. The code's choice C can be over-ridden by setting JOB(5) False, setting NUMX C appropriately and supplying a mesh in XEF(*). C (4) If any of the coefficient functions p,q, or r (or their first C few derivatives) have finite jump discontinuities at points C in the interior of (A,B), then it is advantageous to have C these points in SLEDGE's mesh. Currently, this can only be C accomplished by setting JOB(5) False and supplying an C appropriate mesh using NUMX and XEF(*). C (5) In general, eigenvalue convergence is observed to be more C rapid than eigenfunction convergence; hence, it is C recommended that JOB(2) be False unless eigenfunction C information really is necessary. C (6) When eigenfunction output is sought, unless some knowledge C of the eigenfunction is known in advance, it is recommended C that JOB(5) be True so that the code will attempt to choose C a reasonable distribution for the initial mesh points. C (7) Computing the spectral density function for problems having C continuous spectrum can be very expensive; it is recommended C that initially, relatively crude tolerances (0.001 or so) be C used to get some idea of the effort required. C (8) It is recommended that every problem be classified (JOB(4) C False) by the code before any calculation of spectral C quantities occurs. Only if the user is certain as to what C the classification is (and describes it correctly through C INVEC and TYPE) should the classification option be bypassed. C (9) If the code does the classification of singular problems, it C will automatically choose the Friedrichs' boundary condition C at NONOSC endpoints. If another boundary condition is C desired, the user must use interval truncation in the C calling program (see remark (12)). C (10) While all parts of the code should function on machines C with a fairly narrow exponent range (such as IEEE single C precision), it is better to have a relatively wide exponent C range (IEEE double precision). The classification algorithm, C in particular, is far more reliable if done on a machine with C a fairly wide exponent range. C (11) Care must be taken in writing the subroutine COEFF for the C evaluation of p(x), q(x), and r(x) to avoid arithmetic C exceptions such as overflow and underflow (or trig function C arguments too large). This can be especially delicate on C machines with a small exponent range. C (12) In some cases `interval truncation' is recommended. By this C is meant the user should call SLEDGE several times using a C sequence of regular endpoints (with appropriate boundary C conditions) converging to the singular endpoint. The eigen- C values of the regular problems selected by the user should be C arranged so as to converge to those of the desired singular C problem. For example, if the user wishes to compute eigen- C values associated with a non-Friedrichs' boundary condition C for problems in spectral category 1, the user can experiment C with choosing a sequence of regular approximating intervals, C and vary the boundary conditions appropriately by means of a C `boundary condition function' or known solution of the C equation for a real value of EV on the sequence of regular C intervals until convergence of the regular eigenvalues to C the desired singular one is observed. Similarly, for C problems in spectral category 3 or 5 which involve one or two C endpoints which are LC and OSC, the (necessarily discrete) C spectrum is known to be unbounded below and above. To C implement a given LC boundary condition at a singular LC C endpoint one may choose a `boundary condition function' or C known solution of the equation for a real value of EV and C make use of it on a sequence of regular approximating C intervals to vary the boundary condition on successive calls C to SLEDGE for the sequence of regular intervals until C convergence to the desired singular eigenvalue is observed. C At present these methods are highly experimental and problem- C dependent as good heuristics for the choice of the rate of C convergence of the regular intervals to the singular one C which work well over a wide class of problems are not known. C (The only case in which SLEDGE automatically selects regular C approximating subintervals is for spectral density function C calculations for problems in spectral category 2; but C here the singular endpoint is of LP type, so no singular C boundary condition is required to be implemented.) C (13) Problems of slow convergence can sometimes be avoided by a C judicious change of either dependent or independent variable C (or both). C (14) If the Liouville normal form potential Q(t) has a minimum C far from zero, then the heuristics for generating the initial C mesh may well miss it. In this case, it is advisable to C shift the independent variable. C (15) The determination of the total number of eigenvalues is the C most difficult part of the classification process. When the C theory provides this number, of course, there is no problem; C otherwise, it should be viewed with some skepticism. A more C reliable count of the eigenvalues below the cutoff point of C the essential spectrum can be gained (at some expense) by C trying to compute many eigenvalues near that point. C/////////////////////////////////////////////////////////////////////// SUBROUTINE SLEDGE(JOB,CONS,ENDFIN,INVEC,TOL,TYPE,EV,NUMX,XEF,EF, & PDEF,T,RHO,IFLAG,STORE) INTEGER INVEC(*),NUMX,IFLAG(*) LOGICAL JOB(*),ENDFIN(*),TYPE(4,*) DOUBLE PRECISION CONS(*),TOL(*),EV(*),XEF(*),EF(*),PDEF(*),T(*), & RHO(*),STORE(*) C C This is the interface routine between the user and other routines C which carry out most of the actual calculations. C INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT, & KCLASS(2),I,IBASE,IEV,IPRINT,J,JTOL,K,KCL1,KCL2,LASTEV, & MAXITS,MAXT,MU1,MU2,NEV,NEXTRP,NUMEV,NUMT DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & CP(2),CR(2),CUTOFF,D(4,2),EMU(2),EP(2),EQLNF(2),ER(2), & ETA(2,2),PNU(2), & AA,ALPHA,BB,CEV(2),DENS,DENSHI,DENSLO,DENSOP,ENDFAC, & ENDI(5),ERROR,FZ,HMIN,RHOTOL,SGN,TOL1,TOLMAX,XTOL,ZETA, & ZETAI(5), & ZERO,HALF,ONE,TWO,FOUR LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2), & AAFIN,BBFIN,CSPEC(2),DOMESH,DONE,EDONE,JOBST(3), & LBASE,LMESH,LPLC(2),OSCILL COMMON /SLCLSS/CP,CR,CUTOFF,D,EMU,EP,EQLNF,ER,ETA,PNU,KCLASS COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, HALF = 0.5D0, ONE = 1.0, TWO = 2.0, & FOUR = 4.0, TOLMAX = 1.D-4) DATA DENSLO,DENSOP,DENSHI/4.0, 6.0, 12.0/ DATA ENDI/12.0, 20.0, 85.0, 240.0, 500.0/ DATA ZETAI/2.2, 2.0, 1.5, 1.4, 1.3/ C C Initialize. C AFIN = ENDFIN(1) BFIN = ENDFIN(2) IPRINT = INVEC(1) NUMT = INVEC(2) NEV = INVEC(3) LNF = .FALSE. DOMESH = .TRUE. LMESH = .FALSE. FLAG = 0 IFLAG(1) = 0 IBASE = 1 LBASE = .FALSE. DO 5 I = 1,6 LFLAG(I) = .FALSE. 5 CONTINUE TOL1 = MIN(TOL(1)+TOL(2),TOLMAX) JOBST(1) = JOB(2) IF ((NUMX .GT. 0) .AND. (.NOT. JOB(5))) THEN JOBST(2) = .TRUE. ELSE JOBST(2) = .FALSE. ENDIF JOBST(3) = JOB(3) IF ((.NOT. JOB(1)) .AND. (.NOT. JOB(2))) NEV = 0 CALL START(JOBST,CONS,TOL,NEV,INVEC(4),NUMX,XEF,NUMT,T, & NEXTRP,STORE) IF (JOB(4)) THEN ALPHA = A2*A1P-A1*A2P IF ((A1P .NE. ZERO) .OR. (A2P .NE. ZERO)) THEN IF (ALPHA .LE. ZERO) FLAG = -32 ELSE IF ((A1 .EQ. ZERO) .AND. (A2 .EQ. ZERO)) FLAG = -33 ENDIF IF ((B1 .EQ. ZERO) .AND. (B2 .EQ. ZERO)) FLAG = -31 ENDIF IF (FLAG .LT. 0) GOTO 120 IF (JOB(1) .OR. JOB(2)) THEN DO 10 K = 1,NEV EV(K) = ZERO 10 CONTINUE ENDIF IF (JOB(3)) THEN DO 15 K = 1,NUMT RHO(K) = ZERO 15 CONTINUE ENDIF IF ((.NOT. JOB(4)) .OR. JOB(3)) THEN CALL CLASS(IPRINT,TOL1,JOBST(2),CSPEC,CEV,LASTEV,LPLC,STORE, & JOB(5),HMIN,DOMESH) ALPHA = A2*A1P-A1*A2P IF ((A1P .NE. ZERO) .OR. (A2P .NE. ZERO)) THEN IF (ALPHA .LE. ZERO) FLAG = -32 ELSE IF ((A1 .EQ. ZERO) .AND. (A2 .EQ. ZERO)) FLAG = -33 ENDIF IF ((B1 .EQ. ZERO) .AND. (B2 .EQ. ZERO)) FLAG = -31 IF (FLAG .LT. 0) GOTO 120 DO 20 K = 1,2 TYPE(1,K) = REG(K) TYPE(2,K) = LC(K) TYPE(3,K) = .NOT. OSC(K) TYPE(4,K) = OSC(K) IF (CSPEC(K)) THEN TYPE(3,K) = .FALSE. TYPE(4,K) = .FALSE. ENDIF 20 CONTINUE IF (IPRINT .GT. 2) & CALL DSCRIP(LC,LPLC,TYPE,REG,CSPEC,CEV,CUTOFF,LASTEV, & A1,A1P,A2,A2P,B1,B2) ELSE LNF = .FALSE. KCLASS(1) = 0 KCLASS(2) = 0 DO 25 K = 1,2 REG(K) = TYPE(1,K) LC(K) = TYPE(2,K) OSC(K) = TYPE(4,K) CSPEC(K) = .NOT. (TYPE(3,K) .OR. TYPE(4,K)) 25 CONTINUE IF (.NOT. AFIN) STORE(1) = -99999.0 IF (.NOT. BFIN) STORE(NXINIT) = 99999.0 ENDIF C C Use NSGNF to hold the sign of F when EV is large negative. C SGN = A2P*B2 IF (SGN .NE. ZERO) THEN NSGNF = SIGN(ONE,SGN) ELSE SGN = A1P*B2+A2P*B1 IF (SGN .NE. ZERO) THEN NSGNF = SIGN(ONE,SGN) ELSE SGN = A1P*B1+A2*B2 IF (SGN .NE. ZERO) THEN NSGNF = SIGN(ONE,SGN) ELSE SGN = A1*B2+A2*B1 IF (SGN .NE. ZERO) THEN NSGNF = SIGN(ONE,SGN) ELSE NSGNF = SIGN(ONE,A1*B1) ENDIF ENDIF ENDIF ENDIF OSCILL = ((.NOT. TYPE(1,1)) .AND. TYPE(4,1)) .OR. & ((.NOT. TYPE(1,2)) .AND. TYPE(4,2)) TOL1 = TOL(1)+TOL(2) IF (JOB(1) .OR. JOB(2)) THEN C C Set up approximating regular problems for eigenvalues. C IF (OSCILL) THEN IF (IPRINT .GE. 1) THEN WRITE (*,30) WRITE (21,30) 30 FORMAT(' This problem is oscillatory, you must use ', & 'interval truncation.') ENDIF FLAG = -20 GOTO 120 ENDIF IF (DOMESH) THEN C C Calculate the initial mesh. C K = NXINIT+16 CALL MESH(JOB(5),-1,STORE,STORE(K),STORE(2*K+1), & STORE(3*K+1),STORE(4*K+1),TOL1,HMIN) IF (FLAG .LT. 0) GOTO 120 ENDIF IF (((KCLASS(1) .EQ. 3) .OR. (KCLASS(2) .EQ. 3)) .AND. JOB(5)) & LMESH = .TRUE. IF ((.NOT. LMESH) .AND. (IPRINT .GE. 1)) THEN WRITE (*,35) (STORE(I),I=1,NXINIT) WRITE (21,35) (STORE(I),I=1,NXINIT) 35 FORMAT(' Level 0 mesh:',/,(5G15.6)) ENDIF IF (JOB(5)) NUMX = NXINIT C C Set MAXLVL, the maximum number of levels (mesh bisections). C C IMPORTANT NOTE: the size of various fixed arrays in this C package depends on the value of MAXLVL in this FORTRAN77 C implementation. If MAXLVL is increased, then more storage C may have to be allocated to these arrays. In particular, C check RATIO(*), R(*,*), and W(*,*) in EXTRAP; EVEXT(*) C in REGULR. C MAXLVL = 10 C DO 45 K = 1,NEV EV(K) = ZERO IFLAG(K) = 0 FLAG = 0 CALL REGULR(JOB(2),LMESH,TOL,INVEC(3+K),EV(K),IPRINT, & NEXTRP,XEF,EF(1+NUMX*(K-1)), & PDEF(1+NUMX*(K-1)),HMIN,STORE) IF ((CSPEC(1) .OR. CSPEC(2)) .AND. (IPRINT .GE. 1) .AND. & (.NOT. JOB(4)) .AND. (FLAG .GT. -5)) THEN IF ((EV(K) .GE. CUTOFF) .OR. ((LASTEV .NE. -5) .AND. & (INVEC(3+K) .GE. LASTEV))) THEN WRITE (*,40) INVEC(3+K) WRITE (21,40) INVEC(3+K) 40 FORMAT(' WARNING: Requested eigenvalue ',I6, & ' may not be below the continuous spectrum.') ENDIF ENDIF IF (LFLAG(1)) THEN IFLAG(K) = IFLAG(K)+IBASE LFLAG(1) = .FALSE. LBASE = .TRUE. ENDIF IF (FLAG .LT. 0) IFLAG(K) = FLAG 45 CONTINUE IF (LBASE) THEN IBASE = 10*IBASE LBASE = .FALSE. ENDIF ENDIF IF (JOB(3)) THEN IF (CSPEC(1) .AND. CSPEC(2)) THEN IFLAG(1) = -5 IF (IPRINT .GT. 0) WRITE (*,50) 50 FORMAT(' This problem has continuous spectrum generated by', &' both endpoints. The',/,' calculation of the spectral density', &' function has not yet been implemented',/,' for such cases.',/) GOTO 120 ENDIF IF (.NOT. (CSPEC(1) .OR. CSPEC(2))) THEN IFLAG(1) = -4 IF (IPRINT .GT. 0) WRITE (*,55) 55 FORMAT(' This problem has no continuous spectrum.') GOTO 120 ENDIF IF ((CSPEC(1) .AND. ((KCLASS(2) .EQ. 5) .OR. & (KCLASS(2) .EQ. 9))) .OR. (CSPEC(2) .AND. & ((KCLASS(1) .EQ. 5) .OR. (KCLASS(1) .EQ. 9)))) THEN IFLAG(1) = -6 IF (IPRINT .GT. 0) WRITE (*,60) 60 FORMAT(' The normalization of the spectral density function' & ,' is unknown for this problem.') GOTO 120 ENDIF IF ((CSPEC(1) .AND. (.NOT. BFIN)) .OR. (CSPEC(2) .AND. & (.NOT. AFIN))) THEN IFLAG(1) = -6 IF (IPRINT .GT. 0) WRITE (*,60) GOTO 120 ENDIF IF (OSCILL) THEN FLAG = -25 IF (IPRINT .GT. 0) THEN WRITE (*,30) WRITE (21,30) ENDIF GOTO 120 ENDIF XTOL = -LOG10(MAX(TOL(1),TOL(2))) JTOL = XTOL-HALF JTOL = MIN(MAX(JTOL,1),5) DENSOP = 3*JTOL MAXITS = (15-JTOL)/3 C C Set Maxlvl for the density function calculation; see above C "IMPORTANT NOTE" if this is to be increased. C MAXLVL = (7+JTOL)/2 AAFIN = AFIN AA = A KCL1 = KCLASS(1) BBFIN = BFIN BB = B KCL2 = KCLASS(2) IF (JOB(5)) THEN C C Use interval truncation in this oscillatory regime. C OSCILL = .FALSE. IF ((.NOT. JOB(4)) .AND. ((KCLASS(1) .EQ. 1) .OR. & (KCLASS(2) .EQ. 1))) OSCILL = .TRUE. IF (.NOT. OSCILL) THEN NXINIT = 4*JTOL+5 ENDFAC = ENDI(JTOL) ELSE NXINIT = 24*JTOL+36 ENDFAC = 48.0 ENDIF IF (CSPEC(1)) THEN IF (AFIN) THEN KCLASS(1) = 7 ENDFAC = 4.0*ENDFAC IF (BFIN) THEN A = AA+(BB-AA)/ENDFAC ELSE A = AA+ABS(AA)/ENDFAC ENDIF ELSE AFIN = .TRUE. KCLASS(1) = 0 IF (BFIN) THEN A = -ENDFAC-MIN(-B,ZERO) ELSE A = -ENDFAC ENDIF ENDIF ELSE IF (BFIN) THEN KCLASS(2) = 7 ENDFAC = 4.0*ENDFAC IF (AFIN) THEN B = BB-(BB-AA)/ENDFAC ELSE B = BB-ABS(BB)/ENDFAC ENDIF ELSE BFIN = .TRUE. KCLASS(2) = 0 IF (AFIN) THEN B = ENDFAC+MAX(A,ZERO) ELSE B = ENDFAC ENDIF ENDIF ENDIF ELSE IF (CSPEC(1)) AFIN = .TRUE. IF (CSPEC(2)) BFIN = .TRUE. IF (NUMX .EQ. 0) THEN IFLAG(1) = -30 RETURN ENDIF ENDIF MAXT = NUMT C C Loop over the choices of intervals. C NUMEV = 0 DO 105 K = 1,MAXITS STORE(1) = A STORE(NXINIT) = B LFLAG(3) = .FALSE. FLAG = 0 IF (IPRINT .GE. 1) THEN WRITE (*,65) K WRITE (21,65) K 65 FORMAT(60('-'),/,' Iteration ',I2) WRITE (21,70) A,B,NXINIT 70 FORMAT(/,' For a, b =',2F15.8,/,' Nxinit = ',I4,/) ENDIF IF (JOB(5)) THEN I = NXINIT+16 CALL MESH(.TRUE.,-1,STORE,STORE(I+1),STORE(2*I+1), & STORE(3*I+1),STORE(4*I+1),TOL1,HMIN) ENDIF IF (IPRINT .GE. 3) THEN WRITE (*,75) (STORE(I),I=1,NXINIT) WRITE (21,75) (STORE(I),I=1,NXINIT) 75 FORMAT(' Level 0 mesh:',/,(5G15.6)) ENDIF CALL DENSEF(TOL,CSPEC,IPRINT,K,NEXTRP,MAXT,T,RHO,IEV, & HMIN,NUMEV,STORE) IF (FLAG .EQ. -3) THEN LFLAG(6) = .FALSE. FLAG = 0 ENDIF IF (FLAG .LT. 0) GOTO 120 IF (.NOT. JOB(5)) GOTO 110 IF (K .GT. 1) THEN DONE = .TRUE. J = MAXT DO 80 I = 1,J RHOTOL = TWO*ZETA*MAX(TOL(1),TOL(2)*RHO(I)) ERROR = RHO(I)-STORE(2320+(MAXLVL+2)*NUMT+I) IF (ABS(ERROR) .LE. RHOTOL) THEN EDONE = .TRUE. ELSE EDONE = .FALSE. MAXT = I ENDIF DONE = DONE .AND. EDONE 80 CONTINUE IF (DONE) GOTO 110 ENDIF IF (IPRINT .GE. 2) THEN WRITE (*,85) WRITE (21,85) 85 FORMAT(9X,'t',15X,'Truncated Rho(t)') DO 95 I = 1,NUMT WRITE (*,90) T(I),RHO(I) WRITE (21,90) T(I),RHO(I) 90 FORMAT(F12.4,D31.15) 95 CONTINUE ENDIF DO 100 I = 1,MAXT STORE(2320+(MAXLVL+2)*NUMT+I) = RHO(I) 100 CONTINUE COUNTZ = .TRUE. CALL SHOOT(CUTOFF,STORE,MU1,FZ) CALL SHOOT(T(NUMT),STORE,MU2,FZ) COUNTZ = .FALSE. IF (T(NUMT) .GT. CUTOFF) THEN DENS = (MU2-MU1)/(K*(T(NUMT)-CUTOFF)) ELSE DENS = DENSOP ENDIF IF (.NOT. OSCILL) THEN NXINIT = NXINIT+10 ZETA = ZETAI(JTOL) ELSE ZETA = 2.0 NXINIT = ZETA*NXINIT ENDIF IF (CSPEC(1)) THEN IF (AAFIN) THEN ENDFAC = 5.0*ZETA IF (DENS .LT. DENSLO) ENDFAC = 75.0 IF ((DENS .GT. DENSHI) .AND. (.NOT. OSCILL)) & ENDFAC = 8.0 A = AA+(A-AA)/ENDFAC IF ((AA-A)**2 .LT. U) THEN FLAG = -9 GOTO 110 ENDIF ELSE IF (DENS .LT. DENSLO) ZETA = 2.0 IF ((DENS .GT. DENSHI) .AND. (.NOT. OSCILL))ZETA = 1.4 A = ZETA*A ENDIF ENDIF IF (CSPEC(2)) THEN IF (BBFIN) THEN ENDFAC = 5.0*ZETA IF (DENS .LT. DENSLO) ENDFAC = 75.0 IF ((DENS .GT. DENSHI) .AND. (.NOT. OSCILL)) & ENDFAC = 8.0 B = BB-(BB-B)/ENDFAC IF ((B-BB)**2 .LT. U) THEN FLAG = -9 GOTO 110 ENDIF ELSE IF (DENS .LT. DENSLO) ZETA = 2.0 IF ((DENS .GT. DENSHI) .AND. (.NOT. OSCILL))ZETA = 1.4 B = ZETA*B ENDIF ENDIF IF (MOD(NXINIT,2) .EQ. 0) NXINIT = NXINIT+1 NXINIT = MIN(464,NXINIT) 105 CONTINUE FLAG = -3 110 IF (IPRINT .GE. 1) WRITE (21,115) NUMEV 115 FORMAT(' The total number of eigenvalues computed was ',I10) IF (CSPEC(1)) THEN IF (AAFIN) THEN A = AA KCLASS(1) = KCL1 ELSE AFIN = .FALSE. ENDIF ENDIF IF (CSPEC(2)) THEN IF (BBFIN) THEN B = BB KCLASS(2) = KCL2 ELSE BFIN = .FALSE. ENDIF ENDIF ENDIF C C Set fatal output flags. C 120 IF (FLAG .LT. -9) THEN DO 125 K = 1,MAX(NEV,1) IFLAG(K) = FLAG 125 CONTINUE RETURN ELSE IF ((FLAG .LT. 0) .AND. (.NOT. (JOB(1) .OR. JOB(2)))) & IFLAG(1) = FLAG ENDIF C C Set warning flags. C DO 135 I = 2,5 DO 130 K = 1,MAX(NEV,1) IF (LFLAG(I) .AND. (IFLAG(K) .GE. 0)) IFLAG(K) = & IFLAG(K)+I*IBASE 130 CONTINUE IF (LFLAG(I)) IBASE = 10*IBASE 135 CONTINUE RETURN END C======================================================================= SUBROUTINE INTERV(FIRST,ALPHA,BETA,CONS,ENDFIN,NFIRST,NTOTAL, & IFLAG,X) LOGICAL FIRST,ENDFIN(*) INTEGER NFIRST,NTOTAL,IFLAG DOUBLE PRECISION ALPHA,BETA,CONS(*),X(*) *********************************************************************** * * * INTERV calculates the indices of eigenvalues found in a * * specified interval. * * * *********************************************************************** C Local variables: C INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT, & IDUMMY(1),I1,I2,I3,J1,J2,J3,K,LASTEV,LEVEL0,MU,NEXTRP,NLAST,NUMX DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & CEV(2),CUTOFF,HMIN,TOL(6),V LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2), & CSPEC(2),DOMESH,JOBST(3),LPLC(2) COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER SAVE CUTOFF C NFIRST = -5 NLAST = -5 IFLAG = 0 IF (ALPHA .GE. BETA) THEN IFLAG = -9 RETURN ENDIF TOL(1) = 0.000001 TOL(2) = 0.000001 IF (FIRST) THEN JOBST(1) = .FALSE. JOBST(2) = .FALSE. JOBST(3) = .FALSE. AFIN = ENDFIN(1) BFIN = ENDFIN(2) MAXLVL = 10 NUMX = 0 CALL START(JOBST,CONS,TOL,0,IDUMMY,NUMX,X,0,X,NEXTRP,X) CALL CLASS(0,TOL(1),JOBST(2),CSPEC,CEV,LASTEV,LPLC,X, & .TRUE.,HMIN,DOMESH) CUTOFF = MIN(CEV(1),CEV(2)) IF (FLAG .LT. 0) THEN IFLAG = FLAG RETURN ENDIF IF (OSC(1) .OR. OSC(2)) THEN IFLAG = -25 RETURN ENDIF IF (DOMESH) THEN K = NUMX+16 CALL MESH(.TRUE.,-1,X,X(K),X(2*K+1),X(3*K+1),X(4*K+1), & TOL(1),HMIN) ENDIF ENDIF IF (FLAG .LT. 0) THEN IFLAG = FLAG RETURN ENDIF LEVEL0 = LEVEL COUNTZ = .TRUE. LEVEL = 3 CALL SHOOT(ALPHA,X,MU,V) I1 = MU CALL SHOOT(BETA,X,MU,V) J1 = MU LEVEL = LEVEL+1 CALL SHOOT(ALPHA,X,MU,V) I2 = MU CALL SHOOT(BETA,X,MU,V) J2 = MU 10 LEVEL = LEVEL+1 IF (NFIRST .EQ. -5) THEN CALL SHOOT(ALPHA,X,MU,V) I3 = MU IF ((I1 .EQ. I2) .AND. (I2 .EQ. I3)) THEN NFIRST = I1 GOTO 15 ENDIF I1 = I2 I2 = I3 ENDIF 15 IF (NLAST .EQ. -5) THEN CALL SHOOT(BETA,X,MU,V) J3 = MU IF ((J1 .EQ. J2) .AND. (J2 .EQ. J3)) THEN NLAST = J1-1 IF (NFIRST .NE. -5) GOTO 20 ENDIF J1 = J2 J2 = J3 ENDIF IF (LEVEL .LT. MAXLVL) GOTO 10 20 IF (NFIRST .EQ. -5) THEN NFIRST = I3 IFLAG = 12 ENDIF IF (NLAST .EQ. -5) THEN NLAST = J3 IFLAG = 12 ENDIF NTOTAL = NLAST+1-NFIRST IF (NTOTAL .EQ. 0) IFLAG = 11 IF (BETA .GT. CUTOFF) IFLAG = 13 COUNTZ = .FALSE. LEVEL = LEVEL0 RETURN END C/////////////////////////////////////////////////////////////////////// SUBROUTINE AITKEN(XLIM,TOL,N,X,ERROR) DOUBLE PRECISION XLIM,TOL,X(*),ERROR INTEGER N C C Use Aitken's algorithm to accelerate convergence of the sequence C in X(*). C DOUBLE PRECISION DENOM,XOLD,ZERO,ONE,TWO INTEGER I PARAMETER (ZERO = 0.0, ONE = 1.0, TWO = 2.0) C IF (N .LE. 2) THEN XLIM = X(N) ERROR = ZERO RETURN ENDIF XOLD = 1.D30 DO 10 I = 1,N-2 DENOM = X(I+2)-TWO*X(I+1)+X(I) IF (DENOM .NE. ZERO) THEN XLIM = X(I)-(X(I+1)-X(I))**2/DENOM ERROR = XLIM-XOLD ELSE ERROR = X(I+2)-X(I+1) XLIM = X(I+2) ENDIF IF (ABS(ERROR) .LT. MAX(ONE,ABS(XLIM))*TOL) RETURN XOLD = XLIM 10 CONTINUE RETURN END C---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION ASYMEV(NEV,QINT,RPINT,ALPHA1,ALPHA2, & BETA1,BETA2) INTEGER NEV DOUBLE PRECISION QINT,RPINT,ALPHA1,ALPHA2,BETA1,BETA2 C LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2) DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & FNEV, & ZERO,HALF,ONE,TWO,PI COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, HALF = 0.5D0, ONE = 1.0, TWO = 2.0, & PI = 3.14159265358979324D0) C C Evaluate the asymptotic formula for eigenvalue NEV. C Note: not all cases have been implemented yet. C ASYMEV = -999999.0 FNEV = NEV IF (REG(1)) THEN IF ((A1P .NE. ZERO) .OR. (A2P .NE. ZERO)) THEN IF (A2P .NE. ZERO) THEN ASYMEV = (TWO*A1P/A2P+QINT)/RPINT IF (B2 .NE. ZERO) THEN C Case 1 ASYMEV = ASYMEV+TWO*B1/(B2*RPINT)+((FNEV-ONE)*PI/ & RPINT)**2 ELSE C Case 2 ASYMEV = ASYMEV+((FNEV-HALF)*PI/RPINT)**2 ENDIF ELSE ASYMEV = (TWO*A2/A1P+QINT)/RPINT IF (B2 .NE. ZERO) THEN C Case 3 ASYMEV = ASYMEV+TWO*B1/(B2*RPINT)+((FNEV-HALF)*PI/ & RPINT)**2 ELSE C Case 4 ASYMEV = ASYMEV+(FNEV*PI/RPINT)**2 ENDIF ENDIF ELSE IF (A2 .NE. ZERO) THEN IF (B2 .NE. ZERO) THEN C Case 1 ASYMEV = (FNEV*PI/RPINT)**2+(TWO*(BETA1/BETA2+ & ALPHA1/ALPHA2)+QINT)/RPINT ELSE C Case 2 (Dirichlet at B) ASYMEV = ((FNEV+HALF)*PI/RPINT)**2+(TWO*ALPHA1/ALPHA2 & +QINT)/RPINT ENDIF ELSE IF (B2 .NE. ZERO) THEN C Case 3 (Dirichlet at A) ASYMEV = ((FNEV+HALF)*PI/RPINT)**2+(TWO*BETA1/BETA2 & +QINT)/RPINT ELSE C Case 4 (Dirichlet at A and at B) ASYMEV = ((FNEV+ONE)*PI/RPINT)**2+QINT/RPINT ENDIF ENDIF ENDIF RETURN ENDIF IF (REG(2)) THEN IF (B2 .NE. ZERO) THEN IF (A2 .NE. ZERO) THEN C Case 1 ASYMEV = (FNEV*PI/RPINT)**2+(TWO*(ALPHA1/ALPHA2+ & BETA1/BETA2)+QINT)/RPINT ELSE C Case 2 (Dirichlet at A) ASYMEV = ((FNEV+HALF)*PI/RPINT)**2+(TWO*BETA1/BETA2 & +QINT)/RPINT ENDIF ELSE IF (B2 .NE. ZERO) THEN C Case 3 (Dirichlet at B) ASYMEV = ((FNEV+HALF)*PI/RPINT)**2+(TWO*ALPHA1/ALPHA2 & +QINT)/RPINT ELSE C Case 4 (Dirichlet at A and at B) ASYMEV = ((FNEV+ONE)*PI/RPINT)**2+QINT/RPINT ENDIF ENDIF ENDIF RETURN END C---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION ASYMR(NEV,RPINT,RPATA,RPATB,SCALE) INTEGER NEV DOUBLE PRECISION RPINT,RPATA,RPATB,SCALE C LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2) DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & FNEV, & ZERO,HALF,ONE,TWO,PI COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, HALF = 0.5D0, ONE = 1.0, TWO = 2.0, & PI = 3.14159265358979324D0) C C Evaluate the asymptotic formula for RsubNEV. C Note: not all cases have been implemented yet. See the note C above in ASYMEV. C ASYMR = ZERO FNEV = MAX(NEV,2) IF (REG(1)) THEN IF ((A1P .NE. ZERO) .OR. (A2P .NE. ZERO)) THEN IF (A2P .NE. ZERO) THEN IF (B2 .NE. ZERO) THEN ASYMR = TWO*RPINT**3/(RPATA*A2P**2*((FNEV-ONE)*PI)**4) ELSE ASYMR = TWO*RPINT**3/(RPATA*A2P**2*((FNEV-HALF)*PI)**4) ENDIF ELSE IF (B2 .NE. ZERO) THEN ASYMR = RPATA*TWO*RPINT/(A1P*(FNEV-HALF)*PI)**2 ELSE ASYMR = RPATA*TWO*RPINT/(A1P*FNEV*PI)**2 ENDIF ENDIF ELSE IF (A2 .NE. ZERO) THEN ASYMR = TWO/(RPATA*A2*A2*RPINT) ELSE IF (B2 .NE. ZERO) THEN ASYMR = TWO*RPATA*((FNEV+HALF)*PI/A1)**2/RPINT**3 ELSE ASYMR = TWO*RPATA*((FNEV+ONE)*PI/A1)**2/RPINT**3 ENDIF ENDIF ENDIF RETURN ENDIF IF (REG(2)) THEN IF (A2 .NE. ZERO) THEN ASYMR = TWO/(RPATB*B2*B2*RPINT) ELSE IF (B2 .NE. ZERO) THEN ASYMR = TWO*RPATB*((FNEV+HALF)*PI/B1)**2/RPINT**3 ELSE ASYMR = TWO*RPATB*((FNEV+ONE)*PI/B1)**2/RPINT**3 ENDIF ENDIF ENDIF ASYMR = ASYMR*SCALE**2 RETURN END C----------------------------------------------------------------------- SUBROUTINE BRCKET(N,EVLOW,EVHIGH,FLOW,FHIGH,ABSERR,RELERR,X) INTEGER N DOUBLE PRECISION EVLOW,EVHIGH,FLOW,FHIGH,ABSERR,RELERR,X(*) C C Find values for EVLOW and EVHIGH which bracket the Nth eigenvalue; C in particular, C EV(N-1) < EVLOW < EV(N) < EVHIGH < EV(N+1) . C It is assumed that if U(X,LAMBDA) has NZ zeros in (A,B) then C EV(MU-1) < LAMBDA < EV(MU) C where MU is a function of NZ, LAMBDA, and the constants in the C boundary conditions. The value of MU for a given LAMBDA is C returned by a call to subprogram SHOOT. C INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT, & K,MU DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & DIFF,EV,EVSIGN,FEV, & ZERO,TWO LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2), & LOW,HIGH COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, TWO = 2.0) C C Set COUNTZ so that zeros are counted in SHOOT. C COUNTZ = .TRUE. EVSIGN = NSGNF C C SHOOT with Ev = Evlow should return FEV having sign EVSIGN. C IF (N .NE. 2*(N/2)) EVSIGN = -EVSIGN LOW = .FALSE. HIGH = .FALSE. C C Make EVLOW a lower bound for EV(N). C EV = EVLOW DIFF = ABS(EVHIGH-EVLOW) IF (DIFF .EQ. ZERO) DIFF = ABSERR+RELERR 10 CALL SHOOT(EV,X,MU,FEV) IF (FLAG .LT. 0) THEN COUNTZ = .FALSE. RETURN ENDIF IF (MU .GT. N) THEN EVHIGH = EV FHIGH = FEV EV = EV-DIFF DIFF = TWO*DIFF IF ((MU .EQ. N+1) .AND. (EVSIGN*FEV .LE. ZERO)) HIGH = .TRUE. GOTO 10 ELSE EVLOW = EV FLOW = FEV IF ((MU .EQ. N) .AND. (EVSIGN*FEV .GE. ZERO)) LOW = .TRUE. ENDIF C C Make EVHIGH an upper bound for EV(N). C IF (.NOT. HIGH) THEN EV = EVHIGH DIFF = ABS(EVHIGH-EVLOW) 20 CALL SHOOT(EV,X,MU,FEV) IF (FLAG .LT. 0) RETURN K = NSGNF*(-1)**MOD(MU,2) IF (K*FEV .LT. ZERO) THEN EV = EV+DIFF DIFF = TWO*DIFF GOTO 20 ELSE IF (MU .LE. N) THEN EVLOW = EV FLOW = FEV EV = EV+DIFF DIFF = TWO*DIFF GOTO 20 ELSE EVHIGH = EV FHIGH = FEV IF (MU .EQ. N+1) HIGH = .TRUE. ENDIF ENDIF ENDIF C C Refine the interval [EVLOW,EVHIGH] to include only the Nth C eigenvalue. C 30 IF ((.NOT. LOW) .OR. (.NOT. HIGH)) THEN DIFF = EVHIGH-EVLOW EV = EVLOW+DIFF/TWO C C Check for a cluster of eigenvalues within user's tolerance. C IF (TWO*DIFF .LT. MAX(ABSERR,RELERR*(MAX(ABS(EVLOW), & ABS(EVHIGH))))) THEN LFLAG(1) = .TRUE. COUNTZ = .FALSE. RETURN ENDIF CALL SHOOT(EV,X,MU,FEV) IF (FLAG .LT. 0) THEN COUNTZ = .FALSE. RETURN ENDIF C C Update EVLOW and EVHIGH. C IF (MU .EQ. N) THEN EVLOW = EV LOW = .TRUE. FLOW = FEV ELSE IF (MU .EQ. N+1) THEN EVHIGH = EV HIGH = .TRUE. FHIGH = FEV ELSE IF (MU .LT. N) THEN EVLOW = EV FLOW = FEV ELSE EVHIGH = EV FHIGH = FEV ENDIF ENDIF ENDIF GOTO 30 ENDIF COUNTZ = .FALSE. RETURN END C----------------------------------------------------------------------- SUBROUTINE CLASS(IPRINT,TOL,JOB,CSPEC,CEV,LASTEV,LPLC,X, & JMESH,HMIN,DOMESH) INTEGER IPRINT,LASTEV LOGICAL JOB,CSPEC(*),LPLC(*),JMESH,DOMESH DOUBLE PRECISION TOL,CEV(*),X(*),HMIN C C This routine classifies the Sturm-Liouville problem. Note: C (1) any computational algorithm must be based on a finite C amount of information; hence, there will always be cases that C any algorithm misclassifies. In addition, some problems are C inherently ill-conditioned, in that a small change in the C coefficients can produce a totally different classification. C (2) The maximum number of points sampled for singular problems C is given by the variable KMAX. By increasing this number, the C reliability of the classification may increase; however, the C computing time may also increase. The values we have chosen C seem to be a reasonable balance for most problems. C (3) The algorithms apply standard theorems involving limits of C the Liouville normal form potential. When this is not available, C each coefficient function is approximated by a power function C (c*x^r) and classified according to the properties of the C resulting Liouville approximation. C INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT, & KCLASS(2),IFLAG,J,K,KMAX,KUSED(2),LAST(3),M DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & CP(2),CR(2),CUTOFF,D(4,2),EMU(2),EP(2),EQLNF(2), & ER(2),ETA(2,2),PNU(2), & BASE,BC(2),END,EV,FEV,OVER,PZ(40,2),QZ(40,2), & RZ(40,2),S,SGN,Y(40),Z(40,2), & ZERO,TENTH,ONE,TWO,EIGHT LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2), & ENDFIN(2) COMMON /SLCLSS/CP,CR,CUTOFF,D,EMU,EP,EQLNF,ER,ETA,PNU,KCLASS COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, TENTH = 0.1D0, ONE = 1.0, TWO = 2.0, & EIGHT = 8.0) C C Sample the coefficient functions; determine if the problem as C given is in Liouville normal form. C LNF = .TRUE. DOMESH = .TRUE. DO 40 J = 1,2 IF (J .EQ .1) THEN END = A ENDFIN(1) = AFIN SGN = ONE ELSE END = B ENDFIN(2) = BFIN SGN = -ONE ENDIF K = TENTH-LOG10(TOL) KMAX = MIN(MAX(4*K,10),40) M = 0 BASE = ONE/MIN(MAX(K,4),8) IF (ENDFIN(J)) THEN IF (END .NE. ZERO) KMAX = MIN(-INT(LOG10(U))-1,KMAX) ELSE KMAX = MIN(KMAX,20) BASE = EIGHT ENDIF DO 10 K = 1,KMAX Z(K,J) = BASE**K 10 CONTINUE OVER = UNDER/TWO DO 30 K = 1,KMAX IF (ENDFIN(J)) THEN S = END+SGN*Z(K,J) ELSE S = -SGN*Z(K,J) ENDIF CALL COEFF(S,PZ(K,J),QZ(K,J),RZ(K,J)) NCOEFF = NCOEFF+1 IF ((PZ(K,J) .LE. ZERO) .OR. (RZ(K,J) .LE. ZERO)) THEN FLAG = -15 RETURN ENDIF IF (LNF .AND. (K .GT. 1)) THEN IF (PZ(K,J) .NE. PZ(K-1,J)) LNF = .FALSE. IF (RZ(K,J) .NE. RZ(K-1,J)) LNF = .FALSE. ENDIF S = LOG(PZ(K,J)) IF (ABS(S) .GT. OVER) M = 1 S = LOG(ONE+ABS(QZ(K,J))) IF (ABS(S) .GT. OVER) M = 1 S = LOG(RZ(K,J)) IF (ABS(S) .GT. OVER) M = 1 IF (M .NE. 0) GOTO 35 30 CONTINUE K = KMAX 35 KUSED(J) = K-1 40 CONTINUE DO 50 J = 1,2 CALL CLSEND(Z(1,J),PZ(1,J),QZ(1,J),RZ(1,J),KUSED(J),IPRINT, & END,ENDFIN(J),J,TOL,CEV(J),CSPEC(J),BC,Y,LPLC,IFLAG) IF (.NOT. JOB) THEN IF ((.NOT. AFIN) .AND. (J .EQ. 1)) X(1) = -END IF ((.NOT. BFIN) .AND. (J .EQ. 2)) X(NXINIT) = END ENDIF IF ((CSPEC(J) .OR. (.NOT. OSC(J))) .AND. (.NOT. REG(J))) THEN IF (J .EQ. 1) THEN A1 = BC(1) A1P = ZERO A2 = BC(2) A2P = ZERO ELSE B1 = BC(1) B2 = BC(2) ENDIF ENDIF IF (IFLAG .EQ. 1) LFLAG(2) = .TRUE. 50 CONTINUE CUTOFF = MIN(CEV(1),CEV(2)) C C Find the number of eigenvalues below the start of the C continuous spectrum. C LASTEV = -5 IF ((.NOT. OSC(1)) .AND. (.NOT. LC(2)) .AND. OSC(2)) LASTEV = 0 IF ((.NOT. OSC(2)) .AND. (.NOT. LC(1)) .AND. OSC(1)) LASTEV = 0 IF (OSC(1) .AND. (.NOT. LC(1)) .AND. OSC(2) .AND. LC(2))LASTEV = 0 IF (OSC(2) .AND. (.NOT. LC(2)) .AND. OSC(1) .AND. LC(1))LASTEV = 0 IF ((CSPEC(1) .AND. OSC(2)) .OR. (CSPEC(2) .AND. OSC(1))) LASTEV=0 IF ((CSPEC(1) .AND. (.NOT. OSC(2))) .OR. & (CSPEC(2) .AND. (.NOT. OSC(1))) .OR. & (CSPEC(1) .AND. CSPEC(2))) THEN K = NXINIT+16 CALL MESH(JMESH,-1,X,X(K),X(2*K+1),X(3*K+1),X(4*K+1),TOL,HMIN) DOMESH = .FALSE. COUNTZ = .TRUE. DO 75 J = 1,2 LEVEL = 3*J EV = CUTOFF CALL SHOOT(EV,X,LAST(J),FEV) IF (FLAG .LT. 0) RETURN IF (IPRINT .GE. 5) WRITE (21,70) LEVEL,LAST(J) 70 FORMAT(' When level = ',I2,', Ev index at cutoff is ',I12) 75 CONTINUE COUNTZ = .FALSE. IF (LAST(1) .GE. LAST(2)) THEN LASTEV = LAST(2) IF (LAST(1) .NE. LAST(2)) THEN LFLAG(2) = .TRUE. IF (IPRINT .GE. 3) WRITE (21,80) 80 FORMAT(' The eigenvalue count is uncertain.') IF (LAST(1) .GT. 2*LAST(2)) LASTEV = 0 ENDIF ELSE LASTEV = -5 ENDIF ENDIF RETURN END C----------------------------------------------------------------------- SUBROUTINE CLSEND(Z,PZ,QZ,RZ,KMAX,IPRINT,END,ENDFIN,IEND,TOL,CEV, & CSPEC,BC,Y,LPLC,IFLAG) INTEGER KMAX,IPRINT,IEND,IFLAG LOGICAL ENDFIN,CSPEC,LPLC(*) DOUBLE PRECISION Z(*),PZ(*),QZ(*),RZ(*),END,TOL,CEV,BC(*),Y(*) C C Iflag = 0 if reasonably certain of the classification; C = 1 if not sure. C C Information about the nature of the problem at singular point C IEND is passed through the variable KCLASS(IEND): C KCLASS(*) = 0 normal; C = 1 oscillatory coefficient function; C = 2 regular, but 1/p, q, or r unbounded; C = 3 infinite endpoint, Eqlnf = -1 ; C = 4 finite singular endpoint, Tau unbounded, C (not 8-10); C = 5 not "hard", irregular; C = 6 "hard" irregular with Eta(1) < 0; C = 7 finite end which generates Cspectrum; C = 8 Q is unbounded (< 1/t^2) near a nonoscill- C atory finite end; C = 9 Q is unbounded (like 1/t^2) near a nonosc- C illatory finite end; C = 10 "hard", irregular, Eta(1) > 0. C Note: "hard" means Tau goes to +infinity at a C finite nonoscillatory endpoint. C REG(*) = .True. iff endpoint is regular. C LC(*) = .True. iff endpoint is limit circle. C OSC(*) = .True. iff endpoint is oscillatory for all Ev. C CSPEC = .True. iff endpoint generates continuous spectrum. C LPLC(*)= .True. iff theory yields Lp/Lc classification. C INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT, & KCLASS(2),I,IQLNF,K LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2), & EX,EXACT,IRREG,POSC,QOSC,ROSC DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & CPT(2),CRT(2),CUTOFF,D(4,2),EMU(2),EPT(2),EQLNF(2), & ERT(2),ETA(2,2),PNU(2), & CP,CQ,CR,C1,C2,C3,DELTA,EP,EQ,ER,GAMMA,SGN,TOL4,ZZ, & ZERO,QUART,HALF,QUART3,ONE,TWO,FOUR COMMON /SLCLSS/CPT,CRT,CUTOFF,D,EMU,EPT,EQLNF,ERT,ETA,PNU,KCLASS COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, QUART = 0.25D0, HALF = 0.5D0, & QUART3 = 0.75D0, ONE = 1.0, TWO = 2.0, FOUR = 4.0) C IFLAG = 0 IF (IPRINT .GE. 3) THEN IF (IEND .EQ. 1) WRITE (21,5) 5 FORMAT(/,' For endpoint A') IF (IEND .EQ. 2) WRITE (21,6) 6 FORMAT(/,' For endpoint B:') WRITE (21,7) KMAX 7 FORMAT(' Kmax = ',I3) ENDIF CSPEC = .FALSE. IRREG = .FALSE. LPLC(IEND) = .TRUE. KCLASS(IEND) = 0 PNU(IEND) = ZERO CEV = ONE/U EX = .TRUE. C1 = ZERO C2 = ZERO TOL4 = FOUR*TOL C C Seek monomial approximations to each coefficient function. C CALL POWER(Z,QZ,KMAX,TOL,IPRINT,EQ,CQ,QOSC,EXACT,Y,IFLAG) IF (ABS(CQ) .LE. TOL) CQ = ZERO IF (CQ .EQ. ZERO) EQ = ZERO IF (LNF) THEN EP = ZERO CP = PZ(1) ER = ZERO CR = RZ(1) POSC = .FALSE. ROSC = .FALSE. ELSE CALL POWER(Z,PZ,KMAX,TOL,IPRINT,EP,CP,POSC,EXACT,Y,IFLAG) IF (ABS(CP) .LE. TOL) EP = ZERO EX = EX .AND. EXACT CALL POWER(Z,RZ,KMAX,TOL,IPRINT,ER,CR,ROSC,EXACT,Y,IFLAG) IF (ABS(CR) .LE. TOL) ER = ZERO ENDIF IF (POSC .OR. ROSC) THEN IF (ENDFIN) THEN REG(IEND) = .TRUE. ELSE IFLAG = 1 IF (IPRINT .GE. 3) WRITE (21,10) 10 FORMAT(' WARNING: p(x) or r(x) is not well-approximated by' & ,' a power potential.',/,' Classification is uncertain.') REG(IEND) = .FALSE. KCLASS(IEND) = 1 ENDIF LC(IEND) = .TRUE. OSC(IEND) = .FALSE. ENDIF IF (QOSC) THEN IFLAG = 1 IF (IPRINT .GE. 3) WRITE (21,20) 20 FORMAT(' WARNING: q(x) is not well-approximated by a power ', & 'potential.',/,' Classification is uncertain.') IF (ENDFIN) THEN REG(IEND) = .TRUE. LC(IEND) = .TRUE. OSC(IEND) = .FALSE. ELSE KCLASS(IEND) = 1 REG(IEND) = .FALSE. LC(IEND) = .FALSE. CSPEC = .TRUE. OSC(IEND) = .FALSE. BC(1) = ONE BC(2) = ZERO CEV = QZ(KMAX-1) K = 40 DELTA = (Z(KMAX)-Z(KMAX-1))/(K+1) DO 30 I = 0,K ZZ = Z(KMAX)-I*DELTA CALL COEFF(ZZ,CP,CQ,CR) NCOEFF = NCOEFF+1 CEV = MIN(CEV,CQ) 30 CONTINUE IF (ABS(CEV) .LT. TOL4) CEV = ZERO IF (U*ABS(CEV) .GE. ONE) CSPEC = .FALSE. ENDIF EQLNF(IEND) = ZERO ENDIF IF (IPRINT .GE. 3) THEN WRITE (21,40) CP,EP,CQ,EQ,CR,ER 40 FORMAT(' Cp, Ep; Cq, Eq; Cr, Er =',/,3(2D25.12,/)) ENDIF CPT(IEND) = CP EPT(IEND) = EP C C Analyze this endpoint. C IF ((EP .LT. ONE) .AND. (EQ .GT. -ONE) .AND. & (ER .GT. -ONE) .AND. ENDFIN) THEN REG(IEND) = .TRUE. LC(IEND) = .TRUE. OSC(IEND) = .FALSE. CSPEC = .FALSE. IF ((EP .GT. ZERO) .OR. (EQ .LT. ZERO) .OR. (ER .LT. ZERO)) & KCLASS(IEND) = 2 RETURN ENDIF REG(IEND) = .FALSE. ETA(1,IEND) = HALF*(ER-EP+TWO) IF (ABS(ETA(1,IEND)) .LE. TOL) ETA(1,IEND) = ZERO IF (ETA(1,IEND) .NE. ZERO) THEN EQLNF(IEND) = (EQ-ER)/ETA(1,IEND) IQLNF = EQLNF(IEND)+SIGN(HALF,EQLNF(IEND)) IF (ABS(IQLNF-EQLNF(IEND)) .LT. TOL4) EQLNF(IEND) = IQLNF C1 = (CQ/CR)*(ABS(ETA(1,IEND))*SQRT(CP/CR))**EQLNF(IEND) IF (C1 .EQ. ZERO) EQLNF(IEND) = ZERO C2 = (EP+ER)/(FOUR*ETA(1,IEND)) C2 = C2*(C2-ONE) IF (IPRINT .GE. 5) THEN WRITE (21,50) C1,C2,EQLNF(IEND),ETA(1,IEND) 50 FORMAT(' C1, C2 =',2D20.10,/,' Eqlnf, Eta1 =',2D20.10) ENDIF ELSE C3 = (CP/CR)*(QUART*(EP+ER))**2 IF (IPRINT .GE. 5) THEN WRITE (21,60) C3 60 FORMAT(' C3 = ',D19.10) ENDIF ENDIF IF (.NOT. ENDFIN) THEN C C Make an initial estimate for "infinity" (used in MESH). C IF ((ER .GT. EQ) .OR. (CQ .EQ. ZERO)) THEN IF (ER .NE. EP) THEN GAMMA = ER-EP DELTA = CR/CP ELSE GAMMA = EQ-EP DELTA = ABS(CQ)/CP IF (DELTA .EQ. ZERO) DELTA = ONE ENDIF ELSE IF (EQ .GT. ER) THEN GAMMA = EQ-EP DELTA = ABS(CQ)/CP ELSE IF (ER .GT. EP) THEN GAMMA = ER-EP DELTA = CR/CP ELSE GAMMA = ZERO IF (ER .EQ. EP) THEN DELTA = ABS(CR-CQ)/CP ELSE DELTA = ABS(CQ)/CP ENDIF ENDIF ENDIF ENDIF IF (GAMMA .GT. HALF) THEN IF (GAMMA .LT. TWO) THEN END = 80.0 ELSE END = MIN(MAX(64.0/((TWO*GAMMA-3.0)*DELTA**(ONE/ & (GAMMA+TWO))),ONE),80.D0) IF (GAMMA .GT. 24.0) END = 12.0 ENDIF ELSE IF (GAMMA .LT. -HALF) THEN END = MAX(MIN(600.0*DELTA**(ONE/GAMMA)*5.0**GAMMA, & 120.0D0),ONE) ELSE END = 12.0 ENDIF IF ((GAMMA .EQ. ZERO) .AND. (CQ .NE. ZERO)) END = 40.0 ENDIF ENDIF C C Test for finite irregular singular points. C IF (ENDFIN) THEN SGN = ONE I = ER-EP+SIGN(HALF,ER-EP) K = EQ-EP+SIGN(HALF,EQ-EP) IRREG = .TRUE. IF (CQ .EQ. ZERO) THEN IF ((I .GE. -2) .AND. (ABS(ER-EP-I) .LE. TOL4)) & IRREG = .FALSE. ELSE IF ((ER .LE. EQ) .AND. (I .GE. -2) .AND. & (ABS(ER-EP-I) .LE. TOL4)) IRREG = .FALSE. IF ((ER .GT. EQ) .AND. (K .GE. -2) .AND. & (ABS(EQ-EP-K) .LE. TOL4)) IRREG = .FALSE. ENDIF EMU(IEND) = HALF*(ONE-EP) IF (IRREG) THEN IF (IPRINT .GE. 3) WRITE (21,70) 70 FORMAT(' This is an irregular singular point.') KCLASS(IEND) = 5 ELSE C C Compute the principal Frobenius root. C IF (ETA(1,IEND) .NE. ZERO) THEN IF ((CQ .NE. ZERO) .AND. (ER .GT. EQ) .AND. (K .EQ. -2)) & THEN PNU(IEND) = EMU(IEND)**2+CQ/CP IF (ABS(PNU(IEND)) .LE. TOL4) PNU(IEND) = ZERO IF (PNU(IEND) .GE. ZERO) THEN PNU(IEND) = EMU(IEND)+SQRT(PNU(IEND)) ELSE PNU(IEND) = -EP ENDIF ELSE PNU(IEND) = MAX(ONE-EP,ZERO) ENDIF IF (PNU(IEND) .GT. -EP) THEN IF (IPRINT .GE. 5) WRITE (21,75) PNU(IEND) 75 FORMAT(' The principal Frobenius root is ',E20.8) ENDIF ENDIF ENDIF ELSE SGN = -ONE ENDIF IF (SGN*ETA(1,IEND) .GT. ZERO) THEN C C Carry out the Case 1 tests. C K = 0 IF (EQLNF(IEND) .LT. -TWO) THEN IF (CQ .LT. ZERO) K = 1 IF (CQ .GT. ZERO) K = -1 ENDIF IF (EQLNF(IEND) .EQ. -TWO) THEN IF (ABS(C1+C2+QUART) .LE. TOL4) THEN IF (IPRINT .GE. 3) WRITE (21,80) 80 FORMAT(' WARNING: borderline nonoscillatory/oscillato', & 'ry classification.') K = -1 IFLAG = 1 ELSE IF (C1+C2 .LT. -QUART-TOL4) K = 1 IF (C1+C2 .GT. -QUART) K = -1 ENDIF ENDIF IF (EQLNF(IEND) .GT. -TWO) THEN IF (ABS(C2+QUART) .LE. TOL4) THEN C2 = -QUART IF (IPRINT .GE. 3) WRITE (21,80) IFLAG = 1 ENDIF IF (C2 .GE. -QUART) K = -1 ENDIF IF (K .EQ. 1) THEN OSC(IEND) = .TRUE. ELSE IF (K .EQ. -1) THEN OSC(IEND) = .FALSE. ELSE IF (IPRINT .GE. 3) WRITE (21,85) 85 FORMAT(' NO INFORMATION on osc/nonosc class.') ENDIF ENDIF K = 0 IF (EQLNF(IEND) .LT. -TWO) THEN IF (CQ .GT. ZERO) K = -1 IF (CQ .LT. ZERO) K = 1 ENDIF IF (EQLNF(IEND) .EQ. -TWO) THEN IF (ABS(C1+C2-QUART3) .LE. TOL4) THEN K = -1 IF (IPRINT .GE. 3) WRITE (21,90) 90 FORMAT(' WARNING: borderline Lc/Lp classification.') IFLAG = 1 ENDIF IF (C1+C2 .GE. QUART3) K = -1 IF (ABS(C1+C2) .LT. QUART3-TOL4) K = 1 IF (C1+C2 .LT. -TOL4) K = 1 ENDIF IF (EQLNF(IEND) .GT. -TWO) THEN IF (ABS(C2-QUART3) .LE. TOL4) THEN K = -1 IF (IPRINT .GE. 3) WRITE (21,90) IFLAG = 1 ENDIF IF (C2 .GE. QUART3) K = -1 IF (ABS(C2) .LT. QUART3-TOL4) K = 1 IF (C2 .LT. -TOL4) K = 1 ENDIF IF (K .EQ. 1) THEN LC(IEND) = .TRUE. ELSE IF (K .EQ. -1) THEN LC(IEND) = .FALSE. ELSE WRITE (21,95) 95 FORMAT(' NO INFORMATION on Lp/Lc class.') ENDIF ENDIF ENDIF IF (SGN*ETA(1,IEND) .LT. ZERO) THEN C C Carry out the Case 2 tests. C K = 0 IF ((EQLNF(IEND) .GT. ZERO) .AND. (CQ .LT. ZERO)) K = 1 IF ((EQLNF(IEND) .GT. ZERO) .AND. (CQ .GT. ZERO)) K = -1 IF (EQLNF(IEND) .EQ. ZERO) THEN K = -1 CEV = CQ/CR CSPEC = .TRUE. IF (U*ABS(CEV) .GE. ONE) CSPEC = .FALSE. ENDIF IF (EQLNF(IEND) .LT. ZERO) THEN K = -1 CEV = ZERO CSPEC = .TRUE. ENDIF IF (K .EQ. 1) THEN OSC(IEND) = .TRUE. ELSE IF (K .EQ. -1) THEN OSC(IEND) = .FALSE. ELSE WRITE (21,100) 100 FORMAT(' NO INFORMATION on Osc/Nonosc class.') ENDIF ENDIF K = 0 IF ((EQLNF(IEND) .GT. TWO) .AND. (CQ .GT. ZERO)) K = -1 IF (EQLNF(IEND) .LE. TWO) K = -1 IF ((EQLNF(IEND) .GT. TWO) .AND. (CQ .LT. ZERO)) K = 1 IF (K .EQ. 1) THEN LC(IEND) = .TRUE. ELSE IF (K .EQ. -1) THEN LC(IEND) = .FALSE. ELSE WRITE (21,105) 105 FORMAT(' NO INFORMATION on Lp/Lc class.') ENDIF ENDIF ENDIF IF (ETA(1,IEND) .EQ. ZERO) THEN C C Carry out the Case 3 and 4 tests. C IF ((SGN*(EQ-ER) .LT. ZERO) .AND. (CQ .LT. ZERO)) THEN OSC(IEND) = .TRUE. LC(IEND) = .TRUE. ENDIF IF ((SGN*(EQ-ER) .LT. ZERO) .AND. (CQ .GT. ZERO)) THEN OSC(IEND) = .FALSE. LC(IEND) = .FALSE. ENDIF IF (EQ .EQ. ER) THEN OSC(IEND) = .FALSE. LC(IEND) = .FALSE. CEV = CQ/CR+C3 CSPEC = .TRUE. IF (U*ABS(CEV) .GE. ONE) CSPEC = .FALSE. ENDIF IF ((SGN*(EQ-ER) .GT. ZERO) .OR. (CQ .EQ. ZERO)) THEN OSC(IEND) = .FALSE. LC(IEND) = .FALSE. CEV = C3 CSPEC = .TRUE. IF (U*ABS(CEV) .GE. ONE) CSPEC = .FALSE. ENDIF ENDIF IF (ABS(CEV) .LE. TOL4) CEV = ZERO C C Calculate the Friedrichs boundary condition (if appropriate). C IF (CSPEC) THEN BC(1) = ONE BC(2) = ZERO ENDIF IF ((.NOT. CSPEC) .AND. (.NOT. OSC(IEND))) THEN IF ((SGN*(ER+EP) .GT. ZERO) .AND. (SGN*(EQ+EP) .GT. ZERO)) THEN BC(1) = ZERO BC(2) = ONE ELSE IF ((SGN*(ER+EP) .GT. ZERO) .AND. (EQ+EP .EQ. ZERO)) THEN BC(1) = SQRT(CP*ABS(CQ)) BC(2) = ONE IF (BC(1) .GT. ONE) THEN BC(2) = ONE/BC(1) BC(1) = ONE ENDIF ELSE IF ((SGN*(ER+EP) .LT. ZERO) .OR. & (SGN*(EQ+EP) .LT. ZERO)) THEN BC(1) = ONE BC(2) = ZERO ENDIF ENDIF ENDIF ENDIF IF (.NOT. OSC(IEND)) THEN IF ((.NOT. ENDFIN) .AND. (EQLNF(IEND) .EQ. -ONE)) & KCLASS(IEND) = 3 I = ER-EP IF (CQ .NE. ZERO) THEN K = EQ-EP IF (CQ .GT. ZERO) THEN IF (K .LT. I) I = 0 ELSE I = MIN(I,K) ENDIF ENDIF IF (ENDFIN .AND. (I .LT. 0)) THEN IF (IRREG) KCLASS(IEND) = 6 IF (ETA(1,IEND) .GT. ZERO) THEN C C Transform some nonoscillatory problems for which Tau C is unbounded near a finite singular endpoint. C IF (IRREG) KCLASS(IEND) = 10 CPT(IEND) = CP CRT(IEND) = CR EPT(IEND) = EP ERT(IEND) = ER EMU(IEND) = ZERO D(1,IEND) = (ETA(1,IEND)*SQRT(CP/CR))** & (ONE/ETA(1,IEND)) D(2,IEND) = ONE/SQRT(SQRT(CP*CR*D(1,IEND)**(EP+ER))) IF ((EQLNF(IEND) .EQ. -TWO) .OR. (C1 .EQ. ZERO)) THEN IF (.NOT. IRREG) KCLASS(IEND) = 9 EMU(IEND) = ABS(QUART+C1+C2) IF (EMU(IEND) .LT. TOL4) THEN EMU(IEND) = HALF ELSE EMU(IEND) = HALF+SQRT(EMU(IEND)) ENDIF ELSE IF (.NOT. IRREG) KCLASS(IEND) = 8 ENDIF ETA(2,IEND) = EMU(IEND)-QUART*(EP+ER)/ETA(1,IEND) IF ((KCLASS(IEND) .EQ. 10) .AND. (EMU(IEND) .EQ. ZERO)) & THEN ETA(2,IEND) = HALF*(ONE-EP)/ETA(1,IEND) EMU(IEND) = ETA(2,IEND)+QUART*(EP+ER)/ETA(1,IEND) ENDIF D(3,IEND) = ETA(2,IEND)*(ETA(2,IEND)+(EP-ONE)/ & ETA(1,IEND)) D(4,IEND) = D(3,IEND) IF (EQLNF(IEND) .EQ. -TWO) D(4,IEND) = D(4,IEND)-C1 IF (ABS(D(4,IEND)) .LE. TOL4) D(4,IEND) = ZERO D(4,IEND) = SQRT(ABS(D(4,IEND))) IF (IPRINT .GE. 5) THEN WRITE (21,110) EMU(IEND),ETA(2,IEND) WRITE (21,115) D(3,IEND),D(4,IEND) 110 FORMAT(' Mu =',D20.6,'; Eta2 =',D20.6) 115 FORMAT(' D3 =',D20.6,'; D4 =',D20.6) ENDIF ENDIF IF (KCLASS(IEND) .GE. 9) THEN IF (.NOT. EX) LFLAG(5) = .TRUE. IF (IPRINT .GE. 3) THEN WRITE (21,120) 120 FORMAT(' This problem has unbounded ', & '[Ev*r(x)-q(x)]/p(x).') IF ((EMU(IEND) .GT. ZERO) .OR. (EP+ER .NE. ZERO)) & WRITE (21,125) 125 FORMAT(' A change of variables will be used near', & ' this endpoint.') ENDIF ENDIF ENDIF IF (ENDFIN .AND. (.NOT. REG(IEND)) .AND. (KCLASS(IEND) & .EQ. 0)) KCLASS(IEND) = 4 ENDIF IF ((POSC .OR. QOSC .OR. ROSC) .AND. (.NOT. ENDFIN)) END = 99.0 IF (IPRINT .GE. 5) THEN WRITE (21,130) KCLASS(IEND) 130 FORMAT(' Classification type (KCLASS) is: ',I2) ENDIF RETURN END C----------------------------------------------------------------------- SUBROUTINE DENSEF(TOL,CSPEC,IPRINT,ITER,NEXTRP,NUMT,T,RHO,NEV, & HMIN,NUMEV,STORE) INTEGER IPRINT,ITER,NEXTRP,NUMT,NEV,NUMEV LOGICAL CSPEC(2) DOUBLE PRECISION TOL(*),T(*),RHO(*),HMIN,STORE(*) *********************************************************************** * * * This routine computes the spectral density function rho(t). * * * *********************************************************************** C C Input parameters: C TOL(*) as in SLEDGE. C CSPEC(*) = logical 2-vector; CSPEC(i) = .true. iff endpoint i C (1 = A, 2 = B) generates continuous spectrum. C IPRINT = integer controlling printing. C ITER = iteration from SLEDGE. C NEXTRP = integer giving maximum no. of extrapolations. C NUMT = integer equalling number of T(*) points. C T(*) = real vector of abcissae for spectral function rho(t). C HMIN = minimum stepsize in Level 0 mesh. C C Output parameters: C RHO(*) = real vector of values for spectral density function C rho(t), RHO(I) = rho(T(I)). C NEV = integer pointer to eigenvalue. On a normal return C (FLAG = 0) this is set to the index of the last C eigenvalue computed; if FLAG is not zero, then NEV C gives the index of the eigenvalue where the problem C occurred. C NUMEV = cumulative number of eigenvalues computed. C C Auxiliary storage: C STORE(*) = real vector of auxiliary storage, must be dimensioned C at least 5*Nxinit+(Maxlvl+2)*NUMT. The value of C Nxinit is either the input NUMX or Maxint. Currently, C Maxlvl = 8 and Maxint = 235. C 1 -> Nxinit vector of mesh points X(*), C Nxinit+1 -> 5*Nxinit intermediate RsubN calculations, C 5*Nxinit+1 -> 5*Nxinit+10*NUMT intermediate RHO values. C----------------------------------------------------------------------- C The definition of a spectral density function assumes a certain C normalization on the eigenfunctions. For the case when x = b C generates the continuous spectrum, the normalization used here C (and in routine GETRN below) is: C (1) when x = a is regular u(a) = (A2-A2P*Ev)/SCALE C (pu')(a) = (A1-A1P*Ev)/SCALE, C with SCALE = sqrt(A1**2+A2**2) when A1' = A2' = 0, and C SCALE = sqrt(ALPHA) otherwise. C (2) When x = a is a regular singular point then u(x) is taken to C be asymptotic to the principal Frobenius solution, i.e., C near x = a u(x) ~ (x-a)**Nu with Nu the larger C root of the indicial equation. C Analogous normalizations hold at x = b when the endpoint x = a C generates the continuous spectrum. C---------------------------------------------------------------------- C Local variables: C INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT, & KCLASS(2),I,IASYMP,J,JS,KLVL,KRHO,LPRINT,MAXNEV,NRHO,NSAVE LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2), & DONE,EFIN(2,2),JUMP,RDONE DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER, & CP(2),CR(2),CUTOFF,D(4,2),EMU(2),EP(2),EQLNF(2),ER(2), & ETAT(2,2),PNU(2), & ABSERR,ALPHA,ALPHA1,ALPHA2,ASYMEV,ASYMR,AEV,ARN,AVGRHO, & BETA1,BETA2,BIG,DELTA,DENOM,DX,EFNORM,ERROR,ETA,ETAOLD,EV, & EVHAT,EVHIGH,EVLOW,EVOLD(200),EVSAVE,FLOW,FHIGH,H,HALFH, & PDU,PSRHO,PX,QINT,QLNF,QX,RELERR,RHOSUM,ROLD,RPATA,RPATB, & RPINT,RSAVE(5),RSUBN,RX,SCALE,SQRTRP,TENU,TOL1,TOL2, & TOLMAX,TOLMIN,UX,XLEFT,XTOL,Z,ZABS,ZREL, & ZERO,HALF,ONE,TWO,THREE,FOUR,SIX,EIGHT,TEN COMMON /SLCLSS/CP,CR,CUTOFF,D,EMU,EP,EQLNF,ER,ETAT,PNU,KCLASS COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER PARAMETER (ZERO = 0.0, HALF = 0.5D0, ONE = 1.0, & TWO = 2.0, THREE = 3.0, FOUR = 4.0, SIX = 6.0, & EIGHT = 8.0, TEN = 10.0, TOLMIN = 5.D-3) C C Initialization: C NSAVE = 200 JS = 5*NXINIT TENU = TEN*U BIG = ONE/U AVGRHO = BIG LPRINT = MIN(IPRINT,3) MAXNEV = 0 IF (REG(1)) THEN DX = SQRT(U)*MAX(ONE,ABS(A)) CALL COEFF(A+DX,PX,QX,RX) IF (FLAG .LT. 0) RETURN RPATA = RX*PX ALPHA1 = ONE/RX CALL COEFF(A+TWO*DX,PX,QX,RX) ALPHA1 = ALPHA1*(RPATA-PX*RX)/(DX*FOUR) RPATA = SQRT(RPATA) ALPHA2 = SQRT(RPATA) ALPHA1 = (A1+A2*ALPHA1)/ALPHA2 ALPHA2 = ALPHA2*A2 NCOEFF = NCOEFF+2 ELSE ALPHA1 = ZERO ALPHA2 = ONE RPATA = ONE ENDIF IF (REG(2)) THEN DX = SQRT(U)*MAX(ONE,ABS(B)) CALL COEFF(B-DX,PX,QX,RX) IF (FLAG .LT. 0) RETURN RPATB = RX*PX BETA1 = ONE/RX CALL COEFF(B-TWO*DX,PX,QX,RX) BETA1 = BETA1*(RPATB-PX*RX)/(DX*FOUR) RPATB = SQRT(RPATB) BETA2 = SQRT(RPATB) BETA1 = (B1-B2*BETA1)/BETA2 BETA2 = BETA2*B2 NCOEFF = NCOEFF+2 ELSE BETA1 = ZERO BETA2 = ONE ENDIF ALPHA = A1P*A2-A1*A2P IF (CSPEC(2)) THEN IF (ALPHA .EQ. ZERO) THEN SCALE = SQRT(A1**2+A2**2) ELSE SCALE = SQRT(ALPHA) ENDIF ENDIF IF (CSPEC(1)) SCALE = SQRT(B1**2+B2**2) TOLMAX = MAX(TOL(1),TOL(2)) TOL1 = MIN(TOL(1),TOLMIN) TOL2 = MIN(TOL(2),TOLMIN) ABSERR = TOL1 RELERR = TOL2 KLVL = 1 DELTA = HALF C C Begin the Main loop over LEVEL. C DO 120 LEVEL = 0,MAXLVL IF (IPRINT .GE. 2) THEN WRITE (*,10) LEVEL,ITER WRITE (21,10) LEVEL,ITER 10 FORMAT(' Level',I3,' of iteration',I3) ENDIF NRHO = 1 KRHO = 0 PSRHO = ZERO RHOSUM = ZERO ROLD = ZERO IASYMP = 0 EVSAVE = -BIG NEV = 0 JUMP = .FALSE. C C Compute integrals needed in asymptotic formulas. C QINT = ZERO RPINT = ZERO DO 20 I = 2,NXINIT XLEFT = STORE(I-1) H = (STORE(I)-XLEFT)/KLVL HALFH = HALF*H DO 15 J = 1,KLVL Z = XLEFT+HALFH XLEFT = XLEFT+H CALL PQRINT(Z,SQRTRP,QLNF) IF (FLAG .LT. 0) RETURN QINT = QINT+H*QLNF RPINT = RPINT+H*SQRTRP 15 CONTINUE 20 CONTINUE IF (QINT .GT. ONE/U) QINT = ZERO ZABS = MAX(MIN(ABSERR/100.0,RELERR)/10.0,TENU) ZREL = RELERR/10.0 DELTA = MAX(DELTA/SIX,TOL1+TOL2) C C Begin the secondary loop over NEV. C 25 IF (IASYMP .GE. 2) THEN AEV = ASYMEV(NEV,QINT,RPINT,ALPHA1,ALPHA2,BETA1,BETA2) EVHAT = AEV IF (IASYMP .LE. 3) GOTO 35 RSUBN = ASYMR(NEV,RPINT,RPATA,RPATB,SCALE) GOTO 45 ENDIF IF (HMIN/KLVL .LE. TENU) FLAG = -8 IF (FLAG .LT. 0) RETURN IF ((LEVEL .GT. 0) .AND. (NEV .LT. MIN(MAXNEV,NSAVE))) THEN EV = MAX(HALF*TOL1,DELTA,HALF*TOL2*ABS(EVOLD(NEV+1))) EVLOW = EVOLD(NEV+1)-EV EVHIGH = EVOLD(NEV+1)+EV ELSE IF (LEVEL .EQ. 0) THEN EV = ZERO ELSE EV = ASYMEV(NEV,QINT,RPINT,ALPHA1,ALPHA2,BETA1,BETA2) ENDIF ETA = MAX(HALF*TOL1,DELTA,HALF*TOL2*ABS(EV)) EVLOW = EV-ETA EVHIGH = EV+ETA ENDIF CALL BRCKET(NEV,EVLOW,EVHIGH,FLOW,FHIGH,ZABS,ZREL,STORE) IF ((LEVEL .EQ. 0) .AND. (NEV .EQ. 0)) DELTA = EVHIGH-EVLOW IF (FLAG .LT. 0) RETURN IF (ABS(EVHIGH-EVLOW) .GT. MAX(ZABS,ZREL*ABS(EVHIGH))) & CALL ZZERO(EVLOW,EVHIGH,FLOW,FHIGH,ZABS,ZREL,J,STORE) EVHAT = MIN(EVLOW,EVHIGH) CALL GETEF(EVHAT,EFNORM,LPRINT,STORE,EFIN) IF (FLAG .LT. 0) RETURN IF (CSPEC(2)) THEN IF (REG(1) .OR. (PNU(1) .EQ. ZERO) .OR. & (PNU(1) .EQ. ONE-EP(1))) THEN UX = ABS(STORE(NXINIT+1)) PDU = ABS(STORE(2*NXINIT+1)) IF (A2-A2P*EVHAT .NE. ZERO) THEN DENOM = SCALE*UX/ABS(A2-A2P*EVHAT) ELSE DENOM = SCALE*PDU/ABS(A1-A1P*EVHAT) ENDIF ELSE H = STORE(2)-STORE(1) UX = ABS(STORE(NXINIT+2)) PDU = ABS(STORE(2*NXINIT+2)) IF (UX .GE. PDU) THEN DENOM = UX/H**PNU(1) ELSE DENOM = PDU/(CP(1)*ABS(PNU(1))* & H**(EP(1)+PNU(1)-ONE)) ENDIF ENDIF ELSE IF (REG(2) .OR. (PNU(2) .EQ. ZERO) .OR. & (PNU(2) .EQ. ONE-EP(2))) THEN UX = ABS(STORE(2*NXINIT)) PDU = ABS(STORE(3*NXINIT)) IF (B2 .NE. ZERO) THEN DENOM = SCALE*UX/ABS(B2) ELSE DENOM = SCALE*PDU/ABS(B1) ENDIF ELSE H = STORE(NXINIT)-STORE(NXINIT-1) UX = ABS(STORE(2*NXINIT-1)) PDU = ABS(STORE(3*NXINIT-1)) IF (UX .GE. PDU) THEN DENOM = UX/H**PNU(2) ELSE DENOM = PDU/(CP(2)*ABS(PNU(2))* & H**(EP(2)+PNU(2)-ONE)) ENDIF ENDIF ENDIF IF (BIG*DENOM .GE. ONE) THEN EFNORM = ONE/DENOM**2 ELSE EFNORM = ONE/U**2 ENDIF C C Test for asymptotic EV. C AEV = ASYMEV(NEV,QINT,RPINT,ALPHA1,ALPHA2,BETA1,BETA2) IF (TEN*ABS(AEV-EVHAT) .GT. MAX(ABSERR,RELERR*ABS(EVHAT))) & THEN IASYMP = 0 ELSE IF (IASYMP .LT. 1) THEN IASYMP = 1 ELSE IF (IPRINT .GE. 2) THEN WRITE (*,30) NEV WRITE (21,30) NEV 30 FORMAT(' Switchover to asymptotic eigenvalues at', & ' Nev =',I8) ENDIF IASYMP = 2 ENDIF ENDIF IF (EFNORM .LT. BIG) THEN RSUBN = ONE/(ALPHA+EFNORM) ELSE RSUBN = ZERO ENDIF C C Test for asymptotic RsubN. C 35 ARN = ASYMR(NEV,RPINT,RPATA,RPATB,SCALE) IF (IASYMP .GE. 2) THEN C C Eigenvalues from asymptotic formulas; produce current RsubN. C CALL GETRN(AEV,ALPHA,CSPEC,SCALE,RSUBN,STORE) IF (FLAG .LT. 0) RETURN IF (ABS(ARN-RSUBN) .GT. MAX(ABSERR/100.0,RELERR*ARN)) & THEN IASYMP = 2 ELSE IF (IASYMP .LT. 3) THEN IASYMP = 3 ELSE IF (IPRINT .GE. 2) THEN WRITE (*,40) NEV WRITE (21,40) NEV 40 FORMAT(' Switchover to asymptotic RsubN at ', & 'Nev = ',I8) ENDIF IASYMP = 4 ENDIF ENDIF ENDIF 45 IF (NEV .LT. NSAVE) EVOLD(NEV+1) = EVHAT IF (NEV .GT. 0) ETA = HALF*(EVHAT+EVSAVE) IF (EVHAT .LT. CUTOFF) PSRHO = PSRHO+RSUBN 50 IF (T(NRHO) .LE. CUTOFF) THEN C C Use step functions for Rho(t). C IF (T(NRHO) .LE. EVHAT) THEN RHO(NRHO) = RHOSUM NRHO = NRHO+1 IF (NRHO .LE. NUMT) THEN GOTO 50 ELSE GOTO 85 ENDIF ENDIF ELSE IF (NEV .EQ. 0) GOTO 78 C C Use linear interpolation for Rho(t). C 55 IF (T(NRHO) .LE. ETA) THEN IF (JUMP) THEN 60 IF (T(NRHO) .LE. EVSAVE) THEN RHO(NRHO) = RHOSUM-ROLD NRHO = NRHO+1 IF (NRHO .LE. NUMT) THEN GOTO 60 ELSE JUMP = .FALSE. GOTO 85 ENDIF ELSE JUMP = .FALSE. 65 IF (T(NRHO) .LE. ETA) THEN RHO(NRHO) = RHOSUM NRHO = NRHO+1 IF (NRHO .GT. NUMT) GOTO 85 GOTO 65 ELSE GOTO 70 ENDIF ENDIF ENDIF IF ((NEV .LE. 1) .OR. (CUTOFF .GT. EVSAVE)) THEN UX = CUTOFF ELSE UX = ETAOLD ENDIF RHO(NRHO) = RHOSUM-ROLD*(ETA-T(NRHO))/(ETA-UX) NRHO = NRHO+1 IF (NRHO .LE. NUMT) THEN GOTO 55 ELSE GOTO 85 ENDIF ENDIF 70 IF ((RSUBN .GT. MAX(SIX*ROLD,TEN*TOLMAX,FOUR*AVGRHO)) & .AND. (EVHAT .GT. CUTOFF) .AND. (KRHO .GT. 4)) THEN C C Possible eigenvalue in the continuous spectrum. C LFLAG(3) = .TRUE. IF (IPRINT .GE. 1) THEN WRITE (*,75) EVHAT WRITE (21,75) EVHAT 75 FORMAT(' Large jump in the step spectral density', & ' function at',D17.10) WRITE (*,76) ITER,LEVEL,RSUBN WRITE (21,76) ITER,LEVEL,RSUBN 76 FORMAT(18X,'Iteration =',I2,', level = ',I2, & ', jump = ',D17.10) ENDIF JUMP = .TRUE. ENDIF ENDIF IF (NEV .GT. 0) ETAOLD = ETA 78 RHOSUM = RHOSUM+RSUBN ROLD = RSUBN EVSAVE = EVHAT C C Output requested information. C IF (IPRINT .GE. 3) THEN IF ((NEV .LE. 25) .OR. ((IASYMP .LE. 1) .AND. & (IPRINT .GE. 5))) THEN WRITE (21,80) NEV,EVHAT,RSUBN WRITE (*,80) NEV,EVHAT,RSUBN 80 FORMAT(' Nev =',I7,', EvHat =',D15.6,', RHat =',D15.6) ELSE IF ((NEV .LT. 100) .AND. (10*(NEV/10) .EQ. NEV)) THEN WRITE (21,80) NEV,EVHAT,RSUBN WRITE (*,80) NEV,EVHAT,RSUBN ELSE IF ((NEV .LT. 1000) .AND. (100*(NEV/100) .EQ. NEV)) & THEN WRITE (21,80) NEV,EVHAT,RSUBN WRITE (*,80) NEV,EVHAT,RSUBN ELSE IF (1000*(NEV/1000) .EQ. NEV) THEN WRITE (21,80) NEV,EVHAT,RSUBN WRITE (*,80) NEV,EVHAT,RSUBN ENDIF ENDIF ENDIF ENDIF ENDIF NEV = NEV+1 IF (EVHAT .GE. CUTOFF) THEN IF (KRHO .EQ. 5) THEN RSAVE(1) = RSAVE(2) RSAVE(2) = RSAVE(3) RSAVE(3) = RSAVE(4) RSAVE(4) = RSAVE(5) ELSE KRHO = KRHO+1 ENDIF RSAVE(KRHO) = RSUBN AVGRHO = ZERO DO 84 I = 1,KRHO AVGRHO = AVGRHO+RSAVE(I) 84 CONTINUE IF (KRHO .GT. 0) AVGRHO = AVGRHO/KRHO ENDIF GOTO 25 C C End of Nev loop ------------------ C 85 MAXNEV = NEV NUMEV = NUMEV+MAXNEV IF (IPRINT .GE. 3) THEN WRITE (*,90) MAXNEV WRITE (21,90) MAXNEV 90 FORMAT(' MaxNev = ',I8) ENDIF IF (IPRINT .GE. 4) THEN WRITE (*,95) WRITE (21,95) 95 FORMAT(9X,'t',18X,'RhoHat(t)') DO 105 J = 1,NUMT WRITE (*,100) T(J),RHO(J) WRITE (21,100) T(J),RHO(J) 100 FORMAT(F12.4,E31.15) 105 CONTINUE ENDIF C C Extrapolate interpolated approximations. C DONE = .TRUE. DO 110 J = 1,NUMT XTOL = MAX(TOL1,ABS(RHO(J))*TOL2) CALL EXTRAP(RHO(J),XTOL,LEVEL+1,NEXTRP,.TRUE.,.FALSE.,1, & STORE((MAXLVL+1)*J+JS),IPRINT,ERROR,RDONE) IF (RHO(J) .LT. ZERO) RHO(J) = ZERO IF (J .GT. 1) RHO(J) = MAX(RHO(J),RHO(J-1)) IF (ERROR .LE. HALF*XTOL) RDONE = .TRUE. DONE = RDONE .AND. DONE 110 CONTINUE IF (DONE) RETURN ABSERR = MAX(HALF*ABSERR,TENU) RELERR = MAX(HALF*RELERR,TENU) 120 CONTINUE FLAG = -3 RETURN END C----------------------------------------------------------------------- SUBROUTINE DSCRIP(LC,LPLC,TYPE,REG,CSPEC,CEV,CUTOFF,LASTEV, & A1,A1P,A2,A2P,B1,B2) INTEGER LASTEV LOGICAL LC(*),LPLC(*),TYPE(4,*),REG(*),CSPEC(*) DOUBLE PRECISION CEV(*),CUTOFF,A1,A1P,A2,A2P,B1,B2 C C Output (if requested) a description of the spectrum. C WRITE (21,*) IF (TYPE(3,1) .AND. TYPE(3,2)) THEN C C Category 1 C WRITE (21,123) 1 WRITE (21,100) WRITE (21,121) WRITE (21,102) WRITE (21,110) IF (REG(1)) THEN WRITE (21,112) ELSE WRITE (21,113) WRITE (21,114) IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF WRITE (21,119) A1,A1P,A2,A2P ENDIF WRITE (21,111) IF (REG(2)) THEN WRITE (21,112) ELSE WRITE (21,113) WRITE (21,114) IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF WRITE (21,119) B1,B2 ENDIF ENDIF C IF ((TYPE(3,1) .AND. CSPEC(2)) .OR. & (TYPE(3,2) .AND. CSPEC(1))) THEN C C Category 2 C WRITE (21,123) 2 WRITE (21,100) WRITE (21,103) CUTOFF IF (LASTEV .EQ. -5) THEN WRITE (21,105) WRITE (21,120) ELSE IF (LASTEV .EQ. 0) THEN WRITE (21,107) ELSE IF (LASTEV .EQ. 1) THEN WRITE (21,108) ELSE WRITE (21,109) LASTEV ENDIF ENDIF ENDIF WRITE (21,110) IF (REG(1)) THEN WRITE (21,112) ELSE WRITE (21,113) IF (CSPEC(1)) THEN WRITE (21,116) CEV(1) ELSE WRITE (21,114) WRITE (21,119) A1,A1P,A2,A2P ENDIF IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF ENDIF WRITE (21,111) IF (REG(2)) THEN WRITE (21,112) ELSE WRITE (21,113) IF (CSPEC(2)) THEN WRITE (21,116) CEV(2) ELSE WRITE (21,114) WRITE (21,119) B1,B2 ENDIF IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF ENDIF C IF ((TYPE(3,1) .AND. (TYPE(4,2) .AND. LC(2) .AND. & (.NOT. CSPEC(2)))) .OR. & (TYPE(3,2) .AND. (TYPE(4,1) .AND. LC(1) .AND. & (.NOT. CSPEC(1))))) THEN C C Category 3 C WRITE (21,123) 3 WRITE (21,100) WRITE (21,106) WRITE (21,102) WRITE (21,110) IF (REG(1)) THEN WRITE (21,112) ELSE WRITE (21,113) IF (TYPE(4,1)) THEN WRITE (21,115) ELSE WRITE (21,114) WRITE (21,119) A1,A1P,A2,A2P ENDIF IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF ENDIF WRITE (21,111) IF (REG(2)) THEN WRITE (21,112) ELSE WRITE (21,113) IF (TYPE(4,2)) THEN WRITE (21,115) ELSE WRITE (21,114) WRITE (21,119) B1,B2 ENDIF IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF ENDIF C IF ((TYPE(3,1) .AND. ((.NOT. LC(2)) .AND. TYPE(4,2) .AND. & (.NOT. CSPEC(2)))) .OR. & (TYPE(3,2) .AND. ((.NOT. LC(1)) .AND. TYPE(4,1) .AND. & (.NOT. CSPEC(1))))) THEN C C Category 4 C WRITE (21,123) 4 WRITE (21,100) WRITE (21,104) WRITE (21,110) IF (REG(1)) THEN WRITE (21,112) ELSE WRITE (21,113) IF (TYPE(4,1)) THEN WRITE (21,115) ELSE WRITE (21,114) WRITE (21,119) A1,A1P,A2,A2P ENDIF IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF ENDIF WRITE (21,111) IF (REG(2)) THEN WRITE (21,112) ELSE WRITE (21,113) IF (TYPE(4,2)) THEN WRITE (21,115) ELSE WRITE (21,114) WRITE (21,119) B1,B2 ENDIF IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF ENDIF C IF ((LC(1) .AND. TYPE(4,1)) .AND. (LC(2) .AND. TYPE(4,2))) THEN C C Category 5 C WRITE (21,123) 5 WRITE (21,100) WRITE (21,106) WRITE (21,102) WRITE (21,110) WRITE (21,113) WRITE (21,115) WRITE (21,117) WRITE (21,111) WRITE (21,113) WRITE (21,115) WRITE (21,117) ENDIF C IF ((TYPE(4,1) .AND. (.NOT. LC(1)) .AND. (.NOT. CSPEC(1))) .AND. & (TYPE(4,2) .AND. (.NOT. LC(2)) .AND. (.NOT. CSPEC(2)))) THEN C C Category 6 C WRITE (21,123) 6 WRITE (21,122) WRITE (21,110) WRITE (21,113) WRITE (21,115) WRITE (21,118) WRITE (21,111) WRITE (21,113) WRITE (21,115) WRITE (21,118) ENDIF C IF ((TYPE(4,1) .AND. TYPE(4,2) .AND. .NOT. (CSPEC(1) .OR. & CSPEC(2))) .AND. ((LC(1) .AND. (.NOT. LC(2))) .OR. & ((.NOT. LC(1)) .AND. LC(2)))) THEN C C Category 7 C WRITE (21,123) 7 WRITE (21,104) WRITE (21,110) WRITE (21,113) WRITE (21,115) IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF WRITE (21,111) WRITE (21,113) WRITE (21,115) IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF C IF ((LC(1) .AND. TYPE(4,1) .AND. CSPEC(2)) .OR. & (LC(2) .AND. TYPE(4,2) .AND. CSPEC(1))) THEN C C Category 8 C WRITE (21,123) 8 WRITE (21,100) WRITE (21,103) CUTOFF WRITE (21,110) WRITE (21,113) IF (CSPEC(1)) THEN WRITE (21,116) CEV(1) ELSE WRITE (21,115) ENDIF IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF WRITE (21,111) WRITE (21,113) IF (CSPEC(2)) THEN WRITE (21,116) CEV(2) ELSE WRITE (21,115) ENDIF IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF C IF ((((.NOT. LC(1)) .AND. TYPE(4,1) .AND. (.NOT. CSPEC(1))) & .AND. CSPEC(2)) .OR. & (((.NOT. LC(2)) .AND. TYPE(4,2) .AND. (.NOT. CSPEC(2))) & .AND. CSPEC(1))) THEN C C Category 9 C WRITE (21,123) 9 WRITE (21,101) WRITE (21,110) WRITE (21,113) IF (CSPEC(1)) THEN WRITE (21,116) CEV(1) ELSE WRITE (21,115) ENDIF IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF WRITE (21,111) WRITE (21,113) IF (CSPEC(2)) THEN WRITE (21,116) CEV(2) ELSE WRITE (21,115) ENDIF IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF C IF (CSPEC(1) .AND. CSPEC(2)) THEN C C Category 10 C WRITE (21,123) 10 WRITE (21,101) WRITE (21,103) CUTOFF IF (LASTEV .EQ. -5) THEN WRITE (21,120) ELSE IF (LASTEV .EQ. 0) THEN WRITE (21,107) ELSE IF (LASTEV .EQ. 1) THEN WRITE (21,108) ELSE WRITE (21,109) LASTEV ENDIF ENDIF ENDIF WRITE (21,110) WRITE (21,113) WRITE (21,116) CEV(1) IF (LC(1)) THEN IF (LPLC(1)) WRITE (21,117) ELSE IF (LPLC(1)) WRITE (21,118) ENDIF WRITE (21,111) WRITE (21,113) WRITE (21,116) CEV(2) IF (LC(2)) THEN IF (LPLC(2)) WRITE (21,117) ELSE IF (LPLC(2)) WRITE (21,118) ENDIF ENDIF WRITE (21,*) RETURN 100 FORMAT(' This problem has simple spectrum.') 101 FORMAT(' This problem may have non-simple spectrum.') 102 FORMAT(' There is no continuous spectrum.') 103 FORMAT(' There is continuous spectrum in [Ev, infinity) where' & ,' Ev =',G15.6) 104 FORMAT(' There is continuous spectrum consisting of the entire ', & 'real line.') 105 FORMAT(' The set of eigenvalues is bounded below.') 106 FORMAT(' There are infinitely many negative and infinitely many', & ' positive',/,' eigenvalues (unbounded in either', & ' direction).') 107 FORMAT(' There appear to be no eigenvalues below the start of', & ' the continuous spectrum.') 108 FORMAT(' There appears to be 1 eigenvalue below the start of the', & ' continuous spectrum.') 109 FORMAT(' There appear to be ',I12,' eigenvalues below the start' & ,/,' of the continuous spectrum.') 110 FORMAT(' At endpoint A') 111 FORMAT(' At endpoint B') 112 FORMAT(' the problem is regular;') 113 FORMAT(' the problem is singular;') 114 FORMAT(' it is nonoscillatory for all Ev.') 115 FORMAT(' it is oscillatory for all Ev.') 116 FORMAT(' it is nonoscillatory for Ev <',G15.6,' and oscillatory & otherwise.') 117 FORMAT(' It is limit circle.') 118 FORMAT(' It is limit point.') 119 FORMAT(' The constants for the Friedrichs boundary conditions', & ' are',/,4E18.8) 120 FORMAT(' There appear to be infinitely many eigenvalues below', & ' the start',/,' of the continuous spectrum.') 121 FORMAT(' There are infinitely many eigenvalues, bounded below.') 122 FORMAT(' The nature of the spectrum is unknown; there is likely', & ' to be ',/,' continuous spectrum.') 123 FORMAT(' The spectral category is',I3,'.') END C----------------------------------------------------------------------- SUBROUTINE EXTRAP(V,TOL,IROW,MAXCOL,FULL,TIGHT,MODE,VSAVE, & IPRINT,ERROR,DONE) INTEGER IROW,MAXCOL,MODE,IPRINT LOGICAL FULL,TIGHT,DONE DOUBLE PRECISION V,TOL,VSAVE(*),ERROR C C Use Richardson's h**2 extrapolation (based on doubling) when C suitable, otherwise use Wynn's acceleration scheme. C C Input: C V = real value at current level. C TOL = real tolerance. C IROW = integer giving current