relibc_openlibm/slatec/dbvsup.f

*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
C **********************************************************************
C
      INTEGER ICOCO, IFLAG, IGOFX, IGOFXD, INDPVT, INFO, INHOMO, INTEG,
     1     IS, ISTKOP, IVP, IWORK(*), J, K, K1, K10, K11, K2,
     2     K3, K4, K5, K6, K7, K8, K9, KKKCOE, KKKCOF, KKKG, KKKINT,
     3     KKKS, KKKSTO, KKKSUD, KKKSVC, KKKU, KKKV, KKKWS, KKKYHP,
     4     KKKZPW, KNSWOT, KOP, KPTS, L1, L2, LLLCOF, LLLINT, LLLIP,
     5     LLLIWS, LLLSUD, LLLSVC, LOTJP, LPAR, MNSWOT,
     6     MXNON, MXNONI, MXNONR, NCOMP, NCOMPD, NDEQ, NDISK, NDIW,
     7     NDW, NEEDIW, NEEDW, NEQ, NEQIVD, NEQIVP, NFC, NFCC,
     8     NFCD, NIC, NICD, NITEMP, NON, NOPG, NPS, NROWA, NROWB,
     9     NROWY, NRTEMP, NSWOT, NTAPE, NTP, NUMORT, NXPTS, NXPTSD,
     1     NXPTSM
      DOUBLE PRECISION A(NROWA,*), AE, AED, ALPHA(*),
     1     B(NROWB,*), BETA(*), C, EPS, FOURU, PWCND, PX, RE,
     2     RED, SQOVFL, SRU, TND, TOL, TWOU, URO, WORK(NDW), X, XBEG,
     3     XEND, XOP, XOT, XPTS(*), XSAV, Y(NROWY,*)
      CHARACTER*8 XERN1, XERN2, XERN3, XERN4
C
C     ******************************************************************
C         THE COMMON BLOCK BELOW IS USED TO COMMUNICATE WITH SUBROUTINE
C         DBVDER.  THE USER SHOULD NOT ALTER OR USE THIS COMMON BLOCK IN
C         THE CALLING PROGRAM.
C
      COMMON /DML8SZ/ C,XSAV,IGOFXD,INHOMO,IVP,NCOMPD,NFCD
C
C     ******************************************************************
C         THESE COMMON BLOCKS AID IN REDUCING THE NUMBER OF SUBROUTINE
C         ARGUMENTS PREVALENT IN THIS MODULAR STRUCTURE
C
      COMMON /DML18J/ AED,RED,TOL,NXPTSD,NICD,NOPG,MXNON,NDISK,NTAPE,
     1                NEQ,INDPVT,INTEG,NPS,NTP,NEQIVD,NUMORT,NFCC,
     2                ICOCO
      COMMON /DML17B/ KKKZPW,NEEDW,NEEDIW,K1,K2,K3,K4,K5,K6,K7,K8,K9,
     1                K10,K11,L1,L2,KKKINT,LLLINT
C
C     ******************************************************************
C         THIS COMMON BLOCK IS USED IN SUBROUTINES DBVSUP,DBVPOR,DRKFAB,
C         DREORT, AND DSTWAY. IT CONTAINS INFORMATION NECESSARY
C         FOR THE ORTHONORMALIZATION TESTING PROCEDURE AND A BACKUP
C         RESTARTING CAPABILITY.
C
      COMMON /DML15T/ PX,PWCND,TND,X,XBEG,XEND,XOT,XOP,INFO(15),ISTKOP,
     1                KNSWOT,KOP,LOTJP,MNSWOT,NSWOT
C
C     ******************************************************************
C         THIS COMMON BLOCK CONTAINS THE MACHINE DEPENDENT PARAMETERS
C         USED BY THE CODE
C
      COMMON /DML5MC/ URO,SRU,EPS,SQOVFL,TWOU,FOURU,LPAR
C
C      *****************************************************************
C          SET UP MACHINE DEPENDENT CONSTANTS.
C
C***FIRST EXECUTABLE STATEMENT  DBVSUP
                        CALL DMACON
C
C                       ************************************************
C                           TEST FOR INVALID INPUT
C
                        IF (NROWY .LT. NCOMP) GO TO 80
                        IF (NCOMP .NE. NIC + NFC) GO TO 80
                        IF (NXPTS .LT. 2) GO TO 80
                        IF (NIC .LE. 0) GO TO 80
                        IF (NROWA .LT. NIC) GO TO 80
                        IF (NFC .LE. 0) GO TO 80
                        IF (NROWB .LT. NFC) GO TO 80
                        IF (IGOFX .LT. 0 .OR. IGOFX .GT. 1) GO TO 80
                        IF (RE .LT. 0.0D0) GO TO 80
                        IF (AE .LT. 0.0D0) GO TO 80
                        IF (RE .EQ. 0.0D0 .AND. AE .EQ. 0.0D0) GO TO 80
C                          BEGIN BLOCK PERMITTING ...EXITS TO 70
                              IS = 1
                              IF (XPTS(NXPTS) .LT. XPTS(1)) IS = 2
                              NXPTSM = NXPTS - 1
                              DO 30 K = 1, NXPTSM
                                 IF (IS .EQ. 2) GO TO 10
C                          .........EXIT
                                    IF (XPTS(K+1) .LE. XPTS(K)) GO TO 70
                                 GO TO 20
   10                            CONTINUE
C                          .........EXIT
                                    IF (XPTS(K) .LE. XPTS(K+1)) GO TO 70
   20                            CONTINUE
   30                         CONTINUE
C
C                             ******************************************
C                                 CHECK FOR DISK STORAGE
C
                              KPTS = NXPTS
                              NDISK = 0
                              IF (IWORK(12) .EQ. 0) GO TO 40
                                 NTAPE = IWORK(12)
                                 KPTS = 1
                                 NDISK = 1
   40                         CONTINUE
C
C                             ******************************************
C                                 SET INTEG PARAMETER ACCORDING TO
C                                 CHOICE OF INTEGRATOR.
C
                              INTEG = 1
                              IF (IWORK(9) .EQ. 2) INTEG = 2
C
C                             ******************************************
C                                 COMPUTE INHOMO
C
C                 ............EXIT
                              IF (IGOFX .EQ. 1) GO TO 100
                              DO 50 J = 1, NIC
C                 ...............EXIT
                                 IF (ALPHA(J) .NE. 0.0D0) GO TO 100
   50                         CONTINUE
                              DO 60 J = 1, NFC
C                    ............EXIT
                                 IF (BETA(J) .NE. 0.0D0) GO TO 90
   60                         CONTINUE
                              INHOMO = 3
C              ...............EXIT
                              GO TO 110
   70                      CONTINUE
   80                   CONTINUE
                        IFLAG = -2
C     ..................EXIT
                        GO TO 220
   90                CONTINUE
                     INHOMO = 2
C              ......EXIT
                     GO TO 110
  100             CONTINUE
                  INHOMO = 1
  110          CONTINUE
C
C              *********************************************************
C                  TO TAKE ADVANTAGE OF THE SPECIAL STRUCTURE WHEN
C                  SOLVING A COMPLEX*16 VALUED PROBLEM,WE INTRODUCE
C                  NFCC=NFC WHILE CHANGING THE INTERNAL VALUE OF NFC
C
               NFCC = NFC
               IF (IFLAG .EQ. 13) NFC = NFC/2
C
C              *********************************************************
C                  DETERMINE NECESSARY STORAGE REQUIREMENTS
C
C              FOR BASIC ARRAYS IN DBVPOR
               KKKYHP = NCOMP*(NFC + 1) + NEQIVP
               KKKU = NCOMP*NFC*KPTS
               KKKV = NCOMP*KPTS
               KKKCOE = NFCC
               KKKS = NFC + 1
               KKKSTO = NCOMP*(NFC + 1) + NEQIVP + 1
               KKKG = NCOMP
C
C              FOR ORTHONORMALIZATION RELATED MATTERS
               NTP = (NFCC*(NFCC + 1))/2
               KKKZPW = 1 + NTP + NFCC
               LLLIP = NFCC
C
C              FOR ADDITIONAL REQUIRED WORK SPACE
C                (DLSSUD)
               KKKSUD = 4*NIC + (NROWA + 1)*NCOMP
               LLLSUD = NIC
C              (DVECS)
               KKKSVC = 1 + 4*NFCC + 2*NFCC**2
               LLLSVC = 2*NFCC
C
               NDEQ = NCOMP*NFC + NEQIVP
               IF (INHOMO .EQ. 1) NDEQ = NDEQ + NCOMP
               GO TO (120,130), INTEG
C              (DDERKF)
  120          CONTINUE
                  KKKINT = 33 + 7*NDEQ
                  LLLINT = 34
               GO TO 140
C              (DDEABM)
  130          CONTINUE
                  KKKINT = 130 + 21*NDEQ
                  LLLINT = 51
  140          CONTINUE
C
C              (COEF)
               KKKCOF = 5*NFCC + NFCC**2
               LLLCOF = 3 + NFCC
C
               KKKWS = MAX(KKKSUD,KKKSVC,KKKINT,KKKCOF)
               LLLIWS = MAX(LLLSUD,LLLSVC,LLLINT,LLLCOF)
C
               NEEDW = KKKYHP + KKKU + KKKV + KKKCOE + KKKS + KKKSTO
     1                 + KKKG + KKKZPW + KKKWS
               NEEDIW = 17 + LLLIP + LLLIWS
C              *********************************************************
C                  COMPUTE THE NUMBER OF POSSIBLE ORTHONORMALIZATIONS
C                  WITH THE ALLOTTED STORAGE
C
               IWORK(3) = NEEDW
               IWORK(4) = KKKZPW
               IWORK(5) = NEEDIW
               IWORK(6) = LLLIP
               NRTEMP = NDW - NEEDW
               NITEMP = NDIW - NEEDIW
C           ...EXIT
               IF (NRTEMP .LT. 0) GO TO 180
C           ...EXIT
               IF (NITEMP .LT. 0) GO TO 180
C
               IF (NDISK .EQ. 0) GO TO 150
                  NON = 0
                  MXNON = NRTEMP
               GO TO 160
  150          CONTINUE
C
                  MXNONR = NRTEMP/KKKZPW
                  MXNONI = NITEMP/LLLIP
                  MXNON = MIN(MXNONR,MXNONI)
                  NON = MXNON
  160          CONTINUE
C
               IWORK(2) = MXNON
C
C              *********************************************************
C                  CHECK FOR PRE-ASSIGNED ORTHONORMALIZATION POINTS
C
               NOPG = 0
C        ......EXIT
               IF (IWORK(11) .NE. 1) GO TO 210
               IF (MXNON .LT. IWORK(1)) GO TO 170
                  NOPG = 1
                  MXNON = IWORK(1)
                  WORK(MXNON+1) = 2.0D0*XPTS(NXPTS) - XPTS(1)
C        .........EXIT
                  GO TO 210
  170          CONTINUE
  180       CONTINUE
C
            IFLAG = -1
      IF (NDISK .NE. 1) THEN
         WRITE (XERN1, '(I8)') NEEDW
         WRITE (XERN2, '(I8)') KKKZPW
         WRITE (XERN3, '(I8)') NEEDIW
         WRITE (XERN4, '(I8)') LLLIP
         CALL XERMSG ('SLATEC', 'DBVSUP',
     *      'REQUIRED STORAGE FOR WORK ARRAY IS '  // XERN1 // ' + ' //
     *      XERN2 // '*(EXPECTED NUMBER OF ORTHONORMALIZATIONS) $$'  //
     *      'REQUIRED STORAGE FOR IWORK ARRAY IS ' // XERN3 // ' + ' //
     *      XERN4 // '*(EXPECTED NUMBER OF ORTHONORMALIZATIONS)', 1, 0)
      ELSE
         WRITE (XERN1, '(I8)') NEEDW
         WRITE (XERN2, '(I8)') NEEDIW
         CALL XERMSG ('SLATEC', 'DBVSUP',
     *      'REQUIRED STORAGE FOR WORK ARRAY IS '  // XERN1 //
     *      ' + NUMBER OF ORTHONOMALIZATIONS. $$'  //
     *      'REQUIRED STORAGE FOR IWORK ARRAY IS ' // XERN2, 1, 0)
      ENDIF
      RETURN
C
C        ***************************************************************
C            ALLOCATE STORAGE FROM WORK AND IWORK ARRAYS
C
C         (Z)
  210    K1 = 1 + (MXNON + 1)
C        (P)
         K2 = K1 + NTP*(NON + 1)
C        (W)
         K3 = K2 + NFCC*(NON + 1)
C        (YHP)
         K4 = K3 + KKKYHP
C        (U)
         K5 = K4 + KKKU
C        (V)
         K6 = K5 + KKKV
C        (COEF)
         K7 = K6 + KKKCOE
C        (S)
         K8 = K7 + KKKS
C        (STOWA)
         K9 = K8 + KKKSTO
C        (G)
         K10 = K9 + KKKG
         K11 = K10 + KKKWS
C                  REQUIRED ADDITIONAL DOUBLE PRECISION WORK SPACE
C                  STARTS AT WORK(K10) AND EXTENDS TO WORK(K11-1)
C
C           FIRST 17 LOCATIONS OF IWORK ARE USED FOR OPTIONAL
C           INPUT AND OUTPUT ITEMS
C        (IP)
         L1 = 18 + NFCC*(NON + 1)
         L2 = L1 + LLLIWS
C                   REQUIRED INTEGER WORK SPACE STARTS AT IWORK(L1)
C                   AND EXTENDS TO IWORK(L2-1)
C
C        ***************************************************************
C            SET INDICATOR FOR NORMALIZATION OF PARTICULAR SOLUTION
C
         NPS = 0
         IF (IWORK(10) .EQ. 1) NPS = 1
C
C        ***************************************************************
C            SET PIVOTING PARAMETER
C
         INDPVT = 0
         IF (IWORK(15) .EQ. 1) INDPVT = 1
C
C        ***************************************************************
C            SET OTHER COMMON BLOCK PARAMETERS
C
         NFCD = NFC
         NCOMPD = NCOMP
         IGOFXD = IGOFX
         NXPTSD = NXPTS
         NICD = NIC
         RED = RE
         AED = AE
         NEQIVD = NEQIVP
         MNSWOT = 20
         IF (IWORK(13) .EQ. -1) MNSWOT = MAX(1,IWORK(14))
         XBEG = XPTS(1)
         XEND = XPTS(NXPTS)
         XSAV = XEND
         ICOCO = 1
         IF (INHOMO .EQ. 3 .AND. NOPG .EQ. 1) WORK(MXNON+1) = XEND
C
C        ***************************************************************
C
         CALL DEXBVP(Y,NROWY,XPTS,A,NROWA,ALPHA,B,NROWB,BETA,IFLAG,WORK,
     1               IWORK)
         NFC = NFCC
         IWORK(17) = IWORK(L1)
  220 CONTINUE
      RETURN
      END