*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