***********************************************************************
* *
* 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 row index (1 .le. IROW).
C MAXCOL = integer giving maximum number of columns in table.
C FULL = logical, True iff entire table is to be computed.
C TIGHT = logical, True for conservative convergence tests.
C MODE = integer, value is
C 0 both Richardson and Wynn algorithms can be used;
C 1 only Richardson is used;
C 2 only Wynn is used.
C IPRINT = integer controlling amount of printing.
C Output:
C V = real, best output estimate.
C VSAVE = real vector, holds previous level values.
C DONE = logical, True iff Error is sufficiently small.
C
C If FULL is True, then the entire acceleration array is produced
C (through row IROW); if False, then only the next row is appended.
C Hence, the choice of FULL = True requires more work, but it may
C save some global storage. The vector VSAVE contains the V values
C for levels 0 through max(IROW-1,MAXCOL).
C
INTEGER I,IMIN,J,MAXJ,NCOL
DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER,
& DIFF,EPS,ETEMP,R(12,8),RATIO(12),RHIGH,RLOW,
& RTOL,T,TOL1,TOL2,VTEMP,W(40,11),
& ZERO,TENTH,HALF,ONE,TWO
COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER
SAVE R,W
C
C The local arrays RATIO(*), R(*,*), and W(*,*) must be declared to
C have at least as many rows as the value of MAXLVL initialized in
C routine START.
C
PARAMETER (ZERO = 0.0, TENTH = 0.1D0, HALF = 0.5D0, ONE = 1.0,
& TWO = 2.0)
C
EPS = 0.2
DONE = .FALSE.
ETEMP = ONE/U
VTEMP = V
VSAVE(IROW) = V
TOL1 = TOL
TOL2 = TOL/3.0
IF (MODE .EQ. 2) THEN
MAXJ = MAXCOL
GOTO 40
ELSE
MAXJ = 11
ENDIF
RLOW = MAX(3.0,4.2-0.2*(IROW-1))
RHIGH = MIN(5.0,3.8+0.2*(IROW-1))
RTOL = U
C
C Analyze the rate of convergence to determine NCOL and tolerances.
C
NCOL = 2
DO 10 I = 1,IROW
R(I,1) = VSAVE(I)
RATIO(I) = ONE/U
IF (I .GE. 3) THEN
T = R(I,1)-R(I-1,1)
IF (T .NE. ZERO) THEN
RATIO(I) = (R(I-1,1)-R(I-2,1))/T
ELSE
V = R(I,1)
DONE = .TRUE.
RETURN
ENDIF
IF (((RATIO(I) .GE. RLOW) .AND. (RATIO(I) .LE. RHIGH))
& .OR. (.NOT. TIGHT)) THEN
RTOL = TOL1
NCOL = MIN(MAXCOL,NCOL+1)
ELSE
RTOL = TOL2
IF (RATIO(I) .LT. ZERO) RTOL = TENTH*TOL2
IF (RATIO(I) .LT. TWO) NCOL = 2
ENDIF
ENDIF
10 CONTINUE
IF (FULL) THEN
IMIN = 2
ELSE
IMIN = IROW
ENDIF
C
C Use Richardson's h^2 extrapolation. The number of columns used
C is a function of the amount of data (IROW), the requested order
C (MAXCOL), the observed rate of convergence (NCOL), and the amount
C of storage allocated to R(*,*).
C
DO 30 I = IMIN,IROW
DO 20 J = 2,MIN(I,NCOL,8)
DIFF = (R(I,J-1)-R(I-1,J-1))/(4**(J-1)-1)
R(I,J) = R(I,J-1)+DIFF
IF ((.NOT. FULL) . OR. (I .EQ. IROW)) THEN
T = ABS(DIFF)
IF (T .LE. ETEMP) THEN
ETEMP = T
VTEMP = R(I,J)
IF (T .LE. RTOL) THEN
DONE = .TRUE.
V = VTEMP
ERROR = ETEMP
RETURN
ENDIF
ENDIF
ENDIF
20 CONTINUE
30 CONTINUE
V = VTEMP
ERROR = ETEMP
IF (IROW .LT. 4) RETURN
C
C Test for rate of convergence other than second order.
C
IF (ABS(RATIO(IROW)-RATIO(IROW-1))+ABS(RATIO(IROW-1)-
& RATIO(IROW-2)) .GT. EPS) RETURN
IF ((3.5 .LT. RATIO(IROW)) .AND. (RATIO(IROW) .LT. 4.5)) RETURN
IF (RATIO(IROW) .LT. ONE) RETURN
IF (MODE .NE. 0) RETURN
C
C Use Wynn's algorithm.
C
40 IF (IPRINT .GE. 4) THEN
DIFF = LOG(ABS(RATIO(IROW)))/LOG(TWO)
WRITE (21,50) DIFF
50 FORMAT(' In EXTRAP: using Wynn`s acceleration; rate = ',F8.5)
ENDIF
W(1,1) = VSAVE(1)
ETEMP = ONE/U
DO 60 I = 2,IROW
W(I,1) = VSAVE(I)
DIFF = W(I,1)-W(I-1,1)
IF ((I .EQ. IROW) .OR. (MODE .EQ. 2)) THEN
T = ABS(DIFF)
IF (T .LE. ETEMP) THEN
ETEMP = T
VTEMP = W(I,1)
IF (T .LE. TOL2) THEN
V = VTEMP
ERROR = ETEMP
DONE = .TRUE.
RETURN
ENDIF
ENDIF
ENDIF
IF (DIFF .NE. ZERO) THEN
W(I,2) = ONE/DIFF
ELSE
V = W(I,1)
ERROR = ZERO
DONE = .TRUE.
RETURN
ENDIF
60 CONTINUE
DO 80 J = 3,MIN(IROW,MAXJ)
DO 70 I = J,IROW
DIFF = W(I,J-1)-W(I-1,J-1)
IF (DIFF .NE. ZERO) THEN
DIFF = ONE/DIFF
ELSE
IF (MOD(J,2) .EQ. 0) THEN
V = W(I,J-1)
ERROR = ZERO
DONE = .TRUE.
RETURN
ELSE
DIFF = ONE/U**2
ENDIF
ENDIF
W(I,J) = W(I-1,J-2)+DIFF
IF ((MOD(J,2) .EQ. 1) .AND. ((I .EQ. IROW) .OR.
& (MODE .EQ. 2))) THEN
T = ABS(DIFF)
IF (T .LE. ETEMP) THEN
ETEMP = T
VTEMP = W(I,J)
IF (T .LE. TOL2) THEN
V = VTEMP
ERROR = ETEMP
DONE = .TRUE.
RETURN
ENDIF
ENDIF
ENDIF
70 CONTINUE
80 CONTINUE
V = VTEMP
ERROR = ETEMP
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE GETEF(EV,EFNORM,IPRINT,X,EFIN)
INTEGER IPRINT
DOUBLE PRECISION EV,EFNORM,X(*)
LOGICAL EFIN(2,2)
C
C Compute an eigenfunction for one fixed mesh.
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& KCLASS(2),I,J,JDU,JLAST,JS,JU,JX,KLVL,MIDDLE(2),MODE
LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2),
& ALLOK,SYMM
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),
& CHI,DPSI,DV,DW,FNORM,FSCALE,FSUM,H,HALFH,HOM,OM,OSCALE,
& PDV,PDW,PN,PROD,PSI,RATIO,RN,RNORM,RSCALE,RSUM,SCALE,
& T,TAU,TAUHH,TAUMAX,V,VNEW,W,WNEW,XLEFT,XRIGHT,Z,
& ZERO,C10M4,HALF,ONE,TWO,THREE,FIVE,C15,C21
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, C10M4 = 1.D-4, HALF = 0.5D0, ONE = 1.0,
& TWO = 2.0, THREE = 3.0, FIVE = 5.0, C15 = 15.0,
& C21 = 21.0)
C
KLVL = 2**LEVEL
C
C For this EV and mesh, calculate FSCALE, RSCALE, and MIDDLE(*).
C FSCALE = sum of logs of scale factors 1 through m
C RSCALE = sum of logs of scale factors m+1 through N
C MIDDLE(1), MIDDLE(2) describe the coordinates of the matching
C point M for the shooting; in particular,
C M = X(MIDDLE(1)-1) + MIDDLE(2)*H(MIDDLE(1)-1)/2**LEVEL
C The matching point is chosen to be roughly (a+b)/2 if either
C a = -b and Tau(x) > 0 near 0, or if Tau(x) > 0 for all x;
C otherwise, it is chosen to roughly maximize Tau(x).
C
FSCALE = ZERO
RSCALE = ZERO
TAUMAX = -ONE/U
ALLOK = .TRUE.
SYMM = .FALSE.
IF (A .EQ. -B) SYMM = .TRUE.
EFIN(1,1) = .TRUE.
EFIN(2,1) = .TRUE.
EFIN(1,2) = .TRUE.
EFIN(2,2) = .TRUE.
JX = (NXINIT+1)/2
DO 20 I = 2,NXINIT
MODE = 0
XLEFT = X(I-1)
H = X(I)-XLEFT
IF (I .EQ. 2) THEN
IF (KCLASS(1) .GE. 9) MODE = 1
IF (.NOT. AFIN) THEN
MODE = 3
XLEFT = ZERO
H = -ONE/X(2)
ENDIF
ENDIF
IF (I .EQ. NXINIT) THEN
IF (KCLASS(2) .GE. 9) MODE = 2
IF (.NOT. BFIN) THEN
MODE = 4
H = ONE/X(I-1)
XLEFT = -H
ENDIF
ENDIF
H = H/KLVL
HALFH = HALF*H
DO 10 J = 1,KLVL
Z = XLEFT+HALFH
XLEFT = XLEFT+H
CALL STEP(Z,H,EV,PN,RN,TAU,OM,HOM,PSI,DPSI,SCALE,MODE)
IF (TAU .LT. ZERO) ALLOK = .FALSE.
IF (TAU .GT. TAUMAX) THEN
TAUMAX = TAU
MIDDLE(1) = I
MIDDLE(2) = J
FSCALE = FSCALE+RSCALE+SCALE
OSCALE = SCALE
RSCALE = ZERO
ELSE
RSCALE = RSCALE+SCALE
ENDIF
10 CONTINUE
IF ((I .EQ. JX) .AND. (TAU .LT. ZERO)) SYMM = .FALSE.
20 CONTINUE
IF (ALLOK .OR. SYMM) THEN
MIDDLE(1) = JX
MIDDLE(2) = MAX(KLVL-1,1)
ENDIF
IF ((.NOT. AFIN) .OR. (.NOT. BFIN)) THEN
C
C Don't split near infinity!
C
IF ((.NOT. AFIN) .AND. (MIDDLE(1) .EQ. 2)) THEN
IF (REG(2)) THEN
MIDDLE(1) = NXINIT
ELSE
MIDDLE(1) = JX
ENDIF
MIDDLE(2) = MAX(KLVL-1,1)
ENDIF
IF ((.NOT. BFIN) .AND. (MIDDLE(1) .EQ. NXINIT)) THEN
IF (REG(1)) THEN
MIDDLE(1) = 2
ELSE
MIDDLE(1) = JX
ENDIF
MIDDLE(2) = 1
ENDIF
ENDIF
IF ((LEVEL .GT. 1) .AND. (MIDDLE(2) .EQ. KLVL)) THEN
MIDDLE(2) = KLVL-1
FSCALE = FSCALE-OSCALE
RSCALE = RSCALE+OSCALE
ENDIF
IF (IPRINT .GE. 4) WRITE (21,21) MIDDLE(1),MIDDLE(2)
21 FORMAT(' Coordinates of matching point =',2I6)
JU = NXINIT
JDU = 2*NXINIT
JS = 3*NXINIT
C
C Shoot from x=A to the middle.
C
V = A2-A2P*EV
PDV = A1-A1P*EV
FNORM = ZERO
IF (KCLASS(1) .LT. 9) THEN
SCALE = MAX(ABS(V),ABS(PDV))
V = V/SCALE
PDV = PDV/SCALE
MODE = 0
ELSE
V = ONE
PDV = D(4,1)
MODE = 1
ENDIF
IF (.NOT. AFIN) MODE = 3
FSUM = -FSCALE
IF ((IPRINT .GE. 5) .AND. (MODE .EQ. 0)) WRITE (21,35)
& X(1),V,PDV,FSUM
DO 40 I = 2,MIDDLE(1)
X(JU+I-1) = V
X(JDU+I-1) = PDV
X(JS+I-1) = FSUM
IF (MODE .EQ. 0) THEN
XLEFT = X(I-1)
H = X(I)-XLEFT
ELSE
XLEFT = ZERO
IF (MODE .EQ. 1) THEN
H = ((X(2)-X(1))/D(1,1))**ETA(1,1)
ELSE
H = -ONE/X(2)
ENDIF
ENDIF
H = H/KLVL
HALFH = HALF*H
JLAST = KLVL
IF (I .EQ. MIDDLE(1)) JLAST = MIDDLE(2)
DO 30 J = 1,JLAST
Z = XLEFT+HALFH
XLEFT = XLEFT+H
CALL STEP(Z,H,EV,PN,RN,TAU,OM,HOM,PSI,DPSI,SCALE,MODE)
DV = PDV/PN
IF (ABS(TAU)*H*H .GE. C10M4) THEN
IF (TWO*(FSUM+SCALE) .GT. -UNDER) THEN
FNORM = FNORM+RN*(PSI*DPSI*(V*V-DV*DV/TAU)/TWO+
& V*PSI*DV*PSI)*EXP(TWO*(FSUM+SCALE))
IF (TWO*FSUM .GT. -UNDER) FNORM = FNORM+RN*
& EXP(TWO*FSUM)*H*(V*V+DV*DV/TAU)/TWO
ENDIF
ELSE
IF (TWO*FSUM .GT. -UNDER) THEN
TAUHH = TAU*H*H
FNORM = FNORM+RN*H*EXP(TWO*FSUM)*(
& V*V*(ONE+TAUHH*(TAUHH/FIVE-ONE)/THREE)
& +H*V*DV*(ONE+TAUHH*(TWO*TAUHH/C15-ONE)/THREE)
& +(H*DV)**2*(ONE+TAUHH*(TWO*TAUHH/C21-ONE)/
& FIVE)/THREE )
ENDIF
ENDIF
FSUM = FSUM+SCALE
VNEW = DPSI*V+PSI*DV
PDV = -PN*TAU*PSI*V+DPSI*PDV
V = VNEW
30 CONTINUE
IF (MODE .EQ. 1) THEN
C
C Convert from V(t) to u(x).
C
IF (ETA(2,1) .LT. ZERO) THEN
X(JU+1) = ONE/U
EFIN(1,1) = .FALSE.
ELSE
IF (ETA(2,1) .EQ. ZERO) THEN
X(JU+1) = D(2,1)
ELSE
X(JU+1) = ZERO
ENDIF
ENDIF
Z = ETA(1,1)*ETA(2,1)+EP(1)-ONE
IF (Z .LT. ZERO) THEN
X(JDU+1) = SIGN(ONE/U,ETA(2,1))
EFIN(2,1) = .FALSE.
ELSE
IF (Z .GT. ZERO) THEN
X(JDU+1) = ZERO
ELSE
X(JDU+1) = D(2,1)*CP(1)*ETA(1,1)*ETA(2,1)/
& D(1,1)**(ETA(1,1)*ETA(2,1))
ENDIF
ENDIF
H = X(2)-A
PN = CP(1)*H**(EP(1)-ONE)
T = (H/D(1,1))**ETA(1,1)
CHI = D(2,1)*T**ETA(2,1)
V = V*CHI
PDV = PN*ETA(1,1)*(ETA(2,1)*V+CHI*PDV*
& T**(ONE-TWO*EMU(1)))
IF (IPRINT .GE. 5) WRITE (21,35) X(1),X(JU+1),X(JDU+1)
ENDIF
MODE = 0
IF (IPRINT .GE. 5) WRITE (21,35) XLEFT,V,PDV,FSUM
35 FORMAT(G16.6,3D15.6)
40 CONTINUE
C
C Shoot from x=B to the middle.
C
RNORM = ZERO
MODE = 0
IF (KCLASS(2) .LT. 9) THEN
SCALE = MAX(ABS(B1),ABS(B2))
W = -B2/SCALE
PDW = B1/SCALE
MODE = 0
ELSE
W = -ONE
PDW = -D(4,2)
MODE = 2
ENDIF
IF (.NOT. BFIN) MODE = 4
RSUM = -RSCALE
IF ((IPRINT .GE. 5) .AND. (MODE .EQ. 0)) WRITE (21,35)
& X(NXINIT),W,PDW,RSUM
DO 60 I = NXINIT,MIDDLE(1),-1
X(JU+I) = W
X(JDU+I) = PDW
X(JS+I) = RSUM
IF (MODE .EQ. 0) THEN
XRIGHT = X(I)
H = XRIGHT-X(I-1)
ELSE
XRIGHT = ZERO
IF (MODE .EQ. 2) THEN
H = ((X(NXINIT)-X(NXINIT-1))/D(1,2))**ETA(1,2)
ELSE
H = ONE/X(I-1)
ENDIF
ENDIF
H = H/KLVL
HALFH = HALF*H
JLAST = KLVL
IF (I .EQ. MIDDLE(1)) JLAST = JLAST-MIDDLE(2)
DO 50 J = 1,JLAST
Z = XRIGHT-HALFH
XRIGHT = XRIGHT-H
CALL STEP(Z,H,EV,PN,RN,TAU,OM,HOM,PSI,DPSI,SCALE,MODE)
DW = PDW/PN
IF (ABS(TAU)*H*H .GE. C10M4) THEN
IF (TWO*(RSUM+SCALE) .GT. -UNDER) THEN
RNORM = RNORM+RN*(PSI*DPSI*(W*W-DW*DW/TAU)/TWO-
& W*PSI*DW*PSI)*EXP(TWO*(RSUM+SCALE))
IF (TWO*RSUM .GT. -UNDER) RNORM = RNORM+RN*
& EXP(TWO*RSUM)*H*(W*W+DW*DW/TAU)/TWO
ENDIF
ELSE
IF (TWO*RSUM .GT. -UNDER) THEN
TAUHH = TAU*H*H
RNORM = RNORM+RN*H*EXP(TWO*RSUM)*(
& W*W*(ONE+TAUHH*(TAUHH/FIVE-ONE)/THREE)
& -H*W*DW*(ONE+TAUHH*(TWO*TAUHH/C15-ONE)/THREE)
& +(H*DW)**2*(ONE+TAUHH*(TWO*TAUHH/C21-ONE)/
& FIVE)/THREE )
ENDIF
ENDIF
RSUM = RSUM+SCALE
WNEW = DPSI*W-PSI*DW
PDW = PN*TAU*PSI*W+DPSI*PDW
W = WNEW
50 CONTINUE
IF (MODE .EQ. 2) THEN
C
C Convert from V(t) to u(x).
C
IF (ETA(2,2) .LT. ZERO) THEN
X(JU+NXINIT) = -ONE/U
EFIN(1,2) = .FALSE.
ELSE
IF (ETA(2,2) .EQ. ZERO) THEN
X(JU+NXINIT) = -D(2,2)
ELSE
X(JU+NXINIT) = ZERO
ENDIF
ENDIF
Z = ETA(1,2)*ETA(2,2)+EP(2)-ONE
IF (Z .LT. ZERO) THEN
X(JDU+NXINIT) = SIGN(ONE/U,ETA(2,2))
EFIN(2,2) = .FALSE.
ELSE
IF (Z .GT. ZERO) THEN
X(JDU+NXINIT) = ZERO
ELSE
X(JDU+NXINIT) = D(2,2)*CP(2)*ETA(1,2)*ETA(2,2)/
& D(1,2)**(ETA(1,2)*ETA(2,2))
ENDIF
ENDIF
IF (IPRINT .GE. 5) WRITE(21,35) X(NXINIT),X(JU+NXINIT),
& X(JDU+NXINIT)
H = X(NXINIT)-X(NXINIT-1)
PN = CP(2)*H**(EP(2)-ONE)
T = (H/D(1,2))**ETA(1,2)
CHI = D(2,2)*T**ETA(2,2)
W = CHI*W
PDW = PN*ETA(1,2)*(CHI*PDW*T**(ONE-TWO*EMU(2))-ETA(2,2)*W)
ENDIF
IF (MODE .EQ. 4) THEN
X(JU+NXINIT) = ZERO
X(JDU+NXINIT) = ONE
IF (IPRINT .GE. 5) WRITE(21,35) X(NXINIT),X(JU+NXINIT),
& X(JDU+NXINIT)
ENDIF
MODE = 0
IF ((JLAST .NE. 0) .AND. (IPRINT .GE. 5))
& WRITE (21,35) XRIGHT,W,PDW,RSUM
60 CONTINUE
IF (ABS(W) .GE. ABS(PDW)) THEN
RATIO = V/W
IF (IPRINT .GE. 5) WRITE (21,61) RATIO*PDW-PDV,RATIO
61 FORMAT(' DuHat jump, ratio =',2D24.15)
ELSE
RATIO = PDV/PDW
IF (V*W*RATIO .LT. ZERO) RATIO = -RATIO
IF (IPRINT .GE. 5) WRITE (21,62) RATIO*W-V,RATIO
62 FORMAT(' UHat jump, ratio =',2D24.15)
ENDIF
C
C Calculate weighted 2-norm and scale approximate eigenfunction.
C
FSCALE = EXP(-FSUM)
RSCALE = EXP(-RSUM)
EFNORM = SQRT(FNORM*FSCALE**2+RNORM*(RATIO*RSCALE)**2)
SCALE = LOG(EFNORM)
IF (IPRINT .GE. 5) WRITE (21,65) EFNORM
65 FORMAT(' EFnorm =',D24.15)
DO 70 I = 1,NXINIT
TAU = X(JS+I)-SCALE
IF (TAU .LE. -UNDER) THEN
X(JU+I) = ZERO
X(JDU+I) = ZERO
ELSE
PROD = EXP(TAU)
IF (I .GE. MIDDLE(1)) PROD = PROD*RATIO
X(JU+I) = X(JU+I)*PROD
X(JDU+I) = X(JDU+I)*PROD
ENDIF
70 CONTINUE
IF (IPRINT .GE. 4) THEN
WRITE (21,75)
75 FORMAT(10X,'x',15X,'Uhat(x)',13X,'PUhat`(x)')
DO 85 I = 1,NXINIT
WRITE (21,80) X(I),X(JU+I),X(JDU+I)
80 FORMAT(G16.6,2D20.8)
85 CONTINUE
ENDIF
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE GETRN(EV,ALPHA,CSPEC,DENOM,RSUBN,X)
DOUBLE PRECISION EV,ALPHA,DENOM,RSUBN,X(*)
LOGICAL CSPEC(*)
C
C Compute the RsubN value from the weighted eigenfunction 2-norm
C when standard shooting is stable and an accurate eigenvalue is
C available (from the asymptotic formulas).
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& KCLASS(2), I,J,KLVL
LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2)
DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,URN,UNDER,
& CP(2),CR(2),CUTOFF,D(4,2),EMU(2),EP(2),EQLNF(2),
& ER(2),ETA(2,2),PNU(2),
& DPHI,DPSI,DU,DUSAVE,FNORM,H,HALFH,HOM,HSAVE,OM,PDU,PHI,
& PN,PSI,RN,SCALE,TAU,TAUHH,U,UNEW,USAVE,XLEFT,Z,
& ZERO,C10M4,HALF,ONE,TWO,THREE,FIVE,C15,C21
COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,URN,UNDER
COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG
COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT
COMMON /SLCLSS/CP,CR,CUTOFF,D,EMU,EP,EQLNF,ER,ETA,PNU,KCLASS
PARAMETER (ZERO = 0.0, C10M4 = 1.D-4, HALF = 0.5D0, ONE = 1.D0,
& TWO = 2.0, THREE = 3.0, FIVE = 5.0, C15 = 15.0,
& C21 = 21.0)
C
U = A2-A2P*EV
PDU = A1-A1P*EV
FNORM = ZERO
KLVL = 2**LEVEL
C
C Shoot from x=A to x=B.
C
DO 20 I = 2,NXINIT
XLEFT = X(I-1)
H = (X(I)-XLEFT)/KLVL
HALFH = HALF*H
DO 10 J = 1,KLVL
Z = XLEFT+HALFH
XLEFT = XLEFT+H
CALL STEP(Z,H,EV,PN,RN,TAU,OM,HOM,PSI,DPSI,SCALE,0)
SCALE = EXP(SCALE)
PHI = PSI*SCALE
DPHI = DPSI*SCALE
DU = PDU/PN
IF (ABS(TAU)*H*H .GE. C10M4) THEN
FNORM = FNORM+RN*(PHI*DPHI*(U*U-DU*DU/TAU)/TWO
& +U*PHI*DU*PHI+H*(U*U+DU*DU/TAU)/TWO)
ELSE
TAUHH = TAU*H*H
FNORM = FNORM+RN*H*(U*U*(ONE+TAUHH*(TAUHH/FIVE-ONE)/
& THREE)
& +H*U*DU*(ONE+TAUHH*(TWO*TAUHH/C15-ONE)/THREE)
& +(H*DU)**2*(ONE+TAUHH*(TWO*TAUHH/C21-ONE)/
& FIVE)/THREE )
ENDIF
IF ((I .EQ. NXINIT) .AND. (J .EQ. KLVL) .AND. CSPEC(1)) THEN
HSAVE = H
USAVE = ABS(U)
DUSAVE = ABS(PDU)
ENDIF
UNEW = DPHI*U+PHI*DU
PDU = -PN*TAU*PHI*U+DPHI*PDU
U = UNEW
IF ((I .EQ. 2) .AND. (J .EQ. 1) .AND. CSPEC(2)) THEN
HSAVE = H
USAVE = ABS(U)
DUSAVE = ABS(PDU)
ENDIF
10 CONTINUE
20 CONTINUE
IF (CSPEC(2)) THEN
IF (REG(1) .OR. (PNU(1) .EQ. ZERO) .OR.
& (PNU(1) .EQ. ONE-EP(1))) THEN
PHI = DENOM
ELSE
IF (USAVE .GE. DUSAVE) THEN
PHI = USAVE/HSAVE**PNU(1)
ELSE
PHI = DUSAVE/(CP(1)*ABS(PNU(1))*HSAVE**(EP(1)+PNU(1)-ONE))
ENDIF
ENDIF
ELSE
IF (REG(2) .OR. (PNU(2) .EQ. ZERO) .OR.
& (PNU(2) .EQ. ONE-EP(2))) THEN
PHI = DENOM
ELSE
IF (USAVE .GE. DUSAVE) THEN
PHI = USAVE/HSAVE**PNU(2)
ELSE
PHI = DUSAVE/(CP(2)*ABS(PNU(2))*HSAVE**(EP(2)+PNU(2)-ONE))
ENDIF
ENDIF
ENDIF
RSUBN = ONE/(ALPHA+FNORM/PHI**2)
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE MESH(JOB,NEV,X,G,H,QLNF,Z,TOL,HMIN)
LOGICAL JOB
DOUBLE PRECISION X(*),G(*),H(*),QLNF(*),Z(*),TOL,HMIN
C
C If JOB = True then calculate the initial mesh; redistribute so
C that H(*) is approximately equidistributed. If JOB = False
C then use the mesh input by the user.
C
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& KCLASS(2),I,ITS,J,JTOL,K,MAXITS,N,NADD
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),
& DX,ENDA,ENDB,EPS,EQMAX,EQMIN,EV,GAMMA,P1,P2,P3,QMAX,QMIN,
& Q1,Q2,Q3,R1,R2,R3,WEIGHT,Y,Y1,Y2,Y3,
& ZERO,TENTH,HALF,ONE,TWO,FOUR
LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2),
& DONE
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, HALF = 0.5D0, ONE = 1.0,
& TWO = 2.0, FOUR = 4.0)
SAVE ENDA,ENDB
DATA EPS/0.0001/
C
IF (.NOT. JOB) THEN
HMIN = B-A
DO 5 I = 2,NXINIT
HMIN = MIN(HMIN,X(I)-X(I-1))
5 CONTINUE
RETURN
ENDIF
N = NXINIT-1
EV = NEV
C
C Find an appropriate initial mesh.
C
IF (AFIN) THEN
X(1) = A
IF (BFIN) THEN
X(NXINIT) = B
DX = (B-A)/N
DO 10 I = 2,N
X(I) = X(1)+(I-1)*DX
10 CONTINUE
ELSE
IF (NEV .LT. 0) THEN
ENDB = X(NXINIT)
ELSE
X(NXINIT) = (1+4*NEV)*ENDB
ENDIF
Y1 = X(1)/(ONE+ABS(X(1)))
DX = (X(NXINIT)/(ONE+X(NXINIT))-Y1)/N
DO 11 I = 2,N
Y = Y1+(I-1)*DX
X(I) = Y/(ONE-ABS(Y))
11 CONTINUE
ENDIF
ELSE
IF (NEV .LT. 0) THEN
ENDA = X(1)
ELSE
X(1) = (1+4*NEV)*ENDA
ENDIF
IF (BFIN) THEN
X(NXINIT) = B
Y1 = X(NXINIT)/(ONE+ABS(X(NXINIT)))
DX = (Y1-X(1)/(ONE-X(1)))/N
DO 12 I = 2,N
Y = Y1-(I-1)*DX
X(NXINIT+1-I) = Y/(ONE-ABS(Y))
12 CONTINUE
ELSE
Y1 = X(1)/(ONE-X(1))
IF (NEV .LT. 0) THEN
ENDB = X(NXINIT)
ELSE
X(NXINIT) = (1+4*NEV)*ENDB
ENDIF
Y2 = X(NXINIT)/(ONE+X(NXINIT))
DX = (Y2-Y1)/N
DO 13 I = 2,N
Y = Y1+(I-1)*DX
X(I) = Y/(ONE-ABS(Y))
13 CONTINUE
IF (ABS(X((NXINIT+1)/2)) .LT. TOL) X((NXINIT+1)/2) = ZERO
ENDIF
ENDIF
JTOL = -LOG10(TOL)+HALF
IF (REG(1) .AND. REG(2)) THEN
MAXITS = 6
ELSE
MAXITS = 3
ENDIF
C
C Calculate H(*) and G(*).
C
ITS = 1
IF (.NOT. AFIN) THEN
IF (X(2) .GE. ZERO) X(2) = MAX(HALF*X(1),-ONE)
ENDIF
IF (.NOT. BFIN) THEN
IF (X(N) .LE. ZERO) X(N) = MIN(HALF*X(NXINIT),ONE)
ENDIF
QMIN = 1.E31
QMAX = -QMIN
20 GAMMA = ZERO
C
C Equidistribute { [Qmax - Q]^2 * max[abs(p') , abs(q') , abs(r')] }.
C
EQMAX = ZERO
EQMIN = 1.E31
DO 25 J = 1,N
DX = X(J+1)-X(J)
Y2 = X(J)+HALF*DX
CALL COEFF(Y2,P2,Q2,R2)
IF ((P2 .EQ. ZERO) .OR. (R2 .EQ. ZERO) .OR. (P2*R2 .LT. ZERO))
& THEN
FLAG = -15
RETURN
ENDIF
Y1 = MAX(Y2-EPS,X(1)+TWO*U*ABS(X(1)))
CALL COEFF(Y1,P1,Q1,R1)
Y3 = MIN(Y2+EPS,X(NXINIT)-TWO*U*ABS(X(NXINIT)))
CALL COEFF(Y3,P3,Q3,R3)
IF (LNF) THEN
H(J) = ABS(Q3-Q1)/(Y3-Y1)
QLNF(J) = Q2/R2
ELSE
H(J) = MAX(ABS(P3-P1),ABS(Q3-Q1),ABS(R3-R1))/(Y3-Y1)
Y1 = SQRT(SQRT(R1*P1))
Y2 = SQRT(SQRT(R2*P2))
Y3 = SQRT(SQRT(R3*P3))
Y = SQRT(P2/R2)
QLNF(J) = Q2/R2 + Y*((Y3-Y1)*(SQRT(P3/R3)-SQRT(P1/R1))/FOUR
& +(Y3-TWO*Y2+Y1)*Y)/(Y2*EPS**2)
IF (ABS(QLNF(J)) .LE. EPS) QLNF(J) = ZERO
ENDIF
QMAX = MAX(QMAX,QLNF(J))
QMIN = MIN(QMIN,QLNF(J))
25 CONTINUE
Y = MAX(QMAX-QMIN,ONE)
EV = 100.0+MAX(ZERO,QMIN)
DO 30 J = 1,N
DX = X(J+1)-X(J)
IF (QLNF(J) .LE. EV) THEN
WEIGHT = 3.0*((QMAX-QLNF(J))/Y)**2+ONE
H(J) = MAX(H(J)*WEIGHT,U)
ELSE
Y2 = TWO*DX*SQRT(QLNF(J)-EV)
IF (Y2 .LE. UNDER) THEN
WEIGHT = EXP(-Y2)
H(J) = MAX(WEIGHT*H(J),U)
ELSE
H(J) = U
ENDIF
ENDIF
IF ((.NOT. AFIN) .AND. (X(J+1) .LT. ZERO)) THEN
H(J) = MAX(H(J)*EXP(X(J+1)),U)
IF ((J .EQ. 1) .AND. (H(1) .EQ. U)) X(1) = HALF*(X(1)+X(2))
ENDIF
IF ((.NOT. BFIN) .AND. (X(J) .GT. ZERO)) THEN
H(J) = MAX(H(J)*EXP(-X(J)),U)
IF ((J .EQ. N) .AND. (H(N) .EQ. U)) X(NXINIT) =
& HALF*(X(N)+X(NXINIT))
ENDIF
EQMIN = MIN(EQMIN,H(J)*DX)
EQMAX = MAX(EQMAX,H(J)*DX)
30 CONTINUE
NCOEFF = NCOEFF+3*N
IF (EQMAX-EQMIN .LE. MAX(TENTH*EQMAX,U/TENTH)) GOTO 75
C
C Use a roughly locally quasi-uniform mesh.
C
GAMMA = ZERO
DO 35 I = 1,N
GAMMA = GAMMA+H(I)
35 CONTINUE
GAMMA = ONE/GAMMA
DO 45 I = 1,N
Y = ZERO
DO 40 J = 1,N
Y = MAX( Y,H(J)/(ONE+GAMMA*ABS(X(I)+X(I+1)-X(J)-X(J+1))
& *H(J)) )
40 CONTINUE
Z(I) = Y
45 CONTINUE
DO 50 I = 1,N
H(I) = Z(I)
50 CONTINUE
G(1) = ZERO
DO 55 J = 1,N
G(J+1) = G(J)+H(J)*(X(J+1)-X(J))
55 CONTINUE
GAMMA = G(N+1)/N
C
C Redistribution algorithm:
C
Y = GAMMA
I = 1
DO 65 J = 1,N
60 IF (Y .LE. G(J+1)) THEN
I = I+1
Z(I) = X(J)+(Y-G(J))/H(J)
Y = Y+GAMMA
GOTO 60
ENDIF
65 CONTINUE
Z(1) = X(1)
Z(NXINIT) = X(NXINIT)
DONE = .TRUE.
DO 70 J = 2,N
IF (ABS(Z(J)-X(J)) .GT. TENTH*(Z(J+1)-Z(J-1))) DONE = .FALSE.
X(J) = Z(J)
70 CONTINUE
IF (.NOT. DONE) THEN
IF (ITS .LT. MAXITS) THEN
ITS = ITS+1
GOTO 20
ENDIF
ENDIF
75 IF (.NOT. AFIN) THEN
IF (X(2) .GE. ZERO) X(2) = -ONE
ENDIF
IF (.NOT. BFIN) THEN
IF (X(N) .LE. ZERO) X(N) = ONE
ENDIF
DO 95 K = 1,2
NADD = 0
IF (KCLASS(K) .GT. 0) THEN
C
C Add Nadd extra points near endpoint K.
C
IF (KCLASS(K) .EQ. 1) THEN
NADD = MIN(MAX((JTOL+2)/3,1),4)
Y = TENTH
ENDIF
IF (KCLASS(K) .EQ. 2) THEN
NADD = MIN(MAX(JTOL,3),4)
Y = MIN(MAX(TENTH**(JTOL/3),1.D-3),1.D-2)
ENDIF
IF (KCLASS(K) .EQ. 3) NADD = MIN(MAX(JTOL/2,2),4)
IF (KCLASS(K) .EQ. 4) THEN
NADD = MIN(MAX((5+JTOL)/3,2),5)
Y = TENTH**(NADD-1)
ENDIF
IF (KCLASS(K) .EQ. 5) THEN
NADD = (2**JTOL)**0.4
NADD = MIN(MAX(NADD,3),6)
Y = MIN(MAX((0.1**JTOL)**(0.333),0.001),0.1)
ENDIF
IF (KCLASS(K) .EQ. 6) THEN
NADD = MIN(MAX(JTOL,2),8)
Y = 0.005
ENDIF
IF ((KCLASS(K) .EQ. 7) .OR. (KCLASS(K) .EQ. 10)) THEN
NADD = MIN(MAX((2*JTOL+6)/3,3),8)
Y = MIN(MAX(SQRT(TENTH*TOL),1.D-6),1.D-2)
ENDIF
IF (KCLASS(K) .EQ. 8) THEN
NADD = (2**JTOL)**0.4
NADD = MIN(MAX(NADD,2),6)
Y = MIN(MAX((0.1**JTOL)**(0.333),0.001),0.05)
ENDIF
IF (KCLASS(K) .EQ. 9) THEN
IF(LFLAG(5)) THEN
NADD = MIN(MAX((JTOL+4)/3,2),5)
Y = TENTH**(NADD-1)
ELSE
NADD = MIN(MAX(2+JTOL*(JTOL-3)/40,2),4)
Y = 0.25
ENDIF
ENDIF
IF (K .EQ. 1) THEN
DO 80 I = NXINIT,2,-1
X(I+NADD) = X(I)
80 CONTINUE
NXINIT = NXINIT+NADD
IF (AFIN) DX = X(2)-A
DO 85 I = 1,NADD
IF (AFIN) THEN
X(I+1) = A+DX*Y**((NADD-I+ONE)/NADD)
ELSE
IF (KCLASS(1) .NE. 1) THEN
X(NADD+2-I) = X(NADD+3-I)-
& (X(NADD+4-I)-X(NADD+3-I))*2.4
ELSE
X(NADD+2-I) = X(NADD+3-I)-
& (X(NADD+4-I)-X(NADD+3-I))
ENDIF
ENDIF
85 CONTINUE
ELSE
IF (BFIN) DX = B-X(NXINIT-1)
N = NXINIT-1
NXINIT = NXINIT+NADD
X(NXINIT) = B
DO 90 I = 1,NADD
IF (BFIN) THEN
X(NXINIT-I) = B-DX*Y**((NADD-I+ONE)/NADD)
ELSE
IF (KCLASS(2) .NE. 1) THEN
X(N+I) = X(N+I-1)+(X(N+I-1)-X(N+I-2))*2.4
ELSE
X(N+I) = X(N+I-1)+(X(N+I-1)-X(N+I-2))
ENDIF
ENDIF
90 CONTINUE
ENDIF
ENDIF
95 CONTINUE
IF (.NOT. AFIN) X(1) = -ONE/U
IF (.NOT. BFIN) X(NXINIT) = ONE/U
HMIN = X(NXINIT)-X(1)
DO 100 I = 2,NXINIT
HMIN = MIN(HMIN,X(I)-X(I-1))
100 CONTINUE
RETURN
END
C-----------------------------------------------------------------------------
SUBROUTINE POWER(X,F,N,TOL,IPRINT,EF,CF,OSC,EXACT,Y,IFLAG)
DOUBLE PRECISION X(*),F(*),TOL,EF,CF,Y(*)
INTEGER N,IPRINT,IFLAG
LOGICAL OSC,EXACT
C
C Find the power function which "dominates" the tabled
C coefficient function. The output is Cf and Ef such that
C f(x) is asymptotic to Cf*x^Ef .
C The vectors X(*) and F(*) hold the N input points:
C F(I) = f(X(I)) I = 1,...,N.
C Set IFLAG = 0 for normal return; 1 for uncertainty in Ef;
C 2 if uncertain about Cf (oscillatory).
C
DOUBLE PRECISION ERROR,TOLAIT,TOLMIN, ZERO,HALF,ONE
INTEGER K,NY
PARAMETER (ZERO = 0.0, HALF = 0.5D0, ONE = 1.0)
DATA TOLMIN/1.D-6/
C
C Estimate the exponent.
C
OSC = .FALSE.
NY = N-1
ERROR = 1.E30
TOLAIT = MIN(TOLMIN,TOL)
DO 10 K = 1,NY
IF ((F(K) .NE. ZERO) .AND. (F(K+1) .NE. ZERO)) THEN
Y(K) = LOG(ABS(F(K+1)/F(K)))/LOG(ABS(X(K+1)/X(K)))
ELSE
Y(K) = ZERO
ENDIF
10 CONTINUE
EF = Y(NY)
IF (IPRINT .GE. 5) THEN
WRITE (21,*) ' From POWER; E_k and c_k sequences:'
WRITE (21,15) (Y(K),K=1,NY)
15 FORMAT(4D19.10)
WRITE (21,*)
ENDIF
CALL AITKEN(EF,TOLAIT,NY,Y,ERROR)
K = EF+SIGN(HALF,EF)
IF (ABS(K-EF) .LE. SQRT(TOL)) EF = K
IF (ABS(ERROR) .GT. TOL*MAX(ONE,ABS(EF))) THEN
C
C There is uncertainty in the exponent.
C
IFLAG = 1
ENDIF
IF (ABS(EF) .LE. TOL) EF = ZERO
C
C Estimate the coefficient.
C
DO 20 K = 1,N-1
Y(K) = F(K)/ABS(X(K))**EF
20 CONTINUE
CF = Y(N-1)
IF (IPRINT .GE. 5) WRITE (21,15) (Y(K),K=1,N-1)
CALL AITKEN(CF,TOLAIT,N-1,Y,ERROR)
IF ((EF .GT. 20.) .AND. (ABS(CF) .LE. TOL)) THEN
C
C Coefficient probably has exponential behavior.
C
CF = SIGN(ONE,Y(N-1))
ELSE
IF ((ABS(ERROR) .GT. TOL*MAX(ONE,ABS(CF))) .OR.
& ((ABS(F(N)-CF*X(N)**EF) .GT. 20.0*TOL*ABS(F(N))) .AND.
& (EF .NE. ZERO ))) THEN
C
C There is uncertainty in the coefficient; call such
C cases oscillatory.
C
IFLAG = 2
OSC = .TRUE.
ENDIF
ENDIF
IF (ABS(CF) .GT. 1.D7) THEN
EXACT = .FALSE.
RETURN
ENDIF
K = CF+HALF
IF ((ABS(K-CF) .LE. SQRT(TOL)) .AND. (K .NE. 0)) CF = K
EXACT = .TRUE.
DO 30 K = 1,N
IF (ABS(F(K)-CF*X(K)**EF) .GT. TOL*ABS(F(K))) EXACT = .FALSE.
30 CONTINUE
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE PQRINT(X,SQRTRP,QLNF)
DOUBLE PRECISION X,SQRTRP,QLNF
C
C Evaluate the integrands needed for the asymptotic formulas.
C (1) The Liouville normal form potential Qlnf:
C Qlnf(t) = q/r + f"(t)/f with f = (pr)**.25 .
C (2) The term in the change of independent variable:
C sqrt(r/p) .
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT
DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER,
& EPS,FL,FM,FR,XDOTL,XDOTR,PX,QX,RX,Z,
& ZERO,TWO,FOUR
LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2)
COMMON /SLINT/FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT
COMMON /SLREAL/A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER
COMMON /SLLOG/AFIN,BFIN,COUNTZ,LFLAG,LNF,LC,OSC,REG
PARAMETER (ZERO = 0.0, TWO = 2.0, FOUR = 4.0)
IF (LNF) THEN
CALL COEFF(X,PX,QX,RX)
IF ((PX .EQ. ZERO) .OR. (RX .EQ. ZERO) .OR. (PX*RX .LT. ZERO))
& THEN
FLAG = -15
RETURN
ENDIF
NCOEFF = NCOEFF+1
QLNF = QX/RX
SQRTRP = SQRT(RX/PX)
ELSE
EPS = MIN(1.D-4,MIN(ABS(B-X),ABS(X-A))/TWO)
Z = X-EPS
CALL COEFF(Z,PX,QX,RX)
IF ((PX .EQ. ZERO) .OR. (RX .EQ. ZERO) .OR. (PX*RX .LT. ZERO))
& THEN
FLAG = -15
RETURN
ENDIF
XDOTL = SQRT(PX/RX)
FL = SQRT(SQRT(PX*RX))
Z = X+EPS
CALL COEFF(Z,PX,QX,RX)
XDOTR = SQRT(PX/RX)
FR = SQRT(SQRT(PX*RX))
CALL COEFF(X,PX,QX,RX)
SQRTRP = SQRT(RX/PX)
FM = SQRT(SQRT(PX*RX))
NCOEFF = NCOEFF+3
QLNF = QX/RX + ((FR-FM)*(XDOTR-XDOTL)/FOUR+
& (FR-TWO*FM+FL)/SQRTRP)/(EPS*EPS*FM*SQRTRP)
IF (ABS(QLNF) .LE. EPS) QLNF = ZERO
ENDIF
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE REGULR(JOB,JOBMSH,TOL,NEV,EV,IPRINT,NEXTRP,XEF,EF,PDEF,
& HMIN,STORE)
INTEGER NEV,IPRINT,NEXTRP
LOGICAL JOB,JOBMSH
DOUBLE PRECISION TOL(*),EV,XEF(*),EF(*),PDEF(*),HMIN,STORE(*)
***********************************************************************
* *
* REGULR calculates Sturm-Liouville eigenvalue and (optionally) *
* eigenfunction estimates for the problem described initially. *
* *
***********************************************************************
C
C Input parameters:
C JOB = logical variable describing tasks to be carried out.
C JOB = .True. iff an eigenfunction is to be calculated.
C JOBMSH = logical variable, JOBMSH = .True. iff initial mesh
C is a function of the eigenvalue index.
C CONS(*) = real vector of 8 input constants: A1, A1', A2, A2',
C B1, B2, A, B.
C TOL(*) = real vector of 6 tolerances.
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 e-function
C derivatives,
C TOL(6) is the rel. error tolerance for e-function
C derivatives.
C Eigenfunction tolerances need not be set if JOB is
C False. All absolute error tolerances must be
C positive; all relative must be at least 100 times
C the unit roundoff.
C NEV = integer index for the eigenvalue sought; NEV .GE. 0 .
C EV = real initial guess for eigenvalue NEV; accuracy is
C not at all critical, but if a good estimate is
C available some time may be saved.
C IPRINT = integer controlling amount of internal printing done.
C
C Output parameters:
C EV = real computed approximation to NEVth eigenvalue.
C XEF(*) = real vector of points for eigenfunction output.
C EF(*) = real vector of eigenfunction values: EF(i) is the
C estimate of u(XEF(i)). If JOB is False then this
C vector is not referenced.
C PDEF(*) = real vector of eigenfunction derivative values:
C PDEF(i) is the estimate of (pu')(XEF(i)). If JOB is
C False then this vector is not referenced.
C
C Auxiliary storage:
C STORE(*) = real vector of auxiliary storage, must be dimensioned
C at least max[100,26N]. (N the number of mesh points)
C
C Storage allocation in auxiliary vector (currently Maxlvl = 10):
C STORE(*)
C 1 -> N vector of mesh points X(*),
C N+1 -> 2N best current eigenfunction values,
C 2N+1 -> 3N best current derivative values,
C 3N+1 -> 4N scale factors in GETEF,
C 4N+1 -> (6+2*Maxlvl)N intermediate eigenfunction values.
C-----------------------------------------------------------------------
C Local variables:
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& I,II,J,JDU,JU,KDU,KK,KU
DOUBLE PRECISION A1,A1P,A2,A2P,B1,B2,A,B,U,UNDER,
& ABSERR,ALPHA1,ALPHA2,BETA1,BETA2,DELTA,EFNORM,ERROR,
& EVEXT(20),EVHAT,EVHIGH,EVLOW,FHIGH,FLOW,H,PDUMAX,QINT,QLNF,
& RELERR,RPINT,SQRTRP,TOLEXT,TOLMIN,TOLPDU,TOLSUM,TOL1,TOL2,
& TOL3,TOL4,TOL5,TOL6,UMAX,Z, ASYMEV,
& ZERO,HALF,ONE,THREE,FIVE,TEN
LOGICAL AFIN,BFIN,COUNTZ,LFLAG(6),LNF,LC(2),OSC(2),REG(2),
& DONE,EFDONE,EFIN(2,2),EVDONE,EXFULL
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, THREE = 3.0,
& FIVE = 5.0, TEN = 10.0, TOLMIN = 1.D-3)
C
ALPHA1 = ZERO
ALPHA2 = ONE
BETA1 = ZERO
BETA2 = ONE
IF (NEV .LT. 0) THEN
FLAG = -37
RETURN
ENDIF
EVDONE = .FALSE.
TOL1 = MIN(TOL(1),TOLMIN)
TOL2 = MIN(TOL(2),TOLMIN)
IF (.NOT. (REG(1) .AND. REG(2))) THEN
TOL1 = TOL1/THREE
TOL2 = TOL2/THREE
ENDIF
IF (JOB) THEN
TOL3 = MIN(TOL(3),TOLMIN)
TOL4 = MIN(TOL(4),TOLMIN)
TOL5 = MIN(TOL(5),TOLMIN)
TOL6 = MIN(TOL(6),TOLMIN)
ABSERR = TOL1/FIVE
RELERR = TOL2/FIVE
EXFULL = .TRUE.
ELSE
EFDONE = .TRUE.
ABSERR = TOL1/TEN
RELERR = TOL2/TEN
EXFULL = .FALSE.
ENDIF
TOLSUM = TOL1+TOL2
IF (JOBMSH) THEN
IF (NEV .GE. 0) THEN
II = NXINIT+16
CALL MESH(.TRUE.,NEV,STORE,STORE(II),STORE(2*II+1),
& STORE(3*II+1),STORE(4*II+1),TOLSUM,HMIN)
ENDIF
IF (IPRINT .GE. 1) THEN
WRITE (*,15) (STORE(I),I=1,NXINIT)
WRITE (21,15) (STORE(I),I=1,NXINIT)
15 FORMAT(' Level 0 mesh:',/,(5G15.6))
ENDIF
ENDIF
C
C Compute estimates for integrals in asymptotic formulas (accuracy
C is not all that critical).
C
QINT = ZERO
RPINT = ZERO
DO 20 I = 2,NXINIT
H = STORE(I)-STORE(I-1)
Z = STORE(I-1)+HALF*H
IF ((.NOT. AFIN) .AND. (I .EQ. 2)) THEN
H = STORE(3)-STORE(2)
Z = STORE(2)
ENDIF
IF ((.NOT. BFIN) .AND. (I .EQ. NXINIT)) THEN
H = STORE(NXINIT-1)-STORE(NXINIT-2)
Z = STORE(NXINIT-1)
ENDIF
CALL PQRINT(Z,SQRTRP,QLNF)
IF (FLAG .LT. 0) RETURN
QINT = QINT+H*QLNF
RPINT = RPINT+H*SQRTRP
20 CONTINUE
IF (QINT .GT. ONE/U) QINT = ZERO
IF (RPINT .GT. ONE/U) RPINT = ZERO
C
C Loop over the levels.
C
DO 60 LEVEL = 0,MAXLVL
IF (HMIN/2**LEVEL .LE. TEN*U) THEN
FLAG = -8
GOTO 70
ENDIF
C
C Find a bracket for the Nevth eigenvalue.
C
IF (LEVEL .EQ. 0) THEN
EV = ASYMEV(NEV,QINT,RPINT,ALPHA1,ALPHA2,BETA1,BETA2)
EV = MAX(EV,ZERO)
DELTA = HALF
ELSE
DELTA = MAX(TOLSUM*ABS(EVHAT),HALF*HALF*DELTA)
IF (LEVEL .GT. 1) THEN
ERROR = (EVEXT(LEVEL)-EVEXT(LEVEL-1))/THREE
IF (ABS(ERROR) .LE. 100.) THEN
EV = EVHAT+ERROR
ELSE
DELTA = ONE
EV = EVHAT
ENDIF
ELSE
EV = EVHAT
ENDIF
ENDIF
EVLOW = EV-DELTA
EVHIGH = EV+DELTA
IF (IPRINT .GE. 4) WRITE (21,25) EVLOW,EVHIGH
25 FORMAT(' In bracket:',2D24.15)
CALL BRCKET(NEV,EVLOW,EVHIGH,FLOW,FHIGH,ABSERR,RELERR,STORE)
IF (IPRINT .GE. 4) WRITE (21,30) EVLOW,EVHIGH
30 FORMAT(' Out bracket:',2D24.15)
DELTA = HALF*(EVHIGH-EVLOW)
IF (FLAG .LT. 0) RETURN
IF (ABS(EVHIGH-EVLOW) .GT. MAX(ABSERR,RELERR*ABS(EVHIGH)))
& THEN
CALL ZZERO(EVLOW,EVHIGH,FLOW,FHIGH,ABSERR,RELERR,J,STORE)
IF (J .NE. 0) THEN
FLAG = -7
RETURN
ENDIF
ENDIF
EVHAT = MIN(EVLOW,EVHIGH)
IF (IPRINT .GE. 1) THEN
WRITE (*,40) LEVEL,EVHAT
WRITE (21,40) LEVEL,EVHAT
40 FORMAT(' Level ',I3,' ; EvHat = ',D24.15)
ENDIF
EV = EVHAT
TOLEXT = MAX(TOL1,ABS(EV)*TOL2)
CALL EXTRAP(EV,TOLEXT,LEVEL+1,NEXTRP,EXFULL,.TRUE.,0,EVEXT,
& IPRINT,ERROR,EVDONE)
IF (JOB) THEN
CALL GETEF(EVHAT,EFNORM,IPRINT,STORE,efin)
IF (LEVEL .EQ. 0) THEN
UMAX = ONE
PDUMAX = ONE
C
C Set pointers to STORE(*).
C
JU = NXINIT
JDU = 2*NXINIT
KU = 4*NXINIT
KDU = (MAXLVL+5)*NXINIT
KK = MAXLVL+1
ENDIF
C
C Extrapolate eigenfunction values.
C
TOLEXT = MAX(TOL3,UMAX*TOL4)
TOLPDU = MAX(TOL5,PDUMAX*TOL6)
EFDONE = .TRUE.
IF (AFIN) THEN
UMAX = ABS(STORE(JU+1))
PDUMAX = ABS(STORE(JDU+1))
ELSE
UMAX = ZERO
PDUMAX = ZERO
ENDIF
IF (EFIN(1,1)) THEN
CALL EXTRAP(STORE(JU+1),TOLEXT,LEVEL+1,NEXTRP,EXFULL,
& .TRUE.,0,STORE(KU+1),0,ERROR,DONE)
EFDONE = EFDONE .AND. DONE
ENDIF
IF (EFIN(2,1)) THEN
CALL EXTRAP(STORE(JDU+1),TOLPDU,LEVEL+1,NEXTRP,EXFULL,
& .TRUE.,0,STORE(KDU+1),0,ERROR,DONE)
EFDONE = EFDONE .AND. DONE
ENDIF
DO 50 I = 2,NXINIT-1
II = KK*(I-1)+1
CALL EXTRAP(STORE(JU+I),TOLEXT,LEVEL+1,NEXTRP,EXFULL,
& .TRUE.,0,STORE(KU+II),0,ERROR,DONE)
EFDONE = EFDONE .AND. DONE
CALL EXTRAP(STORE(JDU+I),TOLPDU,LEVEL+1,NEXTRP,EXFULL,
& .TRUE.,0,STORE(KDU+II),0,ERROR,DONE)
EFDONE = EFDONE .AND. DONE
UMAX = MAX(UMAX,ABS(STORE(JU+I)))
PDUMAX = MAX(PDUMAX,ABS(STORE(JDU+I)))
50 CONTINUE
II = KK*(NXINIT-1)+1
IF (EFIN(1,2)) THEN
CALL EXTRAP(STORE(JU+NXINIT),TOLEXT,LEVEL+1,NEXTRP,
& EXFULL,.TRUE.,0,STORE(KU+II),0,ERROR,DONE)
EFDONE = EFDONE .AND. DONE
ENDIF
IF (EFIN(2,2)) THEN
CALL EXTRAP(STORE(JDU+NXINIT),TOLPDU,LEVEL+1,NEXTRP,
& EXFULL,.TRUE.,0,STORE(KDU+II),0,ERROR,DONE)
EFDONE = EFDONE .AND. DONE
ENDIF
IF (BFIN) THEN
UMAX = MAX(UMAX,ABS(STORE(JU+NXINIT)))
PDUMAX = MAX(PDUMAX,ABS(STORE(JDU+NXINIT)))
ENDIF
ABSERR = MAX(HALF*ABSERR,TEN*U)
RELERR = MAX(HALF*RELERR,TEN*U)
ENDIF
IF (EVDONE .AND. (LEVEL .GE. 2) .AND. EFDONE) GOTO 70
60 CONTINUE
IF (.NOT. EVDONE) THEN
FLAG = -1
RETURN
ENDIF
C
C Unload eigenfunction values.
C
70 IF (JOB) THEN
DO 80 I = 1,NXINIT
XEF(I) = STORE(I)
EF(I) = STORE(JU+I)
PDEF(I) = STORE(JDU+I)
80 CONTINUE
IF ((FLAG .GE. 0) .AND. (.NOT. EFDONE)) FLAG = -2
ENDIF
RETURN
END
C---------------------------------------------------------------------
SUBROUTINE SHOOT(EV,X,MU,FEV)
INTEGER MU
DOUBLE PRECISION EV,X(*),FEV
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& KCLASS(2),I,IMAX,J,K1,K2,KLVL,MODE,NSAVE,NZERO
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,
& CP(2),CR(2),CUTOFF,D(4,2),EMU(2),EP(2),EQLNF(2),
& ER(2),ETA(2,2),PNU(2),
& CHI,DPSI,DV,H,HALFH,HOMEGA,OMEGA,PDV,PHASE,PN,PSI,RN,
& SA1,SA2,SB1,SB2,SCALE,SGN,T,TAU,V,VNEW,XLEFT,X2,Z,
& ZERO,HALF,ONE,TWO,PI
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,
& PI = 3.141592653589793D0)
DATA IMAX/1000000/
C
C Make one shot across (A,B) for the current mesh using the scaled
C variable v(x). Count zeros if COUNTZ is True.
C
C Shoot from x=A to x=B.
C
V = A2-A2P*EV
PDV = A1-A1P*EV
NZERO = 0
NSAVE = NSGNF
SA1 = A1
SA2 = A2
SB1 = B1
SB2 = B2
KLVL = 2**LEVEL
MODE = 0
C
C Modify base sign count if necessary.
C
IF (((KCLASS(1) .GE. 9) .OR. (KCLASS(2) .GE. 9)) .AND.
& (EV .LT. CUTOFF)) THEN
IF (KCLASS(1) .GE. 9) THEN
SA1 = D(4,1)
SA2 = ONE
V = SA2
PDV = SA1
MODE = 1
ENDIF
IF (KCLASS(2) .GE. 9) THEN
SB1 = D(4,2)
SB2 = ONE
ENDIF
SGN = A2*B2
IF (SGN .EQ. ZERO) THEN
SGN = SA1*SB2+SA2*SB1
IF (SGN .EQ. ZERO) SGN = A1*B1
ENDIF
NSGNF = SIGN(ONE,SGN)
ENDIF
IF (.NOT. AFIN) MODE = 3
SCALE = MAX(ABS(V),ABS(PDV))
V = V/SCALE
PDV = PDV/SCALE
DO 20 I = 2,NXINIT
XLEFT = X(I-1)
H = X(I)-XLEFT
IF (MODE .EQ. 1) THEN
H = (H/D(1,1))**ETA(1,1)
XLEFT = ZERO
ENDIF
IF (MODE .EQ. 3) THEN
XLEFT = ZERO
H = -ONE/X(2)
ENDIF
IF ((KCLASS(2) .GE. 9) .AND. (I .EQ. NXINIT) .AND.
& (EV .LT. CUTOFF)) THEN
C
C Convert from u(x) to V(t) near x=b.
C
MODE = 2
T = (H/D(1,2))**ETA(1,2)
CHI = D(2,2)*T**ETA(2,2)
V = V/CHI
PN = CP(2)*H**(EP(2)-ONE)
PDV = (PDV/(PN*CHI*ETA(1,2))+ETA(2,2)*V)*T**(TWO*EMU(2)-ONE)
H = T
XLEFT = -H
ENDIF
IF ((.NOT. BFIN) .AND. (I .EQ. NXINIT)) THEN
MODE = 4
H = ONE/X(I-1)
XLEFT = -H
ENDIF
H = H/KLVL
HALFH = HALF*H
DO 10 J = 1,KLVL
Z = XLEFT+HALFH
CALL STEP(Z,H,EV,PN,RN,TAU,OMEGA,HOMEGA,PSI,DPSI,SCALE,MODE)
IF (FLAG .LT. 0) THEN
FEV = ZERO
RETURN
ENDIF
DV = PDV/PN
VNEW = DPSI*V+PSI*DV
XLEFT = XLEFT+H
IF (COUNTZ) THEN
C
C Count zeros of v(x).
C
IF (TAU .LE. ZERO) THEN
IF (VNEW*V .LT. ZERO) NZERO = NZERO+1
ELSE
IF (DV .EQ. ZERO) THEN
NZERO = NZERO+INT(HALF+HOMEGA/PI)
ELSE
PHASE = ATAN(V*OMEGA/DV)
K1 = PHASE/PI
X2 = (PHASE+HOMEGA)/PI
IF (X2 .LT. IMAX) THEN
K2 = X2
NZERO = NZERO+K2-K1
ELSE
NZERO = IMAX
ENDIF
IF (PHASE*(PHASE+HOMEGA) .LT. ZERO) NZERO = NZERO+1
ENDIF
NZERO = MIN(IMAX,NZERO)
ENDIF
ENDIF
PDV = -PN*TAU*PSI*V+DPSI*PDV
V = VNEW
10 CONTINUE
IF (MODE .EQ. 1) THEN
C
C Convert from V(t) back to u(x) near x=a.
C
PN = CP(1)*(X(2)-A)**(EP(1)-ONE)
T = ((X(2)-A)/D(1,1))**ETA(1,1)
CHI = D(2,1)*T**ETA(2,1)
PDV = PN*ETA(1,1)*CHI*(ETA(2,1)*V+PDV*T**(ONE-TWO*EMU(1)))
V = CHI*V
ENDIF
MODE = 0
20 CONTINUE
FEV = SB1*V+SB2*PDV
IF (COUNTZ) THEN
C
C Adjust zero count.
C
MU = NZERO
IF (A2P .NE. ZERO) THEN
IF (EV .GE. SA2/A2P) MU = NZERO+1
ENDIF
IF (SB2 .NE. ZERO) THEN
SGN = SIGN(ONE,NSGNF*FEV)
IF (MOD(MU,2) .EQ. 1) SGN = -SGN
IF (SGN .LT. ZERO) MU = MU+1
ENDIF
ENDIF
NSGNF = NSAVE
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE START(JOB,CONS,TOL,NEV,INDXEV,N,XEF,NUMT,T,NEXTRP,X)
INTEGER NEV,INDXEV(*),N,NUMT,NEXTRP
LOGICAL JOB(*)
DOUBLE PRECISION CONS(*),TOL(*),XEF(*),T(*),X(*)
C
C This routine tests the input data, initializes the labeled
C common blocks, and generates the first mesh. Check
C eigenfunction tolerances iff JOB(1) is True ,
C XEF(*) iff JOB(2) is True,
C NUMT, T(*) iff JOB(3) is True ,
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& I,K
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,
& UFLOW,URN,
& ZERO,HALF,ONE,PIOVR2,TWO,FIVE,HUNDRD
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,
& PIOVR2 = 1.5707963267948966D0, TWO = 2.0, FIVE = 5.0,
& HUNDRD = 100.0)
************************************************************************
* In the following DATA statement, initialize URN to an estimate for *
* the unit roundoff and UFLOW to a value somewhat less than *
* -ln(underflow). E.g., *
* for IEEE double precision use URN = 2.D-16, UFLOW = 650; *
* for VAXen double precision use URN = 1.D-17, UFLOW = 85; *
* for Crays (single precision) use URN = 7.D-15, UFLOW = 5000. *
* Exact values are not at all critical. *
* Here, we assume IEEE double precision: *
* *
DATA URN/2.D-16/,UFLOW/650.0/
C***********************************************************************
U = URN
UNDER = UFLOW
C
C Initialize.
C
A1 = CONS(1)
A1P = CONS(2)
A2 = CONS(3)
A2P = CONS(4)
A = CONS(7)
B1 = CONS(5)
B2 = CONS(6)
B = CONS(8)
NCOEFF = 0
FLAG = 0
C
C Test input.
C
IF (AFIN .AND. BFIN .AND. (A .GE. B)) FLAG = -34
IF (TOL(1) .LE. ZERO) FLAG = -35
IF (TOL(2) .LT. HUNDRD*U) FLAG = -36
IF (JOB(1)) THEN
IF (TOL(3) .LE. ZERO) FLAG = -35
IF (TOL(4) .LT. HUNDRD*U) FLAG = -36
IF (TOL(5) .LE. ZERO) FLAG = -35
IF (TOL(6) .LT. HUNDRD*U) FLAG = -36
ENDIF
IF (JOB(2)) THEN
IF (N .EQ. 1) THEN
FLAG = -30
ELSE
IF (AFIN .AND. (XEF(2) .LE. A)) FLAG = -39
DO 10 I = 3,N-1
IF (XEF(I-1) .GE. XEF(I)) FLAG = -39
10 CONTINUE
IF (BFIN .AND. (XEF(N-1) .GT. B)) FLAG = -39
ENDIF
IF ((.NOT. AFIN) .AND. (XEF(2) .GE. ZERO)) FLAG = -39
IF ((.NOT. BFIN) .AND. (XEF(N-1) .LE. ZERO)) FLAG = -39
ENDIF
IF ((JOB(2) .OR. JOB(3)) .AND. (NEV .GT. 0)) THEN
DO 20 I = 1,NEV
IF (INDXEV(I) .LT. 0) FLAG = -37
20 CONTINUE
ENDIF
IF (JOB(3)) THEN
IF (NUMT .LE. 0) FLAG = -38
DO 30 I = 2,NUMT
IF (T(I) .LE. T(I-1)) FLAG = -39
30 CONTINUE
ENDIF
IF (FLAG .LT. 0) RETURN
C
C Set MAXEXT, the maximum number of extrapolations allowed, and
C MAXINT, the maximum number of intervals in X(*) allowed when
C mesh is chosen by START.
C
C IMPORTANT NOTE: the size of various fixed arrays in this package
C depends on the value of MAXEXT in this FORTRAN77 implementation.
C If MAXEXT is increased, then more storage may have to be allocated
C to the columns of R(*,*) in EXTRAP.
C
MAXEXT = 6
MAXINT = 31
C
C Calculate maximum number of columns in extrapolation table
C and the maximum number of levels allowed.
C
K = -LOG10(TOL(2))
I = MAX(K+3,0)
NEXTRP = MIN(MAX(3,I/2),MAXEXT)
C
C Calculate the initial mesh.
C
IF (N .GT. 0) THEN
NXINIT = N
ELSE
NXINIT = MIN(2*NEXTRP+3,MAXINT)
IF (JOB(1)) N = NXINIT
ENDIF
IF (JOB(2)) THEN
A = XEF(1)
B = XEF(NXINIT)
DO 40 I = 1,NXINIT
X(I) = XEF(I)
40 CONTINUE
ENDIF
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE STEP(X,H,EV,PX,RX,TAU,OMEGA,HOMEGA,PSI,DPSI,SCLOG,MODE)
INTEGER MODE
DOUBLE PRECISION X,H,EV,PX,RX,TAU,OMEGA,HOMEGA,PSI,DPSI,SCLOG
C
C Evaluate the coefficient functions, the scaled basis function PSI,
C its derivative DPSI, and the log of the scale factor SCLOG.
C
INTEGER FLAG,LEVEL,MAXEXT,MAXINT,MAXLVL,NCOEFF,NSGNF,NXINIT,
& KCLASS(2)
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,
& CP(2),CR(2),CUTOFF,D(4,2),EMU(2),EP(2),EQLNF(2),
& ER(2),ETA(2,2),PNU(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
C
DOUBLE PRECISION DX,FP,FR,OVER,QX,T,TMU,Z,
& ZERO,HNDRTH,QUART,HALF,ONE,TWO,SIX,TWELVE,TWENTY
PARAMETER (ZERO = 0.0, HNDRTH = .01, QUART = 0.25D0,
& HALF = 0.5D0, ONE = 1.0, TWO = 2.0, SIX = 6.0,
& TWELVE = 12.0, TWENTY = 20.0)
DATA OVER/1.D8/
C
C Evaluate the coefficient functions at X and calculate TAU. The
C error flag FLAG is zero for a successful calculation; if p(x)
C or r(x) are zero, then FLAG is set to -15. If the argument for
C a trig function exceeds OVER then FLAG is set to -10.
C Proceed normally when Mode = 0; when Mode = 1 or 2 use the
C change of variable for "hard" problems; when Mode = 3 or 4 use
C the change of variable t = -1/x near infinity.
C
IF (MODE .EQ. 0) THEN
CALL COEFF(X,PX,QX,RX)
ELSE
IF (MODE .GE. 3) THEN
T = X
X = -ONE/T
CALL COEFF(X,PX,QX,RX)
T = T*T
PX = T*PX
QX = QX/T
RX = RX/T
ELSE
T = X
IF (ETA(1,MODE) .EQ. ONE) THEN
DX = D(1,MODE)*ABS(T)
ELSE
DX = D(1,MODE)*(ABS(T))**(ONE/ETA(1,MODE))
ENDIF
IF (MODE .EQ. 1) THEN
Z = A+DX
ELSE
Z = B-DX
ENDIF
CALL COEFF(Z,PX,QX,RX)
ENDIF
ENDIF
NCOEFF = NCOEFF+1
IF ((PX .EQ. ZERO) .OR. (RX .EQ. ZERO)) THEN
FLAG = -15
RETURN
ENDIF
IF ((MODE .EQ. 1) .OR. (MODE .EQ. 2)) THEN
IF (EMU(MODE) .EQ. HALF) THEN
TMU = ABS(T)
ELSE
TMU = (T*T)**EMU(MODE)
ENDIF
FP = CP(MODE)
IF (EP(MODE) .NE. ZERO) FP = FP*DX**EP(MODE)
FR = CR(MODE)
IF (ER(MODE) .NE. ZERO) FR = FR*DX**ER(MODE)
PX = PX*TMU/FP
QX = (QX/FR - D(3,MODE)/T**2)*TMU
RX = RX*TMU/FR
ENDIF
TAU = (EV*RX-QX)/PX
OMEGA = SQRT(ABS(TAU))
HOMEGA = H*OMEGA
SCLOG = ZERO
C
C Evaluate the scaled basis functions.
C
IF (HOMEGA .GT. HNDRTH) THEN
IF (TAU .GT. ZERO) THEN
IF (HOMEGA . GT. OVER) THEN
FLAG = -10
RETURN
ENDIF
DPSI = COS(HOMEGA)
PSI = SIN(HOMEGA)/OMEGA
ELSE
SCLOG = HOMEGA
IF (HOMEGA .LT. UNDER) THEN
T = TANH(HOMEGA)
DPSI = ONE/(ONE+T)
PSI = T*DPSI/OMEGA
ELSE
SCLOG = MIN(SCLOG,TWO*UNDER)
DPSI = HALF
PSI = DPSI/OMEGA
ENDIF
ENDIF
ELSE
T = TAU*H*H
DPSI = ONE+T*(T/TWELVE-ONE)/TWO
PSI = H*(ONE+T*(T/TWENTY-ONE)/SIX)
IF (T .LT. ZERO) THEN
T = MAX(ABS(PSI),ABS(DPSI))
SCLOG = LOG(T)
PSI = PSI/T
DPSI = DPSI/T
ENDIF
ENDIF
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE ZZERO(B,C,FB,FC,ABSERR,RELERR,IFLAG,X)
C
DOUBLE PRECISION B,C,FB,FC,ABSERR,RELERR,X(*)
INTEGER IFLAG
C
C ZZERO computes a root of F. The method used is a combination of
C bisection and the secant rule. This code is adapted from one in
C the text "Foundations of Numerical Computing" written by Allen,
C Pruess, and Shampine.
C
C Input parameters:
C B,C = values of X such that F(B)*F(C) .LE. 0.
C FB,FC = values of F at input B and C, resp.
C ABSERR,RELERR = absolute and relative error tolerances. The
C stopping criterion is:
C ABS(B-C) .LE. 2.0*MAX(ABSERR,ABS(B)*RELERR).
C Output parameters:
C B,C = see IFLAG returns.
C FB = value of final residual F(B).
C IFLAG = 0 for normal return; F(B)*F(C) .LT. 0 and the
C stopping criterion is met (or F(B)=0). B always
C satisfies ABS(F(B)) .LE. ABS(F(C)).
C = 1 if too many function evaluations were made; in this version
C 200 are allowed.
C =-2 if F(B)*F(C) is positive on input.
C
C Local variables:
DOUBLE PRECISION A,ACMB,CMB,FA,P,Q,TOL,WIDTH
INTEGER KOUNT,MAXF,MU,NF
C
C Internal constants
C
DOUBLE PRECISION ZERO,ONE,TWO,EIGHT
PARAMETER (ZERO = 0.0, ONE = 1.0, TWO = 2.0, EIGHT = 8.0,
& MAXF = 200)
C
C Initialization.
C
KOUNT = 0
WIDTH = ABS(B-C)
A = C
FA = FC
IF (SIGN(ONE,FA) .EQ. SIGN(ONE,FB)) THEN
IFLAG = -2
RETURN
ENDIF
FC = FA
NF = 2
20 IF (ABS(FC) .LT. ABS(FB)) THEN
C
C Interchange B and C so that ABS(F(B)) .LE. ABS(F(C)).
C
A = B
FA = FB
B = C
FB = FC
C = A
FC = FA
ENDIF
CMB = (C-B)/TWO
ACMB = ABS(CMB)
TOL = MAX(ABSERR,ABS(B)*RELERR)
C
C Test stopping criterion and function count.
C
IF (ACMB .LE. TOL) THEN
IFLAG = 0
RETURN
ENDIF
IF (NF .GE. MAXF) THEN
IFLAG = 1
RETURN
ENDIF
C
C Calculate new iterate implicitly as B+P/Q where we arrange
C P .GE. 0. The implicit form is used to prevent overflow.
C
P = (B-A)*FB
Q = FA-FB
IF (P .LT. ZERO) THEN
P = -P
Q = -Q
ENDIF
C
C Update A; check if reduction in the size of bracketing interval is
C satisfactory. If not, bisect until it is.
C
A = B
FA = FB
KOUNT = KOUNT+1
IF (KOUNT .GE. 4) THEN
IF (EIGHT*ACMB .GE. WIDTH) THEN
B = B+CMB
GOTO 30
ENDIF
KOUNT = 0
WIDTH = ACMB
ENDIF
C
C Test for too small a change.
C
IF (P .LE. ABS(Q)*TOL) THEN
C
C Increment by tolerance.
C
B = B+SIGN(TOL,CMB)
ELSE
C
C Root ought to be between B and (C+B)/2.
C
IF (P .LT. CMB*Q) THEN
C
C Use secant rule.
C
B = B+P/Q
ELSE
C
C Use bisection.
C
B = B+CMB
ENDIF
ENDIF
C
C Have completed computation for new iterate B.
C
30 CALL SHOOT(B,X,MU,FB)
NF = NF+1
IF (ABS(FB) .EQ. ZERO) THEN
IFLAG = 0
C = B
FC = FB
RETURN
ENDIF
IF (SIGN(ONE,FB) .EQ. SIGN(ONE,FC)) THEN
C = A
FC = FA
ENDIF
GOTO 20
END