# dbvsup.f

          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