relibc_openlibm/slatec/drkfab.f

*DECK DRKFAB
      SUBROUTINE DRKFAB (NCOMP, XPTS, NXPTS, NFC, IFLAG, Z, MXNON, P,
     +   NTP, IP, YHP, NIV, U, V, W, S, STOWA, G, WORK, IWORK, NFCC)
C***BEGIN PROLOGUE  DRKFAB
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DBVSUP
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (RKFAB-S, DRKFAB-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C **********************************************************************
C
C     Subroutine DRKFAB integrates the initial value equations using
C     the variable-step Runge-Kutta-Fehlberg integration scheme or
C     the variable-order Adams method and orthonormalization
C     determined by a linear dependence test.
C
C **********************************************************************
C
C***SEE ALSO  DBVSUP
C***ROUTINES CALLED  DBVDER, DDEABM, DDERKF, DREORT, DSTOR1
C***COMMON BLOCKS    DML15T, DML17B, DML18J, DML8SZ
C***REVISION HISTORY  (YYMMDD)
C   750601  DATE WRITTEN
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  DRKFAB
C
      INTEGER ICOCO, IDID, IFLAG, IGOFX, INDPVT, INFO, INHOMO, INTEG,
     1     IPAR, ISTKOP, IVP, J, JFLAG, JON,
     2     K1, K10, K11, K2, K3, K4, K5, K6, K7, K8, K9, KKKINT,
     3     KKKZPW, KNSWOT, KOD, KOP, KOPP, L1, L2, LLLINT, LOTJP,
     4     MNSWOT, MXNON, MXNOND, NCOMP, NCOMPD, NDISK, NEEDIW, NEEDW,
     5     NEQ, NEQIVP, NFC, NFCC, NFCCD, NFCD, NFCP1, NIC, NIV, NON,
     6     NOPG, NPS, NSWOT, NTAPE, NTP, NTPD, NUMORT, NXPTS, NXPTSD,
     7     IP(NFCC,*), IWORK(*)
      DOUBLE PRECISION AE, C, G(*), P(NTP,*), PWCND, PX, RE,
     1     S(*), STOWA(*), TND, TOL, U(NCOMP,NFC,*),
     2     V(NCOMP,*), W(NFCC,*), WORK(*), X, XBEG, XEND, XOP,
     3     XOT, XPTS(*), XSAV, XXOP, YHP(NCOMP,*), 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,NIC,NOPG,MXNOND,NDISK,NTAPE,NEQ,
     1                INDPVT,INTEG,NPS,NTPD,NEQIVP,NUMORT,NFCCD,
     2                ICOCO
      COMMON /DML17B/ KKKZPW,NEEDW,NEEDIW,K1,K2,K3,K4,K5,K6,K7,K8,K9,
     1                K10,K11,L1,L2,KKKINT,LLLINT
C
      EXTERNAL DBVDER
C
C      *****************************************************************
C       INITIALIZATION OF COUNTERS AND VARIABLES.
C
C     BEGIN BLOCK PERMITTING ...EXITS TO 220
C        BEGIN BLOCK PERMITTING ...EXITS TO 10
C***FIRST EXECUTABLE STATEMENT  DRKFAB
            KOD = 1
            NON = 1
            X = XBEG
            JON = 1
            INFO(1) = 0
            INFO(2) = 0
            INFO(3) = 1
            INFO(4) = 1
            WORK(1) = XEND
C        ...EXIT
            IF (NOPG .EQ. 0) GO TO 10
            INFO(3) = 0
            IF (X .EQ. Z(1)) JON = 2
   10    CONTINUE
         NFCP1 = NFC + 1
C
C        ***************************************************************
C        *****BEGINNING OF INTEGRATION LOOP AT OUTPUT
C        POINTS.******************
C        ***************************************************************
C
         DO 210 KOPP = 2, NXPTS
            KOP = KOPP
            XOP = XPTS(KOP)
            IF (NDISK .EQ. 0) KOD = KOP
C
   20       CONTINUE
C
C              STEP BY STEP INTEGRATION LOOP BETWEEN OUTPUT POINTS.
C
C              BEGIN BLOCK PERMITTING ...EXITS TO 190
C                 BEGIN BLOCK PERMITTING ...EXITS TO 30
                     XXOP = XOP
C                 ...EXIT
                     IF (NOPG .EQ. 0) GO TO 30
                     IF (XEND .GT. XBEG .AND. XOP .GT. Z(JON))
     1                  XXOP = Z(JON)
                     IF (XEND .LT. XBEG .AND. XOP .LT. Z(JON))
     1                  XXOP = Z(JON)
   30             CONTINUE
C
C                 ******************************************************
   40             CONTINUE
C                    BEGIN BLOCK PERMITTING ...EXITS TO 170
                        GO TO (50,60), INTEG
C                       DDERKF INTEGRATOR
C
   50                   CONTINUE
                           CALL DDERKF(DBVDER,NEQ,X,YHP,XXOP,INFO,RE,AE,
     1                                 IDID,WORK,KKKINT,IWORK,LLLINT,G,
     2                                 IPAR)
                        GO TO 70
C                       DDEABM INTEGRATOR
C
   60                   CONTINUE
                           CALL DDEABM(DBVDER,NEQ,X,YHP,XXOP,INFO,RE,AE,
     1                                 IDID,WORK,KKKINT,IWORK,LLLINT,G,
     2                                 IPAR)
   70                   CONTINUE
                        IF (IDID .GE. 1) GO TO 80
                           INFO(1) = 1
C                    ......EXIT
                           IF (IDID .EQ. -1) GO TO 170
                           IFLAG = 20 - IDID
C     .....................EXIT
                           GO TO 220
   80                   CONTINUE
C
C                       ************************************************
C                           GRAM-SCHMIDT ORTHOGONALIZATION TEST FOR
C                           ORTHONORMALIZATION (TEMPORARILY USING U AND
C                           V IN THE TEST)
C
                        IF (NOPG .EQ. 0) GO TO 100
                           IF (XXOP .EQ. Z(JON)) GO TO 90
C
C                             ******************************************
C                                 CONTINUE INTEGRATION IF WE ARE NOT AT
C                                 AN OUTPUT POINT.
C
C           ..................EXIT
                              IF (IDID .NE. 1) GO TO 200
C                    .........EXIT
                              GO TO 170
   90                      CONTINUE
                           JFLAG = 2
                        GO TO 110
  100                   CONTINUE
                           JFLAG = 1
                           IF (INHOMO .EQ. 3 .AND. X .EQ. XEND)
     1                        JFLAG = 3
  110                   CONTINUE
C
                        IF (NDISK .EQ. 0) NON = NUMORT + 1
                        CALL DREORT(NCOMP,U(1,1,KOD),V(1,KOD),YHP,NIV,
     1                              W(1,NON),S,P(1,NON),IP(1,NON),STOWA,
     2                              JFLAG)
C
                        IF (JFLAG .NE. 30) GO TO 120
                           IFLAG = 30
C     .....................EXIT
                           GO TO 220
  120                   CONTINUE
C
                        IF (JFLAG .NE. 10) GO TO 130
                           XOP = XPTS(KOP)
                           IF (NDISK .EQ. 0) KOD = KOP
C              ............EXIT
                           GO TO 190
  130                   CONTINUE
C
                        IF (JFLAG .EQ. 0) GO TO 140
C
C                          *********************************************
C                              CONTINUE INTEGRATION IF WE ARE NOT AT AN
C                              OUTPUT POINT.
C
C           ...............EXIT
                           IF (IDID .NE. 1) GO TO 200
C                    ......EXIT
                           GO TO 170
  140                   CONTINUE
C
C                       ************************************************
C                           STORE ORTHONORMALIZED VECTORS INTO SOLUTION
C                           VECTORS.
C
                        IF (NUMORT .LT. MXNON) GO TO 150
                        IF (X .EQ. XEND) GO TO 150
                           IFLAG = 13
C     .....................EXIT
                           GO TO 220
  150                   CONTINUE
C
                        NUMORT = NUMORT + 1
                        CALL DSTOR1(YHP,U(1,1,KOD),YHP(1,NFCP1),
     1                              V(1,KOD),1,NDISK,NTAPE)
C
C                       ************************************************
C                           STORE ORTHONORMALIZATION INFORMATION,
C                           INITIALIZE INTEGRATION FLAG, AND CONTINUE
C                           INTEGRATION TO THE NEXT ORTHONORMALIZATION
C                           POINT OR OUTPUT POINT.
C
                        Z(NUMORT) = X
                        IF (INHOMO .EQ. 1 .AND. NPS .EQ. 0)
     1                     C = S(NFCP1)*C
                        IF (NDISK .EQ. 0) GO TO 160
                           IF (INHOMO .EQ. 1)
     1                        WRITE (NTAPE) (W(J,1), J = 1, NFCC)
                           WRITE (NTAPE)
     1                           (IP(J,1), J = 1, NFCC),
     2                           (P(J,1), J = 1, NTP)
  160                   CONTINUE
                        INFO(1) = 0
                        JON = JON + 1
C                 ......EXIT
                        IF (NOPG .EQ. 1 .AND. X .NE. XOP) GO TO 180
C
C                       ************************************************
C                           CONTINUE INTEGRATION IF WE ARE NOT AT AN
C                           OUTPUT POINT.
C
C           ............EXIT
                        IF (IDID .NE. 1) GO TO 200
  170                CONTINUE
                  GO TO 40
  180             CONTINUE
  190          CONTINUE
            GO TO 20
  200       CONTINUE
C
C           STORAGE OF HOMOGENEOUS SOLUTIONS IN U AND THE PARTICULAR
C           SOLUTION IN V AT THE OUTPUT POINTS.
C
            CALL DSTOR1(U(1,1,KOD),YHP,V(1,KOD),YHP(1,NFCP1),0,NDISK,
     1                  NTAPE)
  210    CONTINUE
C        ***************************************************************
C        ***************************************************************
C
         IFLAG = 0
  220 CONTINUE
      RETURN
      END