relibc_openlibm/slatec/dbvpor.f

*DECK DBVPOR
      SUBROUTINE DBVPOR (Y, NROWY, NCOMP, XPTS, NXPTS, A, NROWA, ALPHA,
     +   NIC, B, NROWB, BETA, NFC, IFLAG, Z, MXNON, P, NTP, IP, W, NIV,
     +   YHP, U, V, COEF, S, STOWA, G, WORK, IWORK, NFCC)
C***BEGIN PROLOGUE  DBVPOR
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DBVSUP
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (BVPOR-S, DBVPOR-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C **********************************************************************
C     INPUT to DBVPOR    (items not defined in DBVSUP comments)
C **********************************************************************
C
C     NOPG = 0 -- orthonormalization points not pre-assigned
C          = 1 -- orthonormalization points pre-assigned
C
C     MXNON = maximum number of orthogonalizations allowed.
C
C     NDISK = 0 -- in-core storage
C           = 1 -- disk storage.  Value of NTAPE in data statement
C                  is set to 13.  If another value is desired,
C                  the data statement must be changed.
C
C     INTEG = type of integrator and associated test to be used
C             to determine when to orthonormalize.
C
C             1 -- use GRAM-SCHMIDT test and DDERKF
C             2 -- use GRAM-SCHMIDT test and DDEABM
C
C     TOL = tolerance for allowable error in orthogonalization test.
C
C     NPS = 0 normalize particular solution to unit length at each
C             point of orthonormalization.
C         = 1 do not normalize particular solution.
C
C     NTP = must be .GE. NFC*(NFC+1)/2.
C
C     NFCC = 2*NFC for special treatment of a COMPLEX*16 valued problem
C
C     ICOCO = 0 skip final computations (superposition coefficients
C               and, hence, boundary problem solution)
C           = 1 calculate superposition coefficients and obtain
C               solution to the boundary value problem
C
C **********************************************************************
C     OUTPUT from DBVPOR
C **********************************************************************
C
C     Y(NROWY,NXPTS) = solution at specified output points.
C
C     MXNON = number of orthonormalizations performed by DBVPOR.
C
C     Z(MXNON+1) = locations of orthonormalizations performed by DBVPOR.
C
C     NIV = number of independent vectors returned from DMGSBV. Normally
C           this parameter will be meaningful only when DMGSBV returns
C           with MFLAG = 2.
C
C **********************************************************************
C
C     The following variables are in the argument list because of
C     variable dimensioning.  In general, they contain no information of
C     use to the user.  The amount of storage set aside by the user must
C     be greater than or equal to that indicated by the dimension
C     statements.  For the disk storage mode, NON = 0 and KPTS = 1,
C     while for the in-core storage mode, NON = MXNON and KPTS = NXPTS.
C
C     P(NTP,NON+1)
C     IP(NFCC,NON+1)
C     YHP(NCOMP,NFC+1)  plus an additional column of the length  NEQIVP
C     U(NCOMP,NFC,KPTS)
C     V(NCOMP,KPTS)
C     W(NFCC,NON+1)
C     COEF(NFCC)
C     S(NFC+1)
C     STOWA(NCOMP*(NFC+1)+NEQIVP+1)
C     G(NCOMP)
C     WORK(KKKWS)
C     IWORK(LLLIWS)
C
C **********************************************************************
C     SUBROUTINES used by DBVPOR
C         DLSSUD -- solves an underdetermined system of linear
C                   equations.  This routine is used to get a full
C                   set of initial conditions for integration.
C                   Called by DBVPOR.
C
C         DVECS -- obtains starting vectors for special treatment
C                   of COMPLEX*16 valued problems, called by DBVPOR.
C
C         DRKFAB -- routine which conducts integration using DDERKF or
C                   DDEABM.
C
C         DSTWAY -- storage for backup capability, called by
C                   DBVPOR and DREORT.
C
C         DSTOR1 -- storage at output points, called by DBVPOR,
C                   DRKFAB, DREORT and DSTWAY.
C
C         DDOT -- single precision vector inner product routine,
C                   called by DBVPOR, DCOEF, DLSSUD, DMGSBV,
C                   DBKSOL, DREORT and DPRVEC.
C         ** NOTE **
C         a considerable improvement in speed can be achieved if a
C         machine language version is used for DDOT.
C
C         DCOEF -- computes the superposition constants from the
C                   boundary conditions at XFINAL.
C
C         DBKSOL -- solves an upper triangular set of linear equations.
C
C **********************************************************************
C
C***SEE ALSO  DBVSUP
C***ROUTINES CALLED  DBKSOL, DCOEF, DDOT, DLSSUD, DRKFAB, DSTOR1,
C                    DSTWAY, DVECS
C***COMMON BLOCKS    DML15T, DML18J, 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   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C***END PROLOGUE  DBVPOR
C
      DOUBLE PRECISION DDOT
      INTEGER I, I1, I2, IC, ICOCO, IFLAG, IGOFX, INDPVT, INFO, INHOMO,
     1     INTEG, IRA, ISFLG, ISTKOP, IVP, J,
     2     K, KNSWOT, KOD, KOP, KPTS, KWC, KWD, KWS, KWT, L, LOTJP, M,
     3     MNSWOT, MXNON, MXNOND, N, NCOMP, NCOMP2, NCOMPD, NDISK, NDW,
     4     NEQ, NEQIVP, NFC, NFCC, NFCCD, NFCD, NFCP1, NFCP2, NIC,
     5     NICD, NIV, NN, NON, NOPG, NPS, NROWA, NROWB, NROWY, NSWOT,
     6     NTAPE, NTP, NTPD, NUMORT, NXPTS, NXPTSD,
     7     IP(NFCC,*), IWORK(*)
      DOUBLE PRECISION A(NROWA,*), AE, ALPHA(*), B(NROWB,*),
     1     BETA(*), C, COEF(*), G(*), P(NTP,*), PWCND, PX,
     2     RE, S(*), STOWA(*), TND, TOL, U(NCOMP,NFC,*),
     3     V(NCOMP,*), W(NFCC,*), WORK(*), X, XBEG, XEND, XOP,
     4     XOT, XPTS(*), XSAV, Y(NROWY,*), YHP(NCOMP,*),
     5     Z(*)
C
C     ******************************************************************
C
      COMMON /DML8SZ/ C,XSAV,IGOFX,INHOMO,IVP,NCOMPD,NFCD
      COMMON /DML15T/ PX,PWCND,TND,X,XBEG,XEND,XOT,XOP,INFO(15),ISTKOP,
     1                KNSWOT,KOP,LOTJP,MNSWOT,NSWOT
      COMMON /DML18J/ AE,RE,TOL,NXPTSD,NICD,NOPG,MXNOND,NDISK,NTAPE,
     1                NEQ,INDPVT,INTEG,NPS,NTPD,NEQIVP,NUMORT,NFCCD,
     2                ICOCO
C
C      *****************************************************************
C
C***FIRST EXECUTABLE STATEMENT  DBVPOR
      NFCP1 = NFC + 1
      NUMORT = 0
      C = 1.0D0
C
C     ******************************************************************
C         CALCULATE INITIAL CONDITIONS WHICH SATISFY
C                       A*YH(XINITIAL)=0  AND  A*YP(XINITIAL)=ALPHA.
C         WHEN NFC .NE. NFCC DLSSUD DEFINES VALUES YHP IN A MATRIX OF
C         SIZE (NFCC+1)*NCOMP AND ,HENCE, OVERFLOWS THE STORAGE
C         ALLOCATION INTO THE U ARRAY. HOWEVER, THIS IS OKAY SINCE
C         PLENTY OF SPACE IS AVAILABLE IN U AND IT HAS NOT YET BEEN
C         USED.
C
      NDW = NROWA*NCOMP
      KWS = NDW + NIC + 1
      KWD = KWS + NIC
      KWT = KWD + NIC
      KWC = KWT + NIC
      IFLAG = 0
      CALL DLSSUD(A,YHP(1,NFCC+1),ALPHA,NIC,NCOMP,NROWA,YHP,NCOMP,IFLAG,
     1            1,IRA,0,WORK(1),WORK(NDW+1),IWORK,WORK(KWS),WORK(KWD),
     2            WORK(KWT),ISFLG,WORK(KWC))
      IF (IFLAG .EQ. 1) GO TO 10
         IFLAG = -4
      GO TO 200
   10 CONTINUE
         IF (NFC .NE. NFCC)
     1      CALL DVECS(NCOMP,NFC,YHP,WORK,IWORK,INHOMO,IFLAG)
         IF (IFLAG .EQ. 1) GO TO 20
            IFLAG = -5
         GO TO 190
   20    CONTINUE
C
C           ************************************************************
C               DETERMINE THE NUMBER OF DIFFERENTIAL EQUATIONS TO BE
C               INTEGRATED, INITIALIZE VARIABLES FOR AUXILIARY INITIAL
C               VALUE PROBLEM AND STORE INITIAL CONDITIONS.
C
            NEQ = NCOMP*NFC
            IF (INHOMO .EQ. 1) NEQ = NEQ + NCOMP
            IVP = 0
            IF (NEQIVP .EQ. 0) GO TO 40
               IVP = NEQ
               NEQ = NEQ + NEQIVP
               NFCP2 = NFCP1
               IF (INHOMO .EQ. 1) NFCP2 = NFCP1 + 1
               DO 30 K = 1, NEQIVP
                  YHP(K,NFCP2) = ALPHA(NIC+K)
   30          CONTINUE
   40       CONTINUE
            CALL DSTOR1(U,YHP,V,YHP(1,NFCP1),0,NDISK,NTAPE)
C
C           ************************************************************
C               SET UP DATA FOR THE ORTHONORMALIZATION TESTING PROCEDURE
C               AND SAVE INITIAL CONDITIONS IN CASE A RESTART IS
C               NECESSARY.
C
            NSWOT = 1
            KNSWOT = 0
            LOTJP = 1
            TND = LOG10(10.0D0*TOL)
            PWCND = LOG10(SQRT(TOL))
            X = XBEG
            PX = X
            XOT = XEND
            XOP = X
            KOP = 1
            CALL DSTWAY(U,V,YHP,0,STOWA)
C
C           ************************************************************
C           ******** FORWARD INTEGRATION OF ALL INITIAL VALUE EQUATIONS
C           **********
C           ************************************************************
C
            CALL DRKFAB(NCOMP,XPTS,NXPTS,NFC,IFLAG,Z,MXNON,P,NTP,IP,YHP,
     1                  NIV,U,V,W,S,STOWA,G,WORK,IWORK,NFCC)
            IF (IFLAG .NE. 0 .OR. ICOCO .EQ. 0) GO TO 180
C
C              *********************************************************
C              **************** BACKWARD SWEEP TO OBTAIN SOLUTION
C              *******************
C              *********************************************************
C
C                  CALCULATE SUPERPOSITION COEFFICIENTS AT XFINAL.
C
C                FOR THE DISK STORAGE VERSION, IT IS NOT NECESSARY TO
C                READ  U  AND  V AT THE LAST OUTPUT POINT, SINCE THE
C                LOCAL COPY OF EACH STILL EXISTS.
C
               KOD = 1
               IF (NDISK .EQ. 0) KOD = NXPTS
               I1 = 1 + NFCC*NFCC
               I2 = I1 + NFCC
               CALL DCOEF(U(1,1,KOD),V(1,KOD),NCOMP,NROWB,NFC,NIC,B,
     1                     BETA,COEF,INHOMO,RE,AE,WORK,WORK(I1),
     2                     WORK(I2),IWORK,IFLAG,NFCC)
C
C              *********************************************************
C                  CALCULATE SOLUTION AT OUTPUT POINTS BY RECURRING
C                  BACKWARDS.  AS WE RECUR BACKWARDS FROM XFINAL TO
C                  XINITIAL WE MUST CALCULATE NEW SUPERPOSITION
C                  COEFFICIENTS EACH TIME WE CROSS A POINT OF
C                  ORTHONORMALIZATION.
C
               K = NUMORT
               NCOMP2 = NCOMP/2
               IC = 1
               IF (NFC .NE. NFCC) IC = 2
               DO 170 J = 1, NXPTS
                  KPTS = NXPTS - J + 1
                  KOD = KPTS
                  IF (NDISK .EQ. 1) KOD = 1
   50             CONTINUE
C                 ...EXIT
                     IF (K .EQ. 0) GO TO 120
C                 ...EXIT
                     IF (XEND .GT. XBEG .AND. XPTS(KPTS) .GE. Z(K))
     1                  GO TO 120
C                 ...EXIT
                     IF (XEND .LT. XBEG .AND. XPTS(KPTS) .LE. Z(K))
     1                  GO TO 120
                     NON = K
                     IF (NDISK .EQ. 0) GO TO 60
                        NON = 1
                        BACKSPACE NTAPE
                        READ (NTAPE)
     1                       (IP(I,1), I = 1, NFCC),(P(I,1), I = 1, NTP)
                        BACKSPACE NTAPE
   60                CONTINUE
                     IF (INHOMO .NE. 1) GO TO 90
                        IF (NDISK .EQ. 0) GO TO 70
                           BACKSPACE NTAPE
                           READ (NTAPE) (W(I,1), I = 1, NFCC)
                           BACKSPACE NTAPE
   70                   CONTINUE
                        DO 80 N = 1, NFCC
                           COEF(N) = COEF(N) - W(N,NON)
   80                   CONTINUE
   90                CONTINUE
                     CALL DBKSOL(NFCC,P(1,NON),COEF)
                     DO 100 M = 1, NFCC
                        WORK(M) = COEF(M)
  100                CONTINUE
                     DO 110 M = 1, NFCC
                        L = IP(M,NON)
                        COEF(L) = WORK(M)
  110                CONTINUE
                     K = K - 1
                  GO TO 50
  120             CONTINUE
                  IF (NDISK .EQ. 0) GO TO 130
                     BACKSPACE NTAPE
                     READ (NTAPE)
     1                    (V(I,1), I = 1, NCOMP),
     2                    ((U(I,M,1), I = 1, NCOMP), M = 1, NFC)
                     BACKSPACE NTAPE
  130             CONTINUE
                  DO 140 N = 1, NCOMP
                     Y(N,KPTS) = V(N,KOD)
     1                           + DDOT(NFC,U(N,1,KOD),NCOMP,COEF,IC)
  140             CONTINUE
                  IF (NFC .EQ. NFCC) GO TO 160
                     DO 150 N = 1, NCOMP2
                        NN = NCOMP2 + N
                        Y(N,KPTS) = Y(N,KPTS)
     1                              - DDOT(NFC,U(NN,1,KOD),NCOMP,
     2                                     COEF(2),2)
                        Y(NN,KPTS) = Y(NN,KPTS)
     1                               + DDOT(NFC,U(N,1,KOD),NCOMP,
     2                                      COEF(2),2)
  150                CONTINUE
  160             CONTINUE
  170          CONTINUE
  180       CONTINUE
  190    CONTINUE
  200 CONTINUE
C
C     ******************************************************************
C
      MXNON = NUMORT
      RETURN
      END