SLATEC Routines --- DBVSUP ---
*DECK DBVSUP
SUBROUTINE DBVSUP (Y, NROWY, NCOMP, XPTS, NXPTS, A, NROWA, ALPHA,
+ NIC, B, NROWB, BETA, NFC, IGOFX, RE, AE, IFLAG, WORK, NDW,
+ IWORK, NDIW, NEQIVP)
C***BEGIN PROLOGUE DBVSUP
C***PURPOSE Solve a linear two-point boundary value problem using
C superposition coupled with an orthonormalization procedure
C and a variable-step integration scheme.
C***LIBRARY SLATEC
C***CATEGORY I1B1
C***TYPE DOUBLE PRECISION (BVSUP-S, DBVSUP-D)
C***KEYWORDS ORTHONORMALIZATION, SHOOTING,
C TWO-POINT BOUNDARY VALUE PROBLEM
C***AUTHOR Scott, M. R., (SNLA)
C Watts, H. A., (SNLA)
C***DESCRIPTION
C
C **********************************************************************
C
C Subroutine DBVSUP solves a linear two-point boundary-value problem
C of the form
C DY/DX = MATRIX(X,U)*Y(X) + G(X,U)
C A*Y(XINITIAL) = ALPHA , B*Y(XFINAL) = BETA
C
C coupled with the solution of the initial value problem
C
C DU/DX = F(X,U)
C U(XINITIAL) = ETA
C
C **********************************************************************
C ABSTRACT
C The method of solution uses superposition coupled with an
C orthonormalization procedure and a variable-step integration
C scheme. Each time the superposition solutions start to
C lose their numerical linear independence, the vectors are
C reorthonormalized before integration proceeds. The underlying
C principle of the algorithm is then to piece together the
C intermediate (orthogonalized) solutions, defined on the various
C subintervals, to obtain the desired solutions.
C
C **********************************************************************
C INPUT to DBVSUP
C **********************************************************************
C
C NROWY = actual row dimension of Y in calling program.
C NROWY must be .GE. NCOMP
C
C NCOMP = number of components per solution vector.
C NCOMP is equal to number of original differential
C equations. NCOMP = NIC + NFC.
C
C XPTS = desired output points for solution. They must be monotonic.
C XINITIAL = XPTS(1)
C XFINAL = XPTS(NXPTS)
C
C NXPTS = number of output points.
C
C A(NROWA,NCOMP) = boundary condition matrix at XINITIAL
C must be contained in (NIC,NCOMP) sub-matrix.
C
C NROWA = actual row dimension of A in calling program,
C NROWA must be .GE. NIC.
C
C ALPHA(NIC+NEQIVP) = boundary conditions at XINITIAL.
C If NEQIVP .GT. 0 (see below), the boundary
C conditions at XINITIAL for the initial value
C equations must be stored starting in
C position (NIC + 1) of ALPHA.
C Thus, ALPHA(NIC+K) = ETA(K).
C
C NIC = number of boundary conditions at XINITIAL.
C
C B(NROWB,NCOMP) = boundary condition matrix at XFINAL.
C Must be contained in (NFC,NCOMP) sub-matrix.
C
C NROWB = actual row dimension of B in calling program,
C NROWB must be .GE. NFC.
C
C BETA(NFC) = boundary conditions at XFINAL.
C
C NFC = number of boundary conditions at XFINAL.
C
C IGOFX =0 -- The inhomogeneous term G(X) is identically zero.
C =1 -- The inhomogeneous term G(X) is not identically zero.
C (if IGOFX=1, then Subroutine DGVEC (or DUVEC) must be
C supplied).
C
C RE = relative error tolerance used by the integrator.
C (see one of the integrators)
C
C AE = absolute error tolerance used by the integrator.
C (see one of the integrators)
C **NOTE- RE and AE should not both be zero.
C
C IFLAG = a status parameter used principally for output.
C However, for efficient solution of problems which
C are originally defined as COMPLEX*16 valued (but
C converted to double precision systems to use this code),
C the user must set IFLAG=13 on input. See the comment
C below for more information on solving such problems.
C
C WORK(NDW) = floating point array used for internal storage.
C
C NDW = actual dimension of work array allocated by user.
C An estimate for NDW can be computed from the following
C NDW = 130 + NCOMP**2 * (6 + NXPTS/2 + expected number of
C orthonormalizations/8)
C For the disk or tape storage mode,
C NDW = 6 * NCOMP**2 + 10 * NCOMP + 130
C However, when the ADAMS integrator is to be used, the estimates are
C NDW = 130 + NCOMP**2 * (13 + NXPTS/2 + expected number of
C orthonormalizations/8)
C and NDW = 13 * NCOMP**2 + 22 * NCOMP + 130 , respectively.
C
C IWORK(NDIW) = integer array used for internal storage.
C
C NDIW = actual dimension of IWORK array allocated by user.
C An estimate for NDIW can be computed from the following
C NDIW = 68 + NCOMP * (1 + expected number of
C orthonormalizations)
C **NOTE -- the amount of storage required is problem dependent and may
C be difficult to predict in advance. Experience has shown
C that for most problems 20 or fewer orthonormalizations
C should suffice. If the problem cannot be completed with the
C allotted storage, then a message will be printed which
C estimates the amount of storage necessary. In any case, the
C user can examine the IWORK array for the actual storage
C requirements, as described in the output information below.
C
C NEQIVP = number of auxiliary initial value equations being added
C to the boundary value problem.
C **NOTE -- Occasionally the coefficients matrix and/or G may be
C functions which depend on the independent variable X and
C on U, the solution of an auxiliary initial value problem.
C In order to avoid the difficulties associated with
C interpolation, the auxiliary equations may be solved
C simultaneously with the given boundary value problem.
C This initial value problem may be linear or nonlinear.
C See SAND77-1328 for an example.
C
C
C The user must supply subroutines DFMAT, DGVEC, DUIVP and DUVEC,
C when needed (they must be so named), to evaluate the derivatives
C as follows
C
C A. DFMAT must be supplied.
C
C SUBROUTINE DFMAT(X,Y,YP)
C X = independent variable (input to DFMAT)
C Y = dependent variable vector (input to DFMAT)
C YP = DY/DX = derivative vector (output from DFMAT)
C
C Compute the derivatives for the homogeneous problem
C YP(I) = DY(I)/DX = MATRIX(X) * Y(I) , I = 1,...,NCOMP
C
C When (NEQIVP .GT. 0) and matrix is dependent on U as
C well as on X, the following common statement must be
C included in DFMAT
C COMMON /DMLIVP/ NOFST
C for convenience, the U vector is stored at the bottom
C of the Y array. Thus, during any call to DFMAT,
C U(I) is referenced by Y(NOFST + I).
C
C
C Subroutine DBVDER calls DFMAT NFC times to evaluate the
C homogeneous equations and, if necessary, it calls DFMAT
C once in evaluating the particular solution. since X remains
C unchanged in this sequence of calls it is possible to
C realize considerable computational savings for complicated
C and expensive evaluations of the matrix entries. To do this
C the user merely passes a variable, say XS, via common where
C XS is defined in the main program to be any value except
C the initial X. Then the non-constant elements of matrix(x)
C appearing in the differential equations need only be
C computed if X is unequal to XS, whereupon XS is reset to X.
C
C
C B. If NEQIVP .GT. 0 , DUIVP must also be supplied.
C
C SUBROUTINE DUIVP(X,U,UP)
C X = independent variable (input to DUIVP)
C U = dependent variable vector (input to DUIVP)
C UP = DU/DX = derivative vector (output from DUIVP)
C
C Compute the derivatives for the auxiliary initial value eqs
C UP(I) = DU(I)/DX, I = 1,...,NEQIVP.
C
C Subroutine DBVDER calls DUIVP once to evaluate the
C derivatives for the auxiliary initial value equations.
C
C
C C. If NEQIVP = 0 and IGOFX = 1 , DGVEC must be supplied.
C
C SUBROUTINE DGVEC(X,G)
C X = independent variable (input to DGVEC)
C G = vector of inhomogeneous terms G(X) (output from
C DGVEC)
C
C Compute the inhomogeneous terms G(X)
C G(I) = G(X) values for I = 1,...,NCOMP.
C
C Subroutine DBVDER calls DGVEC in evaluating the particular
C solution provided G(X) is not identically zero. Thus, when
C IGOFX=0, the user need not write a DGVEC subroutine. Also,
C the user does not have to bother with the computational
C savings scheme for DGVEC as this is automatically achieved
C via the DBVDER subroutine.
C
C
C D. If NEQIVP .GT. 0 and IGOFX = 1 , DUVEC must be supplied.
C
C SUBROUTINE DUVEC(X,U,G)
C X = independent variable (input to DUVEC)
C U = dependent variable vector from the auxiliary initial
C value problem (input to DUVEC)
C G = array of inhomogeneous terms G(X,U)(output from DUVEC)
C
C Compute the inhomogeneous terms G(X,U)
C G(I) = G(X,U) values for I = 1,...,NCOMP.
C
C Subroutine DBVDER calls DUVEC in evaluating the particular
C solution provided G(X,U) is not identically zero. Thus,
C when IGOFX=0, the user need not write a DUVEC subroutine.
C
C
C
C The following is optional input to DBVSUP to give user more
C flexibility in use of code. See SAND75-0198, SAND77-1328,
C SAND77-1690, SAND78-0522, and SAND78-1501 for more information.
C
C ****CAUTION -- The user must zero out IWORK(1),...,IWORK(15)
C prior to calling DBVSUP. These locations define
C optional input and must be zero unless set to special
C values by the user as described below.
C
C IWORK(1) -- number of orthonormalization points.
C A value need be set only if IWORK(11) = 1
C
C IWORK(9) -- integrator and orthonormalization parameter
C (default value is 1)
C 1 = RUNGE-KUTTA-FEHLBERG code using GRAM-SCHMIDT test.
C 2 = ADAMS code using GRAM-SCHMIDT test.
C
C IWORK(11) -- orthonormalization points parameter
C (default value is 0)
C 0 - orthonormalization points not pre-assigned.
C 1 - orthonormalization points pre-assigned in
C the first IWORK(1) positions of work.
C
C IWORK(12) -- storage parameter
C (default value is 0)
C 0 - all storage in core.
C LUN - homogeneous and inhomogeneous solutions at
C output points and orthonormalization information
C are stored on disk. The logical unit number to
C be used for disk I/O (NTAPE) is set to IWORK(12).
C
C WORK(1),... -- pre-assigned orthonormalization points, stored
C monotonically, corresponding to the direction
C of integration.
C
C
C
C ******************************************************
C *** COMPLEX*16 VALUED PROBLEM ***
C ******************************************************
C **NOTE***
C Suppose the original boundary value problem is NC equations
C of the form
C DW/DX = MAT(X,U)*W(X) + H(X,U)
C R*W(XINITIAL)=GAMMA , S*W(XFINAL)=DELTA
C where all variables are COMPLEX*16 valued. The DBVSUP code can be
C used by converting to a double precision system of size 2*NC. To
C solve the larger dimensioned problem efficiently, the user must
C initialize IFLAG=13 on input and order the vector components
C according to Y(1)=DOUBLE PRECISION(W(1)),...,Y(NC)=DOUBLE
C PRECISION(W(NC)),Y(NC+1)=IMAG(W(1)),...., Y(2*NC)=IMAG(W(NC)).
C Then define
C ...............................................
C . DOUBLE PRECISION(MAT) -IMAG(MAT) .
C MATRIX = . .
C . IMAG(MAT) DOUBLE PRECISION(MAT) .
C ...............................................
C
C The matrices A,B and vectors G,ALPHA,BETA must be defined
C similarly. Further details can be found in SAND78-1501.
C
C
C **********************************************************************
C OUTPUT from DBVSUP
C **********************************************************************
C
C Y(NROWY,NXPTS) = solution at specified output points.
C
C IFLAG Output Values
C =-5 algorithm ,for obtaining starting vectors for the
C special COMPLEX*16 problem structure, was unable to
C obtain the initial vectors satisfying the necessary
C independence criteria.
C =-4 rank of boundary condition matrix A is less than NIC,
C as determined by DLSSUD.
C =-2 invalid input parameters.
C =-1 insufficient number of storage locations allocated for
C WORK or IWORK.
C
C =0 indicates successful solution.
C
C =1 a computed solution is returned but uniqueness of the
C solution of the boundary-value problem is questionable.
C For an eigenvalue problem, this should be treated as a
C successful execution since this is the expected mode
C of return.
C =2 a computed solution is returned but the existence of the
C solution to the boundary-value problem is questionable.
C =3 a nontrivial solution approximation is returned although
C the boundary condition matrix B*Y(XFINAL) is found to be
C nonsingular (to the desired accuracy level) while the
C right hand side vector is zero. To eliminate this type
C of return, the accuracy of the eigenvalue parameter
C must be improved.
C ***NOTE-We attempt to diagnose the correct problem behavior
C and report possible difficulties by the appropriate
C error flag. However, the user should probably resolve
C the problem using smaller error tolerances and/or
C perturbations in the boundary conditions or other
C parameters. This will often reveal the correct
C interpretation for the problem posed.
C
C =13 maximum number of orthonormalizations attained before
C reaching XFINAL.
C =20-flag from integrator (DDERKF or DDEABM) values can
C range from 21 to 25.
C =30 solution vectors form a dependent set.
C
C WORK(1),...,WORK(IWORK(1)) = orthonormalization points
C determined by DBVPOR.
C
C IWORK(1) = number of orthonormalizations performed by DBVPOR.
C
C IWORK(2) = maximum number of orthonormalizations allowed as
C calculated from storage allocated by user.
C
C IWORK(3),IWORK(4),IWORK(5),IWORK(6) give information about
C actual storage requirements for WORK and IWORK
C arrays. In particular,
C required storage for work array is
C IWORK(3) + IWORK(4)*(expected number of orthonormalizations)
C
C required storage for IWORK array is
C IWORK(5) + IWORK(6)*(expected number of orthonormalizations)
C
C IWORK(8) = final value of exponent parameter used in tolerance
C test for orthonormalization.
C
C IWORK(16) = number of independent vectors returned from DMGSBV.
C It is only of interest when IFLAG=30 is obtained.
C
C IWORK(17) = numerically estimated rank of the boundary
C condition matrix defined from B*Y(XFINAL)
C
C **********************************************************************
C
C Necessary machine constants are defined in the Function
C Routine D1MACH. The user must make sure that the values
C set in D1MACH are relevant to the computer being used.
C
C **********************************************************************
C **********************************************************************
C
C***REFERENCES M. R. Scott and H. A. Watts, SUPORT - a computer code
C for two-point boundary-value problems via
C orthonormalization, SIAM Journal of Numerical
C Analysis 14, (1977), pp. 40-70.
C B. L. Darlow, M. R. Scott and H. A. Watts, Modifications
C of SUPORT, a linear boundary value problem solver
C Part I - pre-assigning orthonormalization points,
C auxiliary initial value problem, disk or tape storage,
C Report SAND77-1328, Sandia Laboratories, Albuquerque,
C New Mexico, 1977.
C B. L. Darlow, M. R. Scott and H. A. Watts, Modifications
C of SUPORT, a linear boundary value problem solver
C Part II - inclusion of an Adams integrator, Report
C SAND77-1690, Sandia Laboratories, Albuquerque,
C New Mexico, 1977.
C M. E. Lord and H. A. Watts, Modifications of SUPORT,
C a linear boundary value problem solver Part III -
C orthonormalization improvements, Report SAND78-0522,
C Sandia Laboratories, Albuquerque, New Mexico, 1978.
C H. A. Watts, M. R. Scott and M. E. Lord, Computational
C solution of complex*16 valued boundary problems,
C Report SAND78-1501, Sandia Laboratories,
C Albuquerque, New Mexico, 1978.
C***ROUTINES CALLED DEXBVP, DMACON, XERMSG
C***COMMON BLOCKS DML15T, DML17B, DML18J, DML5MC, DML8SZ
C***REVISION HISTORY (YYMMDD)
C 750601 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890921 Realigned order of variables in certain COMMON blocks.
C (WRB)
C 890921 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900510 Convert XERRWV calls to XERMSG calls, remove some extraneous
C comments. (RWC)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DBVSUP