      PROGRAM PITCON_3DG_N2

C PROGRAMME CALCULANT L'ARBRE DE BIFURCATION DU
C MODELE DE CHARNEY ET DE VORE DANS LE SYSTEME
C A 6 DEGRES DE LIBERTES.
C
C LA METHODE UTILISE EST UNE METHODE DE CONTINUATION,
C CE PROGRAMME UTILISE DES SOUS-PROGRAMMES TIREES
C DE LA BIBLIOTHEQUE SLATEC. 
C
      REAL GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,DELN1,DELN2,EPSN
      REAL K01,K02,KN1,KN2,BETN1,BETN2,H01,H02,HN1,HN2
      COMMON/DV2/GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,DELN1,DELN2,EPSN
     .          ,K01,K02,KN1,KN2,BETN1,BETN2,H01,H02,HN1,HN2 


      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      PARAMETER (NVAR=7)
      PARAMETER (LIW=NVAR+29)
      PARAMETER (LRW=29+(6+NVAR)*NVAR)
C
      EXTERNAL  FXROTH
      EXTERNAL  DFROTH
      EXTERNAL  CDV2
      EXTERNAL  DCDV2
      EXTERNAL  DENSLV
      EXTERNAL  PITCON
C
      REAL      FPAR(1)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(1)
      INTEGER   IWORK(LIW)
      INTEGER   J
      CHARACTER NAME*12
      REAL      RWORK(LRW)
      REAL      XR(NVAR)
      REAL      FJAC(NVAR,NVAR)
      REAL      ARR(NVAR-1,NVAR-1),WR(NVAR-1),WI(NVAR-1)
      INTEGER   MATZ,IERR,IV1(NVAR-1)
      REAL      Z(NVAR-1,NVAR-1),FV1(NVAR-1)

C  FICHIERS DE SORTIES POUR TRACER AVEC XMGRACE:

      OPEN(11,FILE='../xmgr/pitcon_6dg_all.dat',form='formatted')
      OPEN(12,FILE='../xmgr/pitcon_6dg_sta.dat',form='formatted')
      OPEN(13,FILE='../xmgr/pitcon_6dg_pit.dat',form='formatted')
      OPEN(14,FILE='../xmgr/pitcon_6dg_hop.dat',form='formatted')
      GAMN1=  2.4008434   
      GAMN2=  1.92067472 
      GAMN3=  3.84134936  
      ALPN1=  1.66954923 
      ALPN2=  1.53617537  
      DELN1=  2.67127872 
      DELN2=  0.438907236 
      EPSN =  2.4250555  
      K01  =  0.0057075154 
      K02  =0.00210434338
      KN1  =  0.00869251415 
      KN2  =0.00914068893
      BETN1=  0.0869251415 
      BETN2= 0.0571293049
      H01  =  0.137028514 
      H02  = 0.0808351859
      HN1  =  0.0417387309 
      HN2  =  0.0438907221
C
C  Set work arrays to zero:
C
      DO 10 I=1,LIW
        IWORK(I)=0
10      CONTINUE
      DO 20 I=1,LRW
        RWORK(I)=0.0
20      CONTINUE
C
C  Set some entries of work arrays.
C
C  IWORK(1)=0 ; This is a startup
C  IWORK(2)=2 ; Use X(2) for initial parameter
C  IWORK(3)=0 ; Program may choose parameter index
C  IWORK(4)=0 ; Update jacobian every newton step
C  IWORK(5)=3 ; Seek target values for X(3)
C  IWORK(6)=1 ; Seek limit points in X(1)
C  IWORK(7)=1 ; Control amount of output.
C  IWORK(8)=6 ; Output unit
C  IWORK(9)=2 ; Jacobian choice.
C
      IWORK(1)=0
      IWORK(2)=7
      IWORK(3)=0
      IWORK(4)=0
      IWORK(5)=1
      IWORK(6)=1
      IWORK(7)=3
      IWORK(8)=6
      IWORK(9)=0
      IWORK(17)=1000
C
C  RWORK(1)=0.00001; Absolute error tolerance
C  RWORK(2)=0.00001; Relative error tolerance
C  RWORK(3)=0.01   ; Minimum stepsize
C  RWORK(4)=20.0   ; Maximum stepsize
C  RWORK(5)=0.3    ; Starting stepsize
C  RWORK(6)=1.0    ; Starting direction
C  RWORK(7)=1.0    ; Target value (Seek solution with X(3)=1)
C
      RWORK(1)=0.000001
      RWORK(2)=0.000001
      RWORK(3)=0.0001
      RWORK(4)=20.0
      RWORK(4)=0.001 
      RWORK(5)=0.0005
      RWORK(6)=1.0
      RWORK(7)=1.0
C
C  Set starting point.
C

c     XR(1)= 0.5638E-01  
c     XR(2)=-0.3669E-01
c     XR(3)=-0.4050E-02 
c     XR(4)= 0.9460E-03  
c     XR(5)=-0.1680E-02
c     XR(6)=-0.2512E-02
c     XR(7)=0.105
C
      XR(1)= 0.           
      XR(2)= 0.           
      XR(3)= 0.           
      XR(4)= 0.           
      XR(5)= 0.           
      XR(6)= 0.           
      XR(7)= 0.           

      WRITE(6,*)' '
      WRITE(6,*)'PCPRB1:'
      WRITE(6,*)' '
      WRITE(6,*)'PITCON sample program.'
      WRITE(6,*)'CdeV With N=2 '   
      WRITE(6,*)'6 equations, 7 variables.'
      WRITE(6,*)' '
      WRITE(6,*)'Step  Type of Point     '//
     *'X(1)         X(2)         X(3)'
      WRITE(6,*)' '
      I=0
      NAME='Start point  '
      WRITE(6,'(1X,I4,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR)
C
      DO 30 ITER=1,100000
C
        CALL PITCON(DCDV2,FPAR,CDV2,IERROR,IPAR,IWORK,LIW,
     *  NVAR,RWORK,LRW,XR,DENSLV)
C 
        IF(IWORK(1).EQ.1)THEN
          NAME='Corrected    '
        ELSEIF(IWORK(1).EQ.2)THEN
          NAME='Continuation '
        ELSEIF(IWORK(1).EQ.3)THEN
          NAME='Target point '
        ELSEIF(IWORK(1).EQ.4)THEN
          NAME='Limit point  '
          ENDIF
C
        WRITE(6,'(1X,I4,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR)
C
        IF(IWORK(1).EQ.3)THEN
          WRITE(6,*)' '
          WRITE(6,*)'We have reached the point we wanted.'
          WRITE(6,*)'The code may stop now.'
          STOP
          ENDIF
C
        IF(IERROR.NE.0)THEN
          WRITE(6,*)' '
          WRITE(6,*)'An error occurred.'
          WRITE(6,*)'We will terminate the computation now.'
          STOP
          ENDIF
 
C   TEST DE LA STABILITE DE LA SOLUTION:

      CALL DCDV2(NVAR,FPAR,IPAR,XR,FJAC,IERROR)

      DO I=1,NVAR-1
      DO J=1,NVAR-1
      ARR(I,J)=FJAC(I,J)
      ENDDO
      ENDDO

      MATZ=0
      CALL RG (NVAR-1,NVAR-1, ARR, WR, WI, MATZ, Z, IV1, FV1, IERR)
      INSTA=0
      DO I=1,NVAR-1
         IF(WR(I).GT.0..AND.WI(I).NE.0.)INSTA=2
         IF(WR(I).GT.0.AND.INSTA.NE.2) INSTA=1
       ENDDO
                      WRITE(11,301)XR(NVAR)/BETN1*ALPN1
     C                ,XR(1)/BETN1*ALPN1
     C                ,XR(2)*H01/K01,XR(3)*H01/K01       
       IF(INSTA.EQ.0) WRITE(12,301)XR(NVAR)/BETN1*ALPN1
     C                ,XR(1)/BETN1*ALPN1
     C                ,XR(2)*H01/K01,XR(3)*H01/K01       
       IF(INSTA.EQ.1) WRITE(13,301)XR(NVAR)/BETN1*ALPN1
     C                ,XR(1)/BETN1*ALPN1
     C                ,XR(2)*H01/K01,XR(3)*H01/K01       
       IF(INSTA.EQ.2) WRITE(14,301)XR(NVAR)/BETN1*ALPN
     C                1,XR(1)/BETN1*ALPN1
     C                ,XR(2)*H01/K01,XR(3)*H01/K01
       
301    FORMAT(7(1x,E15.8))

30      CONTINUE
C
      WRITE(6,*)' '
      WRITE(6,*)'PITCON did not reach the point of interest.'
      STOP
      END

      SUBROUTINE FXROTH(NVAR,FPAR,IPAR,X,F,IERROR)
C
C  Evaluate the function at X.
C
C  ( X1 - ((X2-5.0)*X2+2.0)*X2 - 13.0 + 34.0*(X3-1.0)  )
C  ( X1 + ((X2+1.0)*X2-14.0)*X2 - 29.0 + 10.0*(X3-1.0) )
C
      REAL      F(*)
      REAL      FPAR(*)
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      X(NVAR)
C
      F(1)=X(1)
     *     -((X(2)-5.0)*X(2)+2.0)*X(2)
     *     -13.0
     *     +34.0*(X(3)-1.0)
C
      F(2)=X(1)
     *    +((X(2)+1.0)*X(2)-14.0)*X(2)
     *    -29.0
     *    +10.0*(X(3)-1.0)
C
      RETURN
      END

      SUBROUTINE DFROTH(NVAR,FPAR,IPAR,X,FJAC,IERROR)
C
C  Evaluate the Jacobian:
C
C  ( 1.0   (-3.0*X(2)+10.0)*X(2)-2.0      34.0   )
C  ( 1.0   (3.0*X(2)+2.0)*X(2)-14.0       10.0   )
C
      REAL      FJAC(NVAR,NVAR)
      REAL      FPAR(*)
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      X(NVAR)
C
      FJAC(1,1)=1.0
      FJAC(1,2)=(-3.0*X(2)+10.0)*X(2)-2.0
      FJAC(1,3)=34.0
C
      FJAC(2,1)=1.0
      FJAC(2,2)=(3.0*X(2)+2.0)*X(2)-14.0
      FJAC(2,3)=10.0
C
      RETURN
      END




      SUBROUTINE CDV2(NVAR,FPAR,IPAR,X,F,IERROR)
C
C  Evaluate the function at X.
C
C  ( X1 - ((X2-5.0)*X2+2.0)*X2 - 13.0 + 34.0*(X3-1.0)  )
C  ( X1 + ((X2+1.0)*X2-14.0)*X2 - 29.0 + 10.0*(X3-1.0) )
C
      REAL      F(*)
      REAL      FPAR(*)
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      X(NVAR)

      REAL GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,DELN1,DELN2,EPSN
      REAL K01,K02,KN1,KN2,BETN1,BETN2,H01,H02,HN1,HN2
      COMMON/DV2/GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,DELN1,DELN2,EPSN
     .          ,K01,K02,KN1,KN2,BETN1,BETN2,H01,H02,HN1,HN2 

      F(1)=K01*(X(7)-X(1))+H01*X(3)

      F(2)=     -KN1*X(2)             +BETN1*X(3)
     *           -ALPN1*X(1)*X(3)
     *           -DELN1*X(4)*X(6)

      F(3)=     -KN1*X(3) -HN1*X(1)-BETN1*X(2)
     *           +ALPN1*X(1)*X(2)
     *           +DELN1*X(4)*X(5)

      F(4)=     -K02*X(4)+H02*X(6)
     *           +EPSN*(X(2)*X(6)-X(3)*X(5))

      F(5)=     -KN2*X(5)            +BETN2*X(6)
     *           -ALPN2*X(1)*X(6)
     *           -DELTN2*X(4)*X(3)

      F(6)=     -KN2*X(6)-HN2*X(4)-BETN2*X(5)
     *           +ALPN2*X(1)*X(5)
     *           +DELTN2*X(4)*X(2)

      RETURN
      END

      SUBROUTINE DCDV2(NVAR,FPAR,IPAR,X,FJAC,IERROR)

      REAL GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,DELN1,DELN2,EPSN
      REAL K01,K02,KN1,KN2,BETN1,BETN2,H01,H02,HN1,HN2
      COMMON/DV2/GAMN1,GAMN2,GAMN3,ALPN1,ALPN2,DELN1,DELN2,EPSN
     .          ,K01,K02,KN1,KN2,BETN1,BETN2,H01,H02,HN1,HN2 
C
C  Evaluate the Jacobian:
C
C  ( 1.0   (-3.0*X(2)+10.0)*X(2)-2.0      34.0   )
C  ( 1.0   (3.0*X(2)+2.0)*X(2)-14.0       10.0   )
C
      REAL      FJAC(NVAR,NVAR)
      REAL      FPAR(*)
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      X(NVAR)

      DO I=1,NVAR
      DO J=1,NVAR
         FJAC(I,J)=0.
      ENDDO
      ENDDO
      
      FJAC(1,1)=-K01
      FJAC(1,3)=H01
      FJAC(1,7)= K01

      FJAC(2,1)= -ALPN1*X(3)
      FJAC(2,2)= -KN1
      FJAC(2,3)= +BETN1-ALPN1*X(1)
      FJAC(2,4)= -DELN1*X(6)
      FJAC(2,6)= -DELN1*X(4)

      FJAC(3,1)=  -HN1 +ALPN1*X(2)
      FJAC(3,2)=  -BETN1+ALPN1*X(1)
      FJAC(3,3)=  -KN1
      FJAC(3,4)=  +DELN1*X(5)
      FJAC(3,5)=  +DELN1*X(4)

      FJAC(4,2)=  +EPSN*X(6)
      FJAC(4,3)=  -EPSN*X(5)
      FJAC(4,4)=  -K02
      FJAC(4,5)=  -EPSN*X(3)
      FJAC(4,6)=  +EPSN*X(2)+H02

      FJAC(5,1)=  -ALPN2*X(6)
      FJAC(5,3)=  -DELTN2*X(4)
      FJAC(5,4)=  -DELTN2*X(3)
      FJAC(5,5)=  -KN2
      FJAC(5,6)=  +BETN2 -ALPN2*X(1)

      FJAC(6,1)= +ALPN2*X(5)
      FJAC(6,2)= +DELTN2*X(4)
      FJAC(6,4)= +DELTN2*X(2)-HN2
      FJAC(6,5)= -BETN2 +ALPN2*X(1)
      FJAC(6,6)= -KN2

      RETURN
      END

C
C
C  Here is the output from the run of this sample problem
C
C
C  PCPRB1:
C
C  PITCON sample program.
C  Freudenstein-Roth function
C  2 equations, 3 variables.
C
C  Step  Type of Point     X(1)         X(2)         X(3)
C
C    0  Start point      15.0000      -2.00000      0.000000E+00
C
C  PITCON 6.1
C  University of Pittsburgh Continuation Code
C
C  Last modification on 25 February 1991
C
C    1  Corrected        15.0000      -2.00000      0.000000E+00
C    2  Continuation     14.7105      -1.94205      0.653814E-01
C    3  Continuation     14.2846      -1.72915      0.268742
C    4  Limit point      14.2831      -1.74138      0.258578
C    5  Continuation     16.9061      -1.20941      0.546845
C    6  Continuation     24.9179     -0.599064      0.555136
C    7  Continuation     44.8783      0.487670      0.595261E-01
C    8  Continuation     60.0889       1.57585     -0.542365
C    9  Continuation    -11.1752       4.23516       1.55667
C   10  Target point     5.00000       4.00000       1.00000
C
C  We have reached the point we wanted.
C  The code may stop now.
C  FORTRAN STOP
C
      SUBROUTINE PITCON(DF,FPAR,FX,IERROR,IPAR,IWORK,LIW,NVAR,RWORK,
     *LRW,XR,SLNAME)

      REAL      EIGHT
      REAL      HALF
      REAL      HUNDRD
      REAL      ONE
      REAL      THREE
      REAL      TWO
      REAL      ZERO
C
      PARAMETER (EIGHT=8.0)
      PARAMETER (HALF=0.5)
      PARAMETER (HUNDRD=100.0)
      PARAMETER (ONE=1.0)
      PARAMETER (THREE=3.0)
      PARAMETER (TWO=2.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  COQUAL
      EXTERNAL  CORECT
      EXTERNAL  DF
      EXTERNAL  FX
      EXTERNAL  ISAMAX
      EXTERNAL  REPS
      EXTERNAL  ROOT
      EXTERNAL  SDOT
      EXTERNAL  SETSTP
      EXTERNAL  SLNAME
      EXTERNAL  SNRM2
      EXTERNAL  SSCAL
      EXTERNAL  START
      EXTERNAL  TANGNT
C
      INTRINSIC ABS
      INTRINSIC ATAN2
      INTRINSIC MAX
      INTRINSIC SIGN
      INTRINSIC SQRT
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      REAL      A
      REAL      ATCIPC
      REAL      ATCJPC
      REAL      B
      REAL      DETS
      REAL      DIRLPC
      REAL      FA
      REAL      FB
      REAL      FPAR(*)
      REAL      HTAN
      INTEGER   I
      INTEGER   ICALL
      INTEGER   ICRIT
      INTEGER   IERROR
      INTEGER   IFLAG
      INTEGER   IHOLD
      INTEGER   IMITL
      INTEGER   IPAR(*)
      INTEGER   IPC
      INTEGER   IPL
      INTEGER   ISAMAX
      INTEGER   IT
      INTEGER   ITMP
      INTEGER   IWORK(LIW)
      INTEGER   IWRITE
      INTEGER   JOB
      INTEGER   JPC
      INTEGER   LIM
      INTEGER   LOUNIT
      INTEGER   LPC
      INTEGER   LTC
      INTEGER   LTC1
      INTEGER   LWK
      INTEGER   LWK1
      INTEGER   LXC
      INTEGER   LXC1
      INTEGER   LXF
      INTEGER   LXF1
      INTEGER   MLSTEP
      INTEGER   MODSAV
      REAL      REPS
      REAL      RWORK(LRW)
      REAL      SCIPL
      REAL      SDOT
      REAL      SKALE
      REAL      SN
      REAL      SNL
      REAL      SNRM2
      REAL      STEPX
      REAL      TCIPC
      REAL      TCIPL
      REAL      TCOS
      REAL      TEMP
      REAL      TLIPC
      REAL      TMAX
      REAL      TMAX2
      REAL      TMP
      REAL      TMP1
      REAL      TSIN
      REAL      TSN
      REAL      XABS
      REAL      XDIF
      REAL      XLOW
      REAL      XR(NVAR)
      REAL      XUP
C
      SAVE      ICALL
C
      DATA      ICALL /0/
C
      ICALL=ICALL+1
C
C**********************************************************************
C
C  1.  Preparations.
C
C  Check entries of IWORK and RWORK, compute some constants.
C  For negative values of IWORK(1), compare user and finite difference 
C    jacobian and return.
C  For IWORK(1)=0, set up starting data, check starting point and return.
C  For positive IWORK(1), prepare to compute next point.
C
C**********************************************************************
C
      IERROR=0
      IF(IWORK(8).LT.1.OR.IWORK(8).GT.99)IWORK(8)=6
      LOUNIT=IWORK(8)
      IF(IWORK(7).LT.0.OR.IWORK(7).GT.3)IWORK(7)=1
      IWRITE=IWORK(7)
C
C  Initialization for first call
C
      IF(IWORK(1).EQ.0.OR.ICALL.EQ.1)THEN
        IF(IWRITE.GE.1)THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'PITCON 6.1'
          WRITE(LOUNIT,*)'University of Pittsburgh Continuation Code'
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'Last modification on 27 February 1991'
          WRITE(LOUNIT,*)' '
          ENDIF
        DO 10 I=11,17
          RWORK(I)=ZERO
10        CONTINUE
        DO 20 I=21,29
          RWORK(I)=ZERO
20        CONTINUE
        IWORK(10)=0
        IWORK(11)=0
        IWORK(12)=0
        DO 30 I=18,28
          IWORK(I)=0
30        CONTINUE
C
C  This is the only place REPS is called.  It might be replaced by
C  a call to the SLATEC machine parameter function R1MACH:
C
C    RWORK(8)=R1MACH(4)
C
C  or by simply setting RWORK(8) to the "machine epsilon", usually
C  taken as the smallest power of 2, EPMACH, such that (1.0+EPMACH) is
C  greater than 1.0, but (1.0+(EPMACH/2.0)) equals 1.0.
C
        RWORK(8)=REPS()
        IF(IWRITE.GE.2)WRITE(LOUNIT,919)RWORK(8)
        TCOS=SQRT(ONE-RWORK(8))
        TSIN=SQRT(RWORK(8))
        RWORK(10)=TWO*ATAN2(TSIN,TCOS)
        IF(NVAR.LE.1)THEN
          IERROR=1
          IF(IWRITE.GE.1)WRITE(LOUNIT,*)
     *    'PITCON - Number of variables must be at least 2.'
          GO TO 900
          ENDIF
        ENDIF
C
      LXC=29
      LXC1=LXC+1
      LXF=LXC+NVAR
      LXF1=LXF+1
      LTC=LXF+NVAR
      LTC1=LTC+1
      LWK=LTC+NVAR
      LWK1=LWK+1
      IWORK(13)=30
      IWORK(14)=LIW
      IWORK(15)=LWK+NVAR+1
      IWORK(16)=LRW
      IF(LIW.LT.IWORK(13))THEN
        WRITE(LOUNIT,*)'PITCON - Input value of LIW=',LIW
        WRITE(LOUNIT,*)'         Minimal acceptable value=',IWORK(13)
        IERROR=1
        GO TO 900
      ELSEIF(LRW.LT.IWORK(15))THEN
        WRITE(LOUNIT,*)'PITCON - Input value of LRW=',LRW
        WRITE(LOUNIT,*)'         Minimal acceptable value=',IWORK(15)
        IERROR=1
        GO TO 900
        ENDIF
C
C  Check entries of IWORK
C
      IF(IWORK(2).LT.1.OR.IWORK(2).GT.NVAR)IWORK(2)=NVAR
      IPC=IWORK(2)
      IF(IWORK(3).NE.1)IWORK(3)=0
      IF(IWORK(4).LT.0.OR.IWORK(4).GT.2)IWORK(4)=0
      IF(IWORK(5).LT.1.OR.IWORK(5).GT.NVAR)IWORK(5)=0
      IT=IWORK(5)
      IF(IWORK(6).LT.1.OR.IWORK(6).GT.NVAR)IWORK(6)=0
      LIM=IWORK(6)
      IF(IWORK(9).LT.0.OR.IWORK(9).GT.2)IWORK(9)=0
      IF(IWORK(17).LT.1)IWORK(17)=10
C
C  Check entries of RWORK
C
      IF(RWORK(1).LE.ZERO)RWORK(1)=SQRT(RWORK(8))
      IF(RWORK(2).LE.ZERO)RWORK(2)=SQRT(RWORK(8))
      IF(RWORK(3).LT.SQRT(RWORK(8)))RWORK(3)=SQRT(RWORK(8))
      IF(RWORK(4).LT.RWORK(3))THEN
        TEMP=NVAR
        RWORK(4)=MAX(RWORK(3),SQRT(TEMP))
        ENDIF
      IF(RWORK(5).LT.ZERO)THEN
        RWORK(5)=-RWORK(5)
        RWORK(6)=-RWORK(6)
        ENDIF
      IF(RWORK(5).LT.RWORK(3).OR.RWORK(5).GT.RWORK(4))THEN
        RWORK(5)=HALF*(RWORK(3)+RWORK(4))
        ENDIF
      HTAN=RWORK(5)
      IF(RWORK(6).NE.(-ONE))RWORK(6)=ONE
      IF(RWORK(20).LT.ONE.OR.RWORK(20).GT.HUNDRD)RWORK(20)=THREE
C
C  Check user Jacobian routine versus finite difference calculation.
C
      IF(IWORK(1).LT.0)THEN
        JOB=3
        CALL SLNAME(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,LIW,
     1  JOB,NVAR,RWORK,LRW,XR,RWORK(LWK1))
        GO TO 900
        ENDIF
C
C  Starting point check
C
      IF(IWORK(1).EQ.0)THEN
        CALL START(DF,FPAR,FX,IERROR,IPAR,IPC,IWORK,IWRITE,LIW,
     *  LOUNIT,LRW,NVAR,RWORK,RWORK(LTC1),RWORK(LWK1),RWORK(LXC1),
     *  RWORK(LXF1),XR,SLNAME)
        GO TO 900
        ENDIF
C
C***********************************************************************
C
C  2.  Target point 
C
C  If (IT.NE.0) target points are sought.  Check to see if target component
C  IT has value XIT lying between XC(IT) and XF(IT).  If so, get linearly
C  interpolated starting point, and use Newton's method to get target point.
C
C***********************************************************************
C
      IF(IT.EQ.0.OR.IWORK(10).LE.1)GO TO 300
      IF( (IWORK(1).EQ.3) .AND. (RWORK(7).EQ.RWORK(28)) .AND. 
     *    (IT.EQ.IWORK(11)) ) GO TO 300
      MODSAV=IWORK(4)
210   CONTINUE
      XLOW=RWORK(LXC+IT)
      XUP=RWORK(LXF+IT)
      IF(RWORK(7).LT.XLOW.AND.RWORK(7).LT.XUP)GO TO 300
      IF(RWORK(7).GT.XLOW.AND.RWORK(7).GT.XUP)GO TO 300
      IF(IWRITE.GE.2)WRITE(LOUNIT,*)
     *'PITCON - Attempt correction of approximate target point.'
C
C  Approximate the target point using the bracketing solutions.
C
      IF(XLOW.NE.XUP)THEN
        SKALE=(RWORK(7)-XLOW)/(XUP-XLOW)
      ELSE
        SKALE=ONE
        ENDIF
      CALL SCOPY(NVAR,RWORK(LXF1),1,XR,1)
      CALL SSCAL(NVAR,SKALE,XR,1)
      SKALE=ONE-SKALE
      CALL SAXPY(NVAR,SKALE,RWORK(LXC1),1,XR,1)
      XR(IT)=RWORK(7)
C
C  Call CORECT to compute the exact target point, holding index IT fixed.
C
      ICRIT=0
      CALL CORECT(DF,FPAR,FX,IERROR,IT,IPAR,IWORK,
     1 NVAR,RWORK,TMP,RWORK(LWK1),XR,LRW,LIW,ICRIT,SLNAME)
      IWORK(24)=IWORK(24)+IWORK(28)
      IF(IERROR.NE.0.AND.IWORK(4).GT.0)THEN
        IERROR=0
        IWORK(4)=IWORK(4)-1
        IF(IWRITE.GE.1)WRITE(LOUNIT,1080)IWORK(4)
        GO TO 210
        ENDIF
      IWORK(4)=MODSAV
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,*)'PITCON - Target point calculation failed.'
        GO TO 900
        ENDIF
C
C  Record the values of IT and XIT, and compute the arclength to the target 
C  point.
C
      IWORK(1)=3
      IWORK(11)=IT
      RWORK(28)=RWORK(7)
      SKALE=-ONE
      CALL SCOPY(NVAR,XR,1,RWORK(LWK1),1)
      CALL SAXPY(NVAR,SKALE,RWORK(LXC1),1,RWORK(LWK1),1)
      RWORK(14)=RWORK(12)+SNRM2(NVAR,RWORK(LWK1),1)
      IWORK(27)=IWORK(27)+1
      GO TO 900
C
C***********************************************************************
C
C  3.  Tangent and local continuation parameter calculation.  
C
C  Unless the tangent and limit point calculations were already performed,
C  (because the loop was interrupted for a limit point calculation),
C  set up and solve the equation for the tangent vector.
C
C  Force the tangent vector to be of unit length, and try to preserve
C  the "sense" or "direction" of the curve by forcing the IPL-th component
C  of the tangent vector to agree in sign with the IPL-th component of the
C  previous secant vector.  (On the first step, however, we have to use
C  the user's input direction to choose the sign).
C
C  Set the local continuation parameter IPC as the index of the
C  component of greatest magnitude in the tangent vector, unless a
C  limit point appears to be coming in that direction, and another
C  choice is available.
C
C***********************************************************************
C
  300 CONTINUE
      IF(IWORK(10).EQ.3) GO TO 600
C
C  STORE OLD TANGENT IN WORK ARRAY RWORK(LWK1)-RWORK(LWK+NVAR), COMPUTE
C  NEW TANGENT FOR CURRENT POINT XF IN RWORK(LXF1)-RWORK(LXF+NVAR).
C  THE TANGENT IS CALLED TC AND STORED IN RWORK(LTC1)-RWORK(LTC+NVAR).
C  NOTE THAT FOR A NORMAL RETURN THE CURRENT POINT XF WILL BE STORED
C  AS XC IN RWORK(LXC1)-RWORK(LXC+NVAR) AND HENCE THAT TC WILL THEN
C  CORRESPOND TO THE POINT XC
C
      IPL=IPC
      CALL SCOPY(NVAR,RWORK(LTC1),1,RWORK(LWK1),1)
      RWORK(29)=RWORK(17)
      CALL TANGNT(RWORK(17),FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,NVAR,
     1 RWORK,RWORK(LTC1),RWORK(LXF1),LIW,LRW,SLNAME)
      IF(IERROR.NE.0)THEN
        IF(IWRITE.GE.1)WRITE(LOUNIT,*)
     *  'PITCON - Tangent calculation failed.'
        GO TO 900
        ENDIF
C
C  FIND LARGEST AND SECOND LARGEST COMPONENTS OF TANGENT
C
      IPC=ISAMAX(NVAR,RWORK(LTC1),1)
      TMAX=RWORK(LTC+IPC)
      RWORK(LTC+IPC)=ZERO
      JPC=ISAMAX(NVAR,RWORK(LTC1),1)
      RWORK(LTC+IPC)=TMAX
      TMAX=ABS(TMAX)
      TMAX2=ABS(RWORK(LTC+JPC))
      IF(JPC.EQ.0)JPC=IPC
      IF(IWORK(3).EQ.1)JPC=IPC
C
C  ADJUST SIGN OF THE TANGENT VECTOR.
C  COMPARE THE SIGN OF THE COMPONENT TC(IPL) WITH THE
C  SIGN OF XF(IPL)-XC(IPL). IF A TARGET OR LIMIT POINT
C  HAS BEEN COMPUTED BETWEEN XC AND XF THEN USE XF(IPL)-XR(IPL)
C  INSTEAD. AT THE FIRST STEP COMPARE WITH THE USER INPUT DIRECTION.
C  IF THESE SIGNS DIFFER, CHANGE THE SIGN OF THE TANGENT VECTOR
C  AND THE SIGN OF THE DETERMINANT.
C
      SCIPL=RWORK(6)
      TCIPL=RWORK(LTC+IPL)
      IF(IWORK(10).GT.1)THEN
        SCIPL=RWORK(LXF+IPL)
        TMP=RWORK(LXC+IPL)
        IF(IWORK(1).EQ.3.OR.IWORK(1).EQ.4)TMP=XR(IPL)
        SCIPL=SCIPL-TMP      
        ENDIF
      IF(SIGN(ONE,TCIPL).NE.SIGN(ONE,SCIPL))THEN
        SKALE=-ONE
        CALL SSCAL(NVAR,SKALE,RWORK(LTC1),1)
        RWORK(17)=-RWORK(17)
        ENDIF
C
C  Unless we are computing a starting point, record the new state.
C
      IF(IWORK(10).GT.1)IWORK(10)=3
C
C  THE LOCATION OF THE LARGEST COMPONENT OF THE TANGENT VECTOR
C  WILL BE USED FOR THE LOCAL PARAMETERIZATION OF THE CURVE UNLESS
C  A LIMIT POINT IN IPC APPEARS TO BE COMING.
C  TO CHECK THIS, WE COMPARE TCIPC=TC(IPC) AND THE
C  SECOND LARGEST COMPONENT TCJPC=TC(JPC).  IF TCJPC IS NO LESS
C  THAN .1 OF TCIPC, AND TC(JPC) IS LARGER THAN TL(JPC),
C  WHEREAS TC(IPC) IS LESS THAN TL(IPC), WE WILL RESET THE
C  LOCAL PARAMETER IPC TO JPC.
C  BUT IF NOT ALLOWED TO SWITCH PARAMETERS, IGNORE THE CHECK
C
      IF(IWORK(3).NE.1)THEN
        ATCIPC=ABS(RWORK(LTC+IPC))
        ATCJPC=ABS(RWORK(LTC+JPC))
        IF(JPC.NE.IPC.AND.IWORK(10).GT.1)THEN
          TLIPC=RWORK(LWK+IPC)
          TCIPC=RWORK(LTC+IPC)
          TMP=ABS(RWORK(LWK+JPC))
          IF((SIGN(ONE,TCIPC).EQ.SIGN(ONE,TLIPC))
     1      .AND.(ATCIPC.LT.ABS(TLIPC))
     2      .AND.(ATCJPC.GE.MAX(0.1*ATCIPC,TMP)))THEN
            IF(IWRITE.GE.3)WRITE(LOUNIT,350)IPC
            ITMP=IPC
            IPC=JPC
            JPC=ITMP
            ENDIF
          ENDIF
        ENDIF
      IF(IWRITE.GE.3.AND.IWORK(3).EQ.0)THEN
        WRITE(LOUNIT,360)IPC
        WRITE(LOUNIT,370)JPC
        ENDIF
C
C  RECORD VALUES OF THE COMPONENT TCIPC=TC(IPC) OF
C  NEW TANGENT VECTOR, AS WELL AS OF THE CORRESPONDING COMPONENT
C  TLIPC=WK(IPC) OF THE OLD TANGENT VECTOR.
C  SET SIGN OF THE DETERMINANT.
C  FOR USE IN THE LIMIT POINT COMPUTATION RECORD THE
C  COMPONENT TC(LIM) OF THE NEW TANGENT VECTOR
C
      IWORK(2)=IPC
      IWORK(12)=JPC
      RWORK(24)=RWORK(LTC+IPC)
      RWORK(25)=RWORK(LWK+IPC)
      RWORK(6)=RWORK(17)
      IF(LIM.GT.0)THEN
        RWORK(27)=RWORK(26)
        RWORK(26)=RWORK(LTC+LIM)
        IF(IWRITE.GE.2)WRITE(LOUNIT,1020)RWORK(26)
        ENDIF
C
C  COMPUTE ALPHLC, THE ANGLE BETWEEN OLD TANGENT AND TANGENT AT XC
C
      IF(IWORK(10).LE.1) GO TO 600
      TCOS=SDOT(NVAR,RWORK(LWK1),1,RWORK(LTC1),1)
      IF(TCOS.GT.ONE)TCOS=ONE
      IF(TCOS.LT.(-ONE))TCOS=-ONE
      TSIN=SQRT(ONE-TCOS*TCOS)
      RWORK(11)=ATAN2(TSIN,TCOS)
      IF((IWORK(10).LE.1).OR.(RWORK(17).EQ.RWORK(29)))GO TO 400
      IF(IWRITE.GE.2)WRITE(LOUNIT,*)
     *'PITCON - Possible bifurcation point detected.'
C
C***********************************************************************
C
C  4.  Limit point check.  
C
C  Skip this section if LIM=0.
C
C  Otherwise, user has requested a search for limit points in a given 
C  index by setting LIM to a nonzero value.
C
C  Compare LIM-th components of previous and current tangent vectors.
C  If a sign difference occurs, we assume a limit point has been passed.
C  Attempt to compute a point XR between the previous and current points,
C  for which the LIM-th component of the tangent is zero.
C
C  This search will be guided by a rootfinder.  The scalar function
C  to be zeroed out is the LIM-th tangent vector component.
C
C***********************************************************************
C
  400 CONTINUE
      IF((LIM.LE.0).OR.
     *   (IWORK(10).LE.1).OR.
     *   (SIGN(ONE,RWORK(26)).EQ.SIGN(ONE,RWORK(27)))) GO TO 600
      IF(IWRITE.GE.2)WRITE(LOUNIT,*)
     *'PITCON - Attempt correction of approximate limit point.'
      MODSAV=IWORK(4)
410   CONTINUE
C
C  SEE IF EITHER ENDPOINT CAN BE ACCEPTED AS LIMIT POINT.
C  IF XC IS LIMIT POINT, WORK ALREADY CONTAINS TANGENT AT XC.
C
      IF(ABS(RWORK(27)).LE.HALF*RWORK(1)) THEN
        CALL SCOPY(NVAR,RWORK(LXC1),1,XR,1)
        GO TO 490
        ENDIF
      IF(ABS(RWORK(26)).LE.HALF*RWORK(1)) THEN
        CALL SCOPY(NVAR,RWORK(LXF1),1,XR,1)
        CALL SCOPY(NVAR,RWORK(LTC1),1,RWORK(LWK1),1)
        GO TO 490
        ENDIF
C
C  If interval is extremely small, simply assign one
C  endpoint of the interval as the answer.
C
      XDIF=ABS(RWORK(LXC+LIM)-RWORK(LXF+LIM))
      XABS=MAX(ABS(RWORK(LXC+LIM)),ABS(RWORK(LXF+LIM)))
      IF(XDIF.LE.EIGHT*RWORK(8)*(ONE+XABS)) THEN
        IF(ABS(RWORK(27)).GT.ABS(RWORK(26))) THEN
          CALL SCOPY(NVAR,RWORK(LXF1),1,XR,1)
          CALL SCOPY(NVAR,RWORK(LTC1),1,RWORK(LWK1),1)
        ELSE
          CALL SCOPY(NVAR,RWORK(LXC1),1,XR,1)
          ENDIF
        GO TO 490
        ENDIF
C
C  BEGIN ROOT-FINDING ITERATION WITH INTERVAL (0,1) AND
C  FUNCTION VALUES TLLIM, TCLIM.
C
      A=ZERO
      FA=RWORK(27)
      B=ONE
      FB=RWORK(26)
C
C  FIND LPC, INDEX OF MAXIMUM ENTRY OF SECANT (EXCEPT THAT
C  LPC MUST NOT EQUAL LIM) AND SAVE SIGN IN DIRLPC
C  SO THAT NEW TANGENTS MAY BE PROPERLY SIGNED.
C
      TEMP=RWORK(LXC+LIM)
      RWORK(LXC+LIM)=RWORK(LXF+LIM)
      CALL SCOPY(NVAR,RWORK(LXF1),1,XR,1)
      SKALE=-ONE
      CALL SAXPY(NVAR,SKALE,RWORK(LXC1),1,XR,1)
      LPC=ISAMAX(NVAR,XR,1)
      RWORK(LXC+LIM)=TEMP
      DIRLPC=SIGN(ONE,RWORK(LXF+LPC)-RWORK(LXC+LPC))
C
C  SET FIRST APPROXIMATION TO ROOT TO WHICHEVER ENDPOINT
C  HAS SMALLEST LIM-TH COMPONENT OF TANGENT
C
      IF(ABS(RWORK(26)).GE.ABS(RWORK(27))) THEN
        SN=ZERO
        TSN=RWORK(27)
        CALL SCOPY(NVAR,RWORK(LXC1),1,XR,1)
      ELSE
        SN=ONE
        TSN=RWORK(26)
        CALL SCOPY(NVAR,RWORK(LXF1),1,XR,1)
        ENDIF
C
C  CALL ROOTFINDER FOR APPROXIMATE ROOT SN.  USE LINEAR COMBINATION
C  OF PREVIOUS ROOT SNL, AND ONE OF 0.0 AND 1.0 TO GET A STARTING
C  POINT FOR CORRECTION.  RETURN TO CURVE ON LINE X(LPC)=CONSTANT,
C  COMPUTE TANGENT THERE, AND SET FUNCTION VALUE TO TANGENT(LIM)
C
      MLSTEP=25
      IMITL=0
      IF(IWRITE.GE.2)THEN
        TMP=ZERO
        WRITE(LOUNIT,1100)TMP,RWORK(27)
        TMP=ONE
        WRITE(LOUNIT,1100)TMP,RWORK(26)
        ENDIF
441   CONTINUE
      SNL=SN
      CALL ROOT(A,FA,B,FB,SN,TSN,IMITL,IFLAG,IERROR,RWORK(8))
      IWORK(23)=IWORK(23)+1
      IF(IERROR.NE.0)THEN
        IF(IWRITE.GE.1)WRITE(LOUNIT,*)
     *  'PITCON - Bad interval for rootfinder'
        GO TO 900
        ENDIF
      IF(IFLAG.EQ.(-1).OR.IFLAG.EQ.0)GO TO 490
C
C  FIND WHETHER SN LIES IN (0.0,SNL) OR (SNL,1.0).
C  USE APPROPRIATE LINEAR COMBINATION TO GET STARTING POINT.
C
      IF(SN.LE.SNL) THEN
C
C  SET X(SN)=(SNL-SN)/(SNL-0.0)*X(0.0)+(SN-0.0)/(SNL-0.0)*X(SNL)
C
        IF(SNL.LE.ZERO)THEN
          SKALE=ZERO
        ELSE
          SKALE=SN/SNL
          ENDIF
        IF(SKALE.LE.ZERO)SKALE=ZERO
        IF(SKALE.GE.ONE)SKALE=ONE
        CALL SSCAL(NVAR,SKALE,XR,1)
        SKALE=ONE-SKALE
        CALL SAXPY(NVAR,SKALE,RWORK(LXC1),1,XR,1)
      ELSE
C
C  SET X(SN)=(SN-SNL)/(1.0-SNL)*X(1.0)+(1.0-SN)/(1.0-SNL)*X(SNL)
C
        IF(SNL.GE.ONE)THEN
          SKALE=ZERO
        ELSE
          SKALE=(ONE-SN)/(ONE-SNL)
          ENDIF
        IF(SKALE.LE.ZERO)SKALE=ZERO
        IF(SKALE.GE.ONE)SKALE=ONE
        CALL SSCAL(NVAR,SKALE,XR,1)
        SKALE=ONE-SKALE
        CALL SAXPY(NVAR,SKALE,RWORK(LXF1),1,XR,1)
        ENDIF
C
582   CONTINUE
      ICRIT=0
      CALL CORECT(DF,FPAR,FX,IERROR,LPC,IPAR,IWORK,
     1 NVAR,RWORK,TMP,RWORK(LWK1),XR,LRW,LIW,ICRIT,SLNAME)
      IF(IERROR.NE.0.AND.IWORK(4).GT.0)THEN
        IWORK(4)=IWORK(4)-1
        CALL SCOPY(NVAR,RWORK(LXC1),1,XR,1)
        SKALE=ONE-SN
        CALL SSCAL(NVAR,SKALE,XR,1)
        SKALE=SN
        CALL SAXPY(NVAR,SKALE,RWORK(LXF1),1,XR,1)
        IF(IWRITE.GE.1)WRITE(LOUNIT,1090)IWORK(4)
        GO TO 582
        ENDIF
      IWORK(4)=MODSAV 
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,*)
     *  'PITCON - Corrector failed on limit point.'
        GO TO 900
        ENDIF
      CALL TANGNT(TEMP,FX,DF,FPAR,IERROR,LPC,IPAR,IWORK,NVAR,
     1 RWORK,RWORK(LWK1),XR,LIW,LRW,SLNAME)
      IF(IERROR.NE.0)THEN
        IF(IWRITE.GE.1)WRITE(LOUNIT,*)
     *  'PITCON - Tangent error in limit point calculation.'
        GO TO 900
        ENDIF
C
C  Adjust the sign of the tangent vector so the LPC-th component
C  has the same sign as the LPC-th component of the secant.
C
      IF(DIRLPC.NE.SIGN(ONE,RWORK(LWK+LPC))) THEN
        SKALE=-ONE
        CALL SSCAL(NVAR,SKALE,RWORK(LWK1),1)
        ENDIF
C
C  SEE IF WE CAN ACCEPT THE NEW POINT BECAUSE TANGENT(LIM) IS SMALL
C  OR MUST STOP BECAUSE 25 ITERATIONS TAKEN,
C  OR IF WE MUST GO ON.
C
      TSN=RWORK(LWK+LIM)
      IF(IWRITE.GE.2)WRITE(LOUNIT,1100)SN,TSN
      IF(ABS(TSN).LE.RWORK(1))GO TO 490
      IF(IMITL.LT.MLSTEP)GO TO 441
C
C  Limit point iteration has failed.  We set an error flag,
C  but we return the partial results of the abortive computation
C  anyway.
C
      IERROR=8
      IF(IWRITE.GE.1)WRITE(LOUNIT,*)
     *'PITCON - Too many steps in limit point calculation.'
C
C  Limit point iteration is over.  
C  Compute and store information.
C
  490 CONTINUE
      CALL SCOPY(NVAR,XR,1,RWORK(LWK1),1)
      SKALE=-ONE
      CALL SAXPY(NVAR,SKALE,RWORK(LXC1),1,RWORK(LWK1),1)
      RWORK(14)=RWORK(12)+SNRM2(NVAR,RWORK(LWK1),1)
      IWORK(27)=IWORK(27)+1
      IWORK(1)=4
      GO TO 900
C
C***********************************************************************
C
C  5.  Compute next predictor step length, HTAN.
C
C***********************************************************************
C
  600 CONTINUE
      IF(IWORK(10).GT.1)CALL SETSTP(HTAN,IWORK,LIW,LRW,RWORK)
C
C***********************************************************************
C
C  6.  Continuation step
C
C  Our current data is the current point XC, its tangent vector TC, and
C  a steplength HTAN.  We predict the location of the next point on the
C  curve using the Euler approximation XR=XC+HTAN*TC.
C
C  Newton iteration is applied to this point, to force it to lie on the
C  curve.  In order to make the system square, an augmenting equation
C  is added to the system, specifying that XR(IPC)=XC(IPC)+HTAN*TC(IPC).
C  (The right hand side is a constant.)
C
C  If the Newton correction process fails, the stepsize is reduced and
C  prediction and correction retried.  Failure will most likely be
C  signaled by repeated step reductions, until the minimum allowable
C  stepsize is reached.  If this occurs, PITCON has failed, and cannot
C  proceed along the curve any more.
C
C***********************************************************************
C
  700 CONTINUE
      IWORK(18)=0
      IPC=IWORK(2)
      JPC=IWORK(12)
      IHOLD=IPC
      MODSAV=IWORK(4)
710   CONTINUE
      IF(IWRITE.GE.2)WRITE(LOUNIT,715)HTAN
      CALL SCOPY(NVAR,RWORK(LXF1),1,XR,1)
      CALL SAXPY(NVAR,HTAN,RWORK(LTC1),1,XR,1)
  740 CONTINUE
      IF(IWRITE.GE.2.AND.IWORK(3).EQ.0)WRITE(LOUNIT,725)IHOLD
      ICRIT=0
      CALL CORECT(DF,FPAR,FX,IERROR,IHOLD,IPAR,IWORK,
     1 NVAR,RWORK,STEPX,RWORK(LWK1),XR,LRW,LIW,ICRIT,SLNAME)
      IWORK(25)=IWORK(25)+IWORK(28)
      IF(IERROR.EQ.0)GO TO 800
C
C  Only VERY fatal errors should abort the correction process.
C  Abort this process only after kicking and screaming.
C
      IF(IERROR.EQ.2)THEN
        WRITE(LOUNIT,*)'PITCON - Fatal error during Newton correction'
        WRITE(LOUNIT,*)'         of a continutation point!'
        GO TO 900
        ENDIF
C
C  FOR CASES WHERE CORRECTOR FAILED TO CONVERGE,
C  AND THE VALUE OF IHOLD WAS NOT JPC,
C  AND THE USER HAS NOT FIXED THE CONTINUATION PARAMETER
C  KEEP THE SAME STEPSIZE BUT USE JPC INSTEAD OF IPC.
C
      IF(IWORK(3).NE.1.AND.IHOLD.NE.JPC.AND.JPC.NE.0)THEN
        IHOLD=JPC
        IF(IERROR.EQ.5)GO TO 740
        GO TO 710
        ENDIF
C
C  IF JPC FAILS AS IHOLD VALUE, RETURN TO IPC.
C
      IHOLD=IPC
C
C  IF WE ARE USING SOME FORM OF MODIFIED NEWTON, BUT
C  WE HAVE REACHED MINIMUM STEPSIZE OR 
C  HAVE HAD TWO FAILURES IN A ROW ON THIS STEP,
C  USE A BETTER NEWTON METHOD
C
      IF(IWORK(4).GT.0.AND.
     *  (HTAN.LT.RWORK(20)*RWORK(3).OR.IWORK(18).GE.2) )THEN
        IWORK(4)=IWORK(4)-1
        WRITE(LOUNIT,1130)IWORK(4)
        IF(IERROR.NE.5)GO TO 710
        GO TO 740
        ENDIF
C
C  NO CONVERGENCE, TRY A SMALLER STEPSIZE
C
      IF(HTAN.LT.RWORK(20)*RWORK(3))THEN
        IERROR=4
        IF(IWRITE.GE.1)WRITE(LOUNIT,*) 
     *  'PITCON - Stepsize fell below minimum.'
        GO TO 900
        ENDIF
      HTAN=HTAN/RWORK(20)
      IWORK(18)=IWORK(18)+1
      IF(IERROR.NE.5)GO TO 710
      SKALE=ONE/RWORK(20)
      CALL SSCAL(NVAR,SKALE,XR,1)
      SKALE=ONE-SKALE
      CALL SAXPY(NVAR,SKALE,RWORK(LXF1),1,XR,1)
      GO TO 740
C
C***********************************************************************
C
C  7.  Successful step.  Update information.
C
C***********************************************************************
C
  800 CONTINUE
      IWORK(1)=2
      IF(IWORK(3).NE.1)IWORK(2)=IHOLD
      IWORK(4)=MODSAV
      IWORK(10)=2
      IWORK(26)=IWORK(26)+IWORK(18)
      IWORK(27)=IWORK(27)+1
      RWORK(5)=HTAN
      RWORK(22)=RWORK(21)
      IF(IWRITE.GE.1.AND.IWORK(18).GT.0)WRITE(LOUNIT,755)IWORK(18)
      CALL SCOPY(NVAR,XR,1,RWORK(LWK1),1)
      SKALE=-ONE
      CALL SAXPY(NVAR,SKALE,RWORK(LXF1),1,RWORK(LWK1),1)
      RWORK(21)=SNRM2(NVAR,RWORK(LWK1),1)
      RWORK(12)=RWORK(13)
      RWORK(13)=RWORK(12)+RWORK(21)
      RWORK(14)=RWORK(13)
C
C  UPDATE VALUES OF COMPUTED POINTS.  STORE OLD POINT IN
C  RWORK(LXC1)-RWORK(LXC+NVAR) AND NEW POINT IN RWORK(LXF1)-
C  RWORK(LXF+NVAR). COMPUTE AND STORE CORDIS, THE MAXIMUM NORM
C  OF THE CORRECTOR DISTANCE.
C
      RWORK(15)=ZERO
      IF(IWORK(28).NE.0)THEN
        DO 810 I=1,NVAR
          TMP1=XR(I)-(RWORK(LXF+I)+HTAN*RWORK(LTC+I))
          IF(ABS(TMP1).GT.RWORK(15))RWORK(15)=ABS(TMP1)
 810      CONTINUE
        ENDIF
      CALL SCOPY(NVAR,RWORK(LXF1),1,RWORK(LXC1),1)
      CALL SCOPY(NVAR,XR,1,RWORK(LXF1),1)
C
C  COMPUTE CONVERGENCE QUALITY
C
      CALL COQUAL(STEPX,IWORK,LIW,RWORK,LRW)
C
C**********************************************************************
C
C  8.  All returns to calling program should exit here.
C
C**********************************************************************
C
  900 CONTINUE
      IF(IERROR.EQ.0)THEN
      ELSEIF(IERROR.EQ.1)THEN
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)'         Storage or data error.'
      ELSEIF(IERROR.EQ.2)THEN
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)
     *  '         User-defined function or jacobian error.'
      ELSEIF(IERROR.EQ.3)THEN
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)'         Solver failed.'
      ELSEIF(IERROR.EQ.4)THEN
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)'         Iteration failed.'
      ELSEIF(IERROR.EQ.5)THEN
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)'         Too many corrector steps.'
      ELSEIF(IERROR.EQ.6)THEN
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)'         Zero tangent vector.'
      ELSEIF(IERROR.EQ.7)THEN
        WRITE(LOUNIT,*)'PITCON - Nonfatal error.'
        WRITE(LOUNIT,*)'         Root finder failed on limit point.'
      ELSEIF(IERROR.EQ.8)THEN
        WRITE(LOUNIT,*)'PITCON - Nonfatal error.'
        WRITE(LOUNIT,*)'         Too many limit point steps.'
      ELSE
        IERROR=9
        WRITE(LOUNIT,*)'PITCON - FATAL ERROR.'
        WRITE(LOUNIT,*)'         Unidentified error.'
        ENDIF
      RETURN
350   FORMAT(' PITCON - Expect limit point in index ',I5)
360   FORMAT(' PITCON - Continuation index: First choice',I5)
370   FORMAT(' PITCON -                    Second choice',I5)
715   FORMAT(' PITCON - Predictor using stepsize ',G14.6)
725   FORMAT(' PITCON - Corrector fixing index',I5)
755   FORMAT(' PITCON - Predictor stepsize was reduced ',I6,' times.')
919   FORMAT(' PITCON - Machine epsilon=',G14.6)
1020  FORMAT(' PITCON - Tangent vector has limit component =',G14.6)
1080  FORMAT(' PITCON - Retry target computation with IWORK(4)=',I6)
1090  FORMAT(' PITCON - Retry limit computation with IWORK(4)=',I6)
1100  FORMAT(' PITCON - For S=',G14.6,' TAN(X(S))(LIM)=',G14.6)
1130  FORMAT(' PITCON - Retrying step with IWORK(4)=',I6)
      END
      SUBROUTINE COQUAL(STEPX,IWORK,LIW,RWORK,LRW)
C
C***********************************************************************
C
C  COQUAL computes the factor QUAL which is based on the 'quality'
C  (number of steps, last step/total step) of the corrector iteration.  
C  This factor is used as part of the stepsize algorithm.  
C
C  See the paper "On Steplength Algorithms for a Class of Continuation
C  Methods", listed in the documentation.
C
C  STEPX  = Input, size of the last correction.
C  CORDIS = Input, distance between first and last point.
C  MAXCOR = Input, maximum number of corrections allowed.
C  MODNEW = Input, defines type of corrector process used.
C  NCOR   = Input, actual number of corrections used.
C  QUAL   = Output, the quality factor, between 1/8 and 8,
C           with 1/8 poor, 1 average, 8 superior.  Used to
C           modify the next stepsize.
C
      REAL      EIGHT
      REAL      ONE
      REAL      THGIE
C
      PARAMETER (EIGHT=8.0)
      PARAMETER (ONE=1.0)
      PARAMETER (THGIE=0.125)
C
      INTEGER   LIW
      INTEGER   LRW
C
      REAL      BASE
      REAL      BOT
      REAL      CORDIS
      REAL      EPSATE
      REAL      ESAB
      REAL      EXPO
      INTEGER   IBOT
      INTEGER   ITOP
      INTEGER   IWORK(LIW)
      INTEGER   LOUNIT
      INTEGER   MAXCOR
      INTEGER   MODNEW
      INTEGER   NAVE
      INTEGER   NMAX
      REAL      QUAL
      REAL      RWORK(LRW)
      REAL      STEPX
      REAL      TERM
      REAL      TEST
      REAL      TOP
C
      CORDIS=RWORK(15)
      EPSATE=EIGHT*RWORK(8)
      MODNEW=IWORK(4)
      MAXCOR=IWORK(17)
C
C  SEE IF MINIMAL NUMBER OF STEPS TAKEN
C
      IF((IWORK(28).LE.1).OR.(CORDIS.LE.EPSATE)) THEN
        QUAL=EIGHT
        GO TO 88
        ENDIF
C
C  SEE IF AVERAGE NUMBER OF STEPS TAKEN
C
      IF(MODNEW.EQ.0)THEN
        NAVE=(MAXCOR-1)/2
      ELSE
        NAVE=MAXCOR
        ENDIF
      IF(IWORK(28).EQ.NAVE)THEN
        QUAL=ONE
        GO TO 88
        ENDIF
C
C  SEE IF MAXIMUM NUMBER OF STEPS TAKEN
C
      IF(MODNEW.EQ.0)THEN
        NMAX=MAXCOR
      ELSE
        NMAX=2*MAXCOR
        ENDIF
      IF(IWORK(28).GE.NMAX)THEN
        QUAL=THGIE
        GO TO 88
        ENDIF
C
C  COMPUTE QUAL
C
C  FOR MODNEW=0, LET W=(STEPX/CORDIS),
C                    IEXP=1/(2**(NCOR-1)-1)
C                    U=W**IEXP
C                    JEXP=2**(NCOR-NAVE)
C              THEN  QUAL=(U+1+(1/U)) / (U**JEXP+1+(1/U**JEXP))
C
C  OTHERWISE,    EXP=(NCOR-NAVE)/(NCOR-1)
C                W=(STEPX/CORDIS)
C                QUAL=W**EXP
C
      IF(MODNEW.EQ.0)THEN
        ITOP=2**(IWORK(28)-NAVE)
        IBOT=2**(IWORK(28)-1)-1
        EXPO=ONE/FLOAT(IBOT)
        BASE=(STEPX/CORDIS)**EXPO
        TOP=BASE+ONE+(ONE/BASE)
        TERM=BASE**ITOP
        BOT=TERM+ONE+(ONE/TERM)
        QUAL=TOP/BOT
        IF(QUAL.GT.EIGHT)QUAL=EIGHT
        IF(QUAL.LT.THGIE)QUAL=THGIE
        GO TO 88
        ENDIF
C
C  TRY TO ANTICIPATE UNDERFLOW AND OVERFLOW
C
      ITOP=IWORK(28)-NAVE
      IBOT=IWORK(28)-1
      EXPO=FLOAT(IBOT)/FLOAT(ITOP)
      TEST=EIGHT**EXPO
      BASE=STEPX/CORDIS
      ESAB=CORDIS/STEPX
      IF((IWORK(28).LT.NAVE.AND.TEST.GT.BASE)
     1  .OR.(IWORK(28).GT.NAVE.AND.TEST.LT.BASE))THEN
        QUAL=EIGHT
        GO TO 88
        ENDIF
      IF((IWORK(28).LT.NAVE.AND.TEST.GT.ESAB)
     1  .OR.(IWORK(28).GT.NAVE.AND.TEST.LT.ESAB))THEN
        QUAL=THGIE
        GO TO 88
        ENDIF
C
C  COMPUTE QUAL
C
      EXPO=FLOAT(ITOP)/FLOAT(IBOT)
      QUAL=BASE**EXPO
      IF(QUAL.GT.EIGHT)QUAL=EIGHT
      IF(QUAL.LT.THGIE)QUAL=THGIE
88    CONTINUE
      RWORK(23)=QUAL
      IF(IWORK(7).GE.3)THEN
        LOUNIT=IWORK(8)
        WRITE(LOUNIT,1000)QUAL
        ENDIF
      RETURN
1000  FORMAT(' COQUAL - Corrector convergence quality factor=',G14.6)
      END
      SUBROUTINE CORECT(DF,FPAR,FX,IERROR,IHOLD,IPAR,IWORK,
     1 NVAR,RWORK,STEPX,WK,XR,LRW,LIW,ICRIT,SLNAME)
C
C***********************************************************************
C
C  CORECT PERFORMS THE CORRECTION OF AN APPROXIMATE
C  SOLUTION OF THE EQUATION F(XR)=0.0. FOR THE CORRECTION  EITHER
C  NEWTON'S METHOD OR THE CHORD NEWTON METHOD IS USED. IN THE LATTER
C  CASE THE JACOBIAN IS  EVALUATED ONLY AT THE STARTING POINT.
C  IF B IS THE VALUE OF XR(IHOLD) FOR THE INPUT STARTING POINT,
C  THEN THE AUGMENTING EQUATION IS XR(IHOLD)=B, THAT IS, THE IHOLD-TH
C  COMPONENT OF XR IS TO BE HELD FIXED.  THE AUGMENTED SYSTEM TO BE
C  SOLVED IS THEN DFA(XR,IHOLD)*(-DELX)=FA(XR)
C
C  ACCEPTANCE AND REJECTION CRITERIA-
C
C  LET NCOR DENOTE THE CURRENT ITERATION INDEX; THAT IS,  NCOR=0 FOR
C  THE PREDICTED POINT WHICH SERVES AS STARTING POINT.
C  After each iteration, the current point will be accepted
C  if any of the following conditions hold:
C
C  STRONG ACCEPTANCE CRITERION
C
C  1.  FNRM.LE.ABSERR AND STEPX.LE.(ABSERR+RELERR*XNRM))
C
C  WEAK ACCEPTANCE CRITERIA
C
C  2.  (NCOR.EQ.0) AND
C      FNRM.LE.(.5*ABSERR)
C  3.  FNRM.LE.EPSATE OR STEPX.LE.EPSATE
C  4.  (NCOR.GE.2) AND
C      (FNRM+FNRML).LE.ABSERR AND
C      STEPX.LE.8*(ABSERR+RELERR*XNRM)
C  5.  (NCOR.GE.2) AND
C      FNRM.LE.8.0*ABSERR AND
C      (STEPX+STEPXL).LE.(ABSERR+RELERR*XNRM)
C
C  AFTER EACH STEP OF THE ITERATION, THE PROCESS IS TO BE ABORTED IF
C  ANY OF THE FOLLOWING CONDITIONS HOLDS
C
C  1.  FNRM.GT.(FMP*FNRML+ABSERR)
C  2.  (NCOR.GE.2) AND (ICRIT.EQ.0) AND
C      (STEPX.GT.(FMP*STEPXL+ABSERR))
C  3.  NCOR.GE.MAXCOR (NUMBER OF ITERATIONS EXCEEDED).
C
C  HERE FMP=2.0 FOR NCOR.EQ.1, FMP=1.05 FOR NCOR.GT.1
C
C  ERROR CONDITIONS
C         IERROR=1  DATA OR STORAGE ERROR
C         IERROR=2  ERROR IN FUNCTION OR DERIVATIVE CALL
C         IERROR=3  SOLVER FAILED
C         IERROR=4  UNSUCCESSFUL ITERATION
C         IERROR=5  TOO MANY CORRECTOR STEPS
C
      REAL      EIGHT
      REAL      HALF
      REAL      ONE
      REAL      ONEFIV
      REAL      TWO
      REAL      ZERO
C
      PARAMETER (EIGHT=8.0)
      PARAMETER (HALF=0.5)
      PARAMETER (ONE=1.0)
      PARAMETER (ONEFIV=1.05)
      PARAMETER (TWO=2.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  DF
      EXTERNAL  FX
      EXTERNAL  ISAMAX
      EXTERNAL  SAXPY
      EXTERNAL  SLNAME
      EXTERNAL  SNRM2
C
      INTRINSIC ABS
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      REAL      ABSERR
      REAL      DETS
      REAL      EPSATE
      REAL      FMP
      REAL      FNRM
      REAL      FNRML
      REAL      FPAR(*)
      INTEGER   I
      INTEGER   ICRIT
      INTEGER   IERROR
      INTEGER   IFMAX
      INTEGER   IHOLD
      INTEGER   IPAR(*)
      INTEGER   ISAMAX
      INTEGER   IWORK(LIW)
      INTEGER   IWRITE
      INTEGER   IXMAX
      INTEGER   JOB
      INTEGER   KSMAX
      INTEGER   LOUNIT
      INTEGER   MAXCOR
      INTEGER   MAXNEW
      INTEGER   MODNEW
      REAL      RELERR
      REAL      RWORK(LRW)
      REAL      SKALE
      REAL      SNRM2
      REAL      STEPX
      REAL      STEPXL
      REAL      TLSTEP
      REAL      WK(NVAR)
      REAL      XNRM
      REAL      XR(NVAR)
      REAL      XVALUE
C
C  Initialize.
C
      ABSERR=RWORK(1)
      RELERR=RWORK(2)
      EPSATE=EIGHT*RWORK(8)
      MODNEW=IWORK(4)
      IWRITE=IWORK(7)
      LOUNIT=IWORK(8)
      MAXCOR=IWORK(17)
      IERROR=0
      IWORK(28)=0
      MAXNEW=MAXCOR
      IF(MODNEW.NE.0)MAXNEW=2*MAXCOR
      FMP=TWO
      STEPX=ZERO
      XVALUE=XR(IHOLD)
C
C  STORE INITIAL FUNCTION VALUE IN THE WORK VECTOR WK
C
      CALL FX(NVAR,FPAR,IPAR,XR,WK,IERROR)
      IWORK(22)=IWORK(22)+1
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,*)
     *  'CORECT - Error flag from user function routine.'
        RETURN
        ENDIF
      WK(NVAR)=XR(IHOLD)-XVALUE
      IFMAX=ISAMAX(NVAR,WK,1)
      FNRM=SNRM2(NVAR-1,WK,1)
      IXMAX=ISAMAX(NVAR,XR,1)
      XNRM=SNRM2(NVAR,XR,1)
      IF(IWRITE.GE.2)WRITE(LOUNIT,70)IWORK(28),FNRM,IFMAX
      IF(FNRM.LE.HALF*ABSERR)RETURN
C
C  NEWTON ITERATION LOOP BEGINS
C
      DO 100 I=1,MAXNEW
        IWORK(28)=I
C
C  SOLVE SYSTEM FPRIME*(-DELX)=FX FOR NEWTON DIRECTION DELX
C
        JOB=0
        IF(I.NE.1.AND.I.NE.MAXCOR.AND.MODNEW.EQ.1)JOB=1
        IF(MODNEW.EQ.2)JOB=1
        CALL SLNAME(DETS,FX,DF,FPAR,IERROR,IHOLD,IPAR,IWORK,LIW,
     1  JOB,NVAR,RWORK,LRW,XR,WK)
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,12)IERROR
          RETURN
          ENDIF
C
C  UPDATE X, COMPUTE SIZE OF STEP AND OF X
C
        SKALE=-ONE
        CALL SAXPY(NVAR,SKALE,WK,1,XR,1)
        STEPXL=STEPX
        KSMAX=ISAMAX(NVAR,WK,1)
        STEPX=ABS(WK(KSMAX))
        IXMAX=ISAMAX(NVAR,XR,1)
        XNRM=SNRM2(NVAR,XR,1)
C
C  Compute F(X) and take its norm.
C
        CALL FX(NVAR,FPAR,IPAR,XR,WK,IERROR)
        IWORK(22)=IWORK(22)+1
        IF(IERROR.NE.0) THEN
          WRITE(LOUNIT,*)
     *    'CORECT - Error flag from user function routine.'
          RETURN
          ENDIF
        WK(NVAR)=XR(IHOLD)-XVALUE
        FNRML=FNRM
        IFMAX=ISAMAX(NVAR,WK,1)
        FNRM=SNRM2(NVAR-1,WK,1)
        IF(IWRITE.GE.2) THEN
          WRITE(LOUNIT,80)IWORK(28),STEPX,KSMAX
          WRITE(LOUNIT,70)IWORK(28),FNRM,IFMAX
          ENDIF
C
C  Check for strong acceptance of function and stepsize.
C
        TLSTEP=ABSERR+RELERR*XNRM
        IF(FNRM.LE.ABSERR.AND.STEPX.LE.TLSTEP)RETURN
C
C  Check for weak acceptance of function and stepsize.
C
        IF(FNRM.LE.EPSATE.OR.STEPX.LE.EPSATE)RETURN
        IF(IWORK(28).GT.1)THEN
          IF((FNRM+FNRML).LE.ABSERR.AND.STEPX.LE.EIGHT*TLSTEP)RETURN
          IF(FNRM.LE.EIGHT*ABSERR.AND.(STEPX+STEPXL).LE.TLSTEP)RETURN
          ENDIF
C
C  Decide if iteration should be aborted
C
        IF(IWORK(28).GT.1)THEN
          IF(ICRIT.LT.1.AND.STEPX.GT.(FMP*STEPXL+ABSERR)) THEN
            IERROR=4
            IF(IWRITE.GE.2)WRITE(LOUNIT,*)
     *      'CORECT - Size of correction not decreasing.'
            RETURN
            ENDIF
          ENDIF
        IF(ICRIT.LT.2.AND.FNRM.GT.(FMP*FNRML+ABSERR)) THEN
           IERROR=4
           IF(IWRITE.GE.2)WRITE(LOUNIT,*)
     *     'CORECT - Residual is not decreasing.'
           RETURN
           ENDIF
        FMP=ONEFIV
  100   CONTINUE
C
C  Reached maximum number of steps without acceptance or rejection.
C
      IERROR=5
      IF(IWRITE.GE.2)WRITE(LOUNIT,*)'CORECT - Convergence too slow.'
      RETURN
70    FORMAT(' CORECT - Residual ',I6,'=',G14.6,' I=',I6)
80    FORMAT(' CORECT - Step     ',I6,15X,G14.6,' I=',I6)
12    FORMAT(' CORECT - Error flag=',I6,' from solver.')
      END
      SUBROUTINE ROOT(A,FA,B,FB,U,FU,KOUNT,IFLAG,IERROR,EPMACH)
C
C***********************************************************************
C
C  ROOT SEEKS A ROOT OF THE EQUATION F(X)=0.0,
C  GIVEN A STARTING INTERVAL (A,B) ON WHICH F CHANGES SIGN.
C  ON FIRST CALL TO ROOT, THE INTERVAL AND FUNCTION VALUES FA AND FB
C  ARE INPUT AND AN APPROXIMATION U FOR THE ROOT IS RETURNED.
C  BEFORE EACH SUBSEQUENT CALL, THE USER EVALUATES FU=F(U), AND THE
C  PROGRAM TRIES TO RETURN A BETTER APPROXIMATION U.
C
C  SEE THE BOOK BY BRENT LISTED in the documentation.
C
C  VARIABLES
C
C  A     = IS ONE ENDPOINT OF AN INTERVAL IN WHICH F CHANGES SIGN.
C  FA    = THE VALUE OF F(A).  THE USER MUST EVALUATE F(A) BEFORE
C          FIRST CALL ONLY.  THEREAFTER THE PROGRAM SETS FA.
C  B     = IS THE OTHER ENDPOINT OF THE INTERVAL IN WHICH
C          F CHANGES SIGN.  NOTE THAT THE PROGRAM WILL RETURN
C          IMMEDIATELY WITH AN ERROR FLAG IF FB*FA.GT.0.0.
C  FB    = THE VALUE OF F(B).  THE USER MUST EVALUATE F(B) BEFORE
C          FIRST CALL ONLY.  THERAFTER THE PROGRAM SETS FB.
C  U     = ON FIRST CALL, U SHOULD NOT BE SET BY THE USER.
C          ON SUBSEQUENT CALLS, U SHOULD NOT BE CHANGED
C          FROM ITS OUTPUT VALUE, THE CURRENT APPROXIMANT
C          TO THE ROOT.
C  FU    = ON FIRST CALL, FU SHOULD NOT BE SET BY THE USER.
C          THEREAFTER, THE USER SHOULD EVALUATE THE FUNCTION
C          AT THE OUTPUT VALUE U, AND RETURN FU=F(U).
C  KOUNT = A COUNTER FOR THE NUMBER OF CALLS TO ROOT.  KOUNT
C          SHOULD BE SET TO ZERO ON THE FIRST CALL FOR A GIVEN
C          ROOT PROBLEM.
C  IFLAG = PROGRAM RETURN FLAG
C          IFLAG=-1  THE CURRENT BRACKETING INTERVAL WITH
C                    ENDPOINTS  STORED IN A AND B IS LESS THAN
C                    4*EPMACH*ABS(U)+EPMACH
C                    HENCE U SHOULD BE ACCEPTED AS THE ROOT.
C                    THE FUNCTION VALUE F(U) IS STORED IN FU.
C          IFLAG= 0  THE INPUT VALUE FU IS EXACTLY ZERO, AND
C                    U SHOULD BE ACCEPTED AS THE ROOT.
C          IFLAG>0   THE CURRENT APPROXIMATION TO THE ROOT IS
C                    CONTAINED IN U.  IF A BETTER APPROXIMATION IS
C                    DESIRED, SET FU=F(U) AND CALL ROOT AGAIN.
C                    THE VALUE OF IFLAG INDICATES
C                    THE METHOD THAT WAS USED TO PRODUCE U.
C
C          IFLAG= 1  BISECTION WAS USED.
C          IFLAG= 2  LINEAR INTERPOLATION (SECANT METHOD).
C          IFLAG= 3  INVERSE QUADRATIC INTERPOLATION.
C
C  IERROR= GLOBAL ERROR FLAG.
C          IERROR=0 NO ERROR FOUND
C          IERROR=7  ON FIRST CALL, FA*FB.GT.0.0. AND HENCE
C                    THE GIVEN INTERVAL IS UNACCEPTABLE.
C  EPMACH= THE RELATIVE MACHINE PRECISION
C
      REAL      EIGHT
      REAL      HALF
      REAL      ONE
      REAL      ONEP5
      REAL      TWO
      REAL      ZERO
C
      PARAMETER (EIGHT=8.0)
      PARAMETER (HALF=0.5)
      PARAMETER (ONE=1.0)
      PARAMETER (ONEP5=1.5)
      PARAMETER (TWO=2.0)
      PARAMETER (ZERO=0.0)
C
      INTRINSIC ABS
      INTRINSIC SIGN
C
      REAL      A
      REAL      B
      REAL      EPMACH
      REAL      FA
      REAL      FB
      REAL      FU
      REAL      HALFUB
      INTEGER   IERROR
      INTEGER   IFLAG
      INTEGER   KOUNT
      REAL      P
      REAL      Q 
      REAL      R
      REAL      S
      REAL      SDEL1
      REAL      SDEL2
      REAL      SDEL3
      REAL      SDEL4
      REAL      STEP
      REAL      TOLER
      REAL      U
C
C  SEGMENT 1   FIRST CALL HANDLED SPECIALLY.  DO BOOKKEEPING.
C
      IF(KOUNT.LE.0)THEN
        IF((FA.GT.ZERO.AND.FB.GT.ZERO) 
     1  .OR.(FA.LT.ZERO.AND.FB.LT.ZERO))THEN
          IERROR=7
          KOUNT=0
          RETURN
          ENDIF
        KOUNT=1
        SDEL1=TWO*ABS(B-A)
        SDEL2=TWO*SDEL1
        SDEL3=TWO*SDEL2
        U=B
        B=A
        FU=FB
        FB=FA
      ELSE
C
C  INCREMENT COUNTER AND CHECK WHETHER F(U)=0.
C
        KOUNT=KOUNT+1
        IF(FU.EQ.ZERO)THEN
          IFLAG=0
          RETURN
          ENDIF
C
C  IF FU AND FB HAVE THE SAME SIGN OVERWRITE B WITH A
C
        IF(SIGN(ONE,FU).EQ.SIGN(ONE,FB))THEN
          B=A
          FB=FA
          ENDIF
        ENDIF
C
C  SEGMENT 2   REARRANGE POINTS AND FUNCTION VALUES IF
C  NECESSARY SO THAT  ABS(FU).LT.ABS(FB)
C
      IF(ABS(FB).LT.ABS(FU)) THEN
        A=U
        U=B
        B=A
        FA=FU
        FU=FB
        FB=FA
        ENDIF
C
C  SEGMENT 3   CHECK FOR ACCEPTANCE BECAUSE OF SMALL INTERVAL
C  CURRENT CHANGE IN SIGN INTERVAL IS (B,U) OR (U,B).
C
      TOLER=TWO*EPMACH*ABS(U)+EPMACH
      HALFUB=HALF*(B-U)
      SDEL4=SDEL3
      SDEL3=SDEL2
      SDEL2=SDEL1
      SDEL1=ABS(B-U)
      IF(ABS(HALFUB).LE.TOLER) THEN
        IFLAG=-1
        A=U
        FA=FU
        RETURN
        ENDIF
C
C  SEGMENT 4   COMPUTE NEW APPROXIMANT TO ROOT OF THE FORM
C  U(NEW)=U(OLD)+STEP.
C  METHODS AVAILABLE ARE LINEAR INTERPOLATION
C  INVERSE QUADRATIC INTERPOLATION
C  AND BISECTION.
C
      IF(ABS(FU).GE.ABS(FA)) THEN
        IFLAG=1
        STEP=HALFUB
        GO TO 80
        ENDIF
C
C  ATTEMPT LINEAR INTERPOLATION IF ONLY TWO POINTS AVAILABLE
C  COMPUTE P AND Q FOR APPROXIMATION U(NEW)=U(OLD)+P/Q
C
      IF(A.EQ.B) THEN
        IFLAG=2
        S=FU/FA
        P=TWO*HALFUB*S
        Q=ONE-S
C    
C  ATTEMPT INVERSE QUADRATIC INTERPOLATION IF THREE POINTS AVAILABLE
C  COMPUTE P AND Q FOR APPROXIMATION U(NEW)=U(OLD)+P/Q
C
      ELSE
        IFLAG=3
        S=FU/FA
        Q=FA/FB
        R=FU/FB
        P=S*(TWO*HALFUB*Q*(Q-R)-(U-A)*(R-ONE))
        Q=(Q-ONE)*(R-ONE)*(S-ONE)
        ENDIF
C
C  CORRECT THE SIGNS OF P AND Q
C
      IF(P.GT.ZERO)Q=-Q
      P=ABS(P)
C
C  IF P/Q IS TOO LARGE, GO BACK TO BISECTION
C
      IF((EIGHT*SDEL1.GT.SDEL4) 
     1 .OR.(P.GE.ONEP5*ABS(HALFUB*Q)-ABS(TOLER*Q))) THEN
        IFLAG=1
        STEP=HALFUB
        GO TO 80
        ENDIF
      STEP=P/Q
C
C  SEGMENT 5   VALUE OF STEP HAS BEEN COMPUTED.
C  UPDATE INFORMATION  A =U, FA =FU, U =U+STEP.
C  CHANGE IN SIGN INTERVAL IS NOW (A,B) OR (B,A).
C
80    CONTINUE
      A=U
      FA=FU
      IF(ABS(STEP).LE.TOLER) STEP=SIGN(TOLER,HALFUB)
      U=U+STEP
      RETURN
      END
      SUBROUTINE SETSTP(HTAN,IWORK,LIW,LRW,RWORK)
C
C***********************************************************************
C
C  Predictor stepsize computation.
C
C  THE FORMULAS UNDERLYING THE ALGORITHM ARE
C
C  ALPHLC=MAXIMUM OF ARCCOS(TL,TC) AND ALFMIN=2*ARCCOS(1-EPMACH)
C  HSEC  =NORM(XC-XF)
C  HSECL =NORM(XL-XC)
C  ABSNLC=ABS(SIN(.5*ALPHLC))
C  CURV  =LAST VALUE OF THE CURVATURE
C  CURVN =2*ABSNLC/HSEC
C  CORDIS=CORRECTOR DISTANCE TO CURRENT CONTINUATION POINT.
C           FORCE CORDIS TO LIE BETWEEN .01*HSEC AND HSEC.
C           UNLESS (CORDIS.EQ.0.0), BECAUSE THE PREDICTED POINT WAS
C           ACCEPTED.  IN SUCH A CASE, SET HTAN=HFACT*HSEC
C           INSTEAD OF USING FIRST ESTIMATE FOR HTAN.
C
C  THEN
C
C  CURVX=CURVN+HSEC*(CURVN-CURV)/(HSEC+HSECL)
C  A SIMPLER FORMULA IS USED IF WE DO NOT HAVE DATA AT TWO PREVIOUS
C  POINTS.
C
C  IF (IWORK(10).GE.2) CURVX MUST BE AT LEAST AS LARGE AS
C  THE MAXIMUM OF 0.001 AND .01/HSEC
C
C  FIRST ESTIMATE   (UNLESS (CORDIS.EQ.0.0) )
C
C  HTAN =SQRT(2*QUAL*CORDIS/CURVX)
C
C  ADJUSTED VALUE
C
C  HTAN=HTAN*(1.0+HTAN*(TC(IPC)-TL(IPC))/(2*HSEC*TC(IPC)))
C
C  READJUSTMENT AND TRUNCATIONS
C
C  IF STEPSIZE REDUCTION OCCURRED DURING LAST CORRECTOR PROCESS,
C  HTAN IS FORCED TO BE LESS THAN (HFACT-1)*HSEC/2.
C
C  HTAN MUST LIE BETWEEN (HSEC/HFACT) AND (HSEC*HFACT).
C  HTAN IS ALWAYS FORCED TO LIE BETWEEN USER SPECIFIED
C  MINIMUM AND MAXIMUM STEPS.
C
      REAL      HALF
      REAL      HUNDTH
      REAL      ONE
      REAL      THOUTH
      REAL      TWO
      REAL      ZERO
C
      PARAMETER (HALF=0.5)
      PARAMETER (HUNDTH=0.01)
      PARAMETER (ONE=1.0)
      PARAMETER (THOUTH=0.001)
      PARAMETER (TWO=2.0)
      PARAMETER (ZERO=0.0)
C
      INTRINSIC ABS
      INTRINSIC MAX
      INTRINSIC MIN
      INTRINSIC SIN
      INTRINSIC SQRT
C
      INTEGER   LIW
      INTEGER   LRW
C
      REAL      CURV
      REAL      CURVX
      REAL      HTAN
      INTEGER   IWORK(LIW)
      REAL      RWORK(LRW)
      REAL      TEMP
      REAL      TMP
C
      RWORK(11)=MAX(RWORK(10),RWORK(11))
C
C  COMPUTE NEW CURVATURE DATA
C
      CURV=RWORK(16)
      RWORK(16)=TWO*ABS(SIN(HALF*RWORK(11)))/RWORK(21)
      IF(CURV.EQ.ZERO)CURV=RWORK(16)
      CURVX=RWORK(16)
      IF(RWORK(22).NE.ZERO)CURVX=
     1 RWORK(16)+RWORK(21)*(RWORK(16)-CURV)/(RWORK(21)+RWORK(22))
      TMP=MAX(THOUTH,(HUNDTH)/RWORK(21))
      CURVX=MAX(CURVX,TMP)
C
C  IF CONVERGENCE DISTANCE IS ZERO, SET ESTIMATE TO HFACT*HSEC.
C  ELSE TRUNCATE QUAL*CORDIS TO LIE BETWEEN .01*HSEC AND HSEC.
C
      HTAN=RWORK(20)*RWORK(21)
      IF(RWORK(15).NE.ZERO)THEN
        TEMP=MAX(RWORK(23)*RWORK(15),HUNDTH*RWORK(21))
        TEMP=MIN(TEMP,RWORK(21))
        HTAN=SQRT(TWO*TEMP/CURVX)
        ENDIF
C
C  Adjust the step to account for estimated curvature in the direction
C  of the parameter.
C
      IF(IWORK(18).GT.0)
     * HTAN=MIN(HTAN,HALF*(RWORK(20)-ONE)*RWORK(21))
      IF(IWORK(3).NE.1)THEN
        TEMP=ONE+(ONE-RWORK(25)/RWORK(24))*(HALF*HTAN)/RWORK(21)
        HTAN=HTAN*TEMP
        ENDIF
C
C  Enforce restrictions on growth and size of the step.
C
      HTAN=MAX(HTAN,RWORK(21)/RWORK(20))
      HTAN=MIN(HTAN,RWORK(21)*RWORK(20))
      HTAN=MAX(HTAN,RWORK(3))
      HTAN=MIN(HTAN,RWORK(4))
      RETURN
      END
      SUBROUTINE START(DF,FPAR,FX,IERROR,IPAR,IPC,IWORK,IWRITE,LIW,
     *LOUNIT,LRW,NVAR,RWORK,TC,WORK1,XC,XF,XR,SLNAME)
C
C***********************************************************************
C
C  On first call, make sure that the input point satisfies the
C  nonlinear equations.  If not, call CORECT and attempt to force this.
C  
      REAL      ONE
      REAL      ZERO
C
      PARAMETER (ONE=1.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  COQUAL
      EXTERNAL  CORECT
      EXTERNAL  DF
      EXTERNAL  FX
      EXTERNAL  ISAMAX
      EXTERNAL  SAXPY
      EXTERNAL  SCOPY
      EXTERNAL  SLNAME
C
      INTRINSIC ABS
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      REAL      DETS
      REAL      FPAR(*)
      INTEGER   I
      INTEGER   ICRIT
      INTEGER   IERROR
      INTEGER   IMAX
      INTEGER   IPAR(*)
      INTEGER   IPC
      INTEGER   ISAMAX
      INTEGER   IWORK(LIW)
      INTEGER   IWRITE
      INTEGER   JOB
      INTEGER   LOUNIT
      INTEGER   MODSAV
      REAL      RWORK(LRW)
      REAL      SKALE
      REAL      STEPX
      REAL      TC(NVAR)
      REAL      WORK1(NVAR)
      REAL      XC(NVAR)
      REAL      XF(NVAR)
      REAL      XR(NVAR)
C
C  If user is requesting that Jacobian be used as long as possible,
C  then go ahead, generate and factor the first one now.
C
      IF(IWORK(4).EQ.2)THEN
        JOB=2
        CALL SLNAME(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,LIW,
     1  JOB,NVAR,RWORK,LRW,XR,WORK1)
        RWORK(17)=DETS
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)
     *    'START  - Could not factor initial jacobian.'
          RETURN
          ENDIF
        ENDIF
      IF(IWRITE.GE.2)WRITE(LOUNIT,1040)IPC
1040  FORMAT(' START  - Correct initial point, fixing index ',I5)
C
C  Set up a pseudo-tangent vector and other data for first time,
C  check that starting point satisfies the equations.
C
      DO 50 I=1,NVAR
        TC(I)=ZERO
50      CONTINUE
      CALL SCOPY(NVAR,XR,1,XC,1)
      TC(IPC)=ONE
      MODSAV=IWORK(4)
      ICRIT=1
C
60    CONTINUE
      CALL SCOPY(NVAR,XC,1,XR,1)
      CALL CORECT(DF,FPAR,FX,IERROR,IPC,IPAR,IWORK,NVAR,RWORK,STEPX,
     *WORK1,XR,LRW,LIW,ICRIT,SLNAME)
      IWORK(25)=IWORK(25)+IWORK(28)
C
C  If possible, retry with ICRIT=2.
C
      IF(IERROR.NE.0.AND.ICRIT.EQ.1)THEN
        IF(IWRITE.GE.1)WRITE(LOUNIT,*)
     *  'START -  Retry starting point correction'
        ICRIT=2
        GO TO 60
      ELSEIF(IERROR.NE.0)THEN
        ICRIT=1
        ENDIF
C
      IF(IERROR.NE.0.AND.IWORK(4).GT.0)THEN
        IWORK(4)=IWORK(4)-1
        IERROR=0
1110  FORMAT(' START  - Retrying starting point with IWORK(4)=',I6)
        IF(IWRITE.GE.1)WRITE(LOUNIT,1110)IWORK(4)
        GO TO 60
        ENDIF
      IWORK(4)=MODSAV
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,*)
     *  'START  - Starting point correction failed.'
        RETURN
        ENDIF
      SKALE=-ONE
      CALL SAXPY(NVAR,SKALE,XR,1,XC,1)
      IMAX=ISAMAX(NVAR,XC,1)
      RWORK(15)=ABS(XC(IMAX))
      CALL SCOPY(NVAR,XR,1,XC,1)
      CALL SCOPY(NVAR,XR,1,XF,1)
      CALL COQUAL(STEPX,IWORK,LIW,RWORK,LRW)
      RWORK(14)=RWORK(13)
      IWORK(27)=IWORK(27)+1
      IWORK(10)=1
      IWORK(1)=1
      RETURN
      END
      SUBROUTINE TANGNT(DETSN,FX,DF,FPAR,IERROR,IP,IPAR,IWORK,
     1 NVAR,RWORK,TAN,XR,LIW,LRW,SLNAME)
C
C***********************************************************************
C
C  TANGNT COMPUTES A UNIT TANGENT VECTOR TO THE SOLUTION
C  CURVE OF THE UNDERDETERMINED NONLINEAR SYSTEM FX = 0.  THE
C  TANGENT VECTOR TAN IS THE SOLUTION OF THE LINEAR SYSTEM
C
C        DFA(X,IPL)*TAN = E(NVAR)
C
C  HERE, E(I) ALWAYS DENOTES THE I-TH BASIS NATURAL BASIS VECTOR
C  IN NVAR-SPACE; THAT IS, THE VECTOR WITH A 1 IN THE I-TH
C  COMPONENT AND ZEROS ELSEWHERE. THEN DFA(X,IPL) IS THE
C  AUGMENTED JACOBIAN WHOSE FIRST NVAR-1 ROWS ARE DFX/DX (X), THE
C  DERIVATIVE OF FX EVELUATED AT X, AND WHOSE LAST ROW IS
C  THE TRANSPOSE OF E(IPL).
C
C  THE TANGENT VECTOR IS THEN NORMALIZED, BUT THE ADJUSTMENT
C  OF ITS SIGN IS PERFORMED OUTSIDE THE ROUTINE.
C
C  ERROR RETURNS  IERROR=1  DATA OR STORAGE ERROR
C                 IERROR=2  ERROR IN DERIVATIVE CALL
C                 IERROR=3  SOLVER FAILED
C                 IERROR=6  TANGENT VECTOR ZERO
C
      REAL      ONE
      REAL      ZERO
C
      PARAMETER (ONE=1.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  FX
      EXTERNAL  DF
      EXTERNAL  SLNAME
      EXTERNAL  SNRM2
      EXTERNAL  SSCAL
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      REAL      DETSN
      REAL      FPAR(*)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IP
      INTEGER   IPAR(*)
      INTEGER   IWORK(LIW)
      INTEGER   JOB
      REAL      RWORK(LRW)
      REAL      SKALE
      REAL      SNRM2
      REAL      TAN(NVAR)
      REAL      TNORM
      REAL      XR(NVAR)
C
C  SET RIGHT HAND SIDE OF TANGENT SYSTEM
C
      DO 10 I=1,NVAR
        TAN(I)=ZERO
   10   CONTINUE
      TAN(NVAR)=ONE
C
C  CALL SOLVER
C
      JOB=0
      IF(IWORK(4).EQ.2)JOB=1
      CALL SLNAME(DETSN,FX,DF,FPAR,IERROR,IP,IPAR,IWORK,LIW,
     1JOB,NVAR,RWORK,LRW,XR,TAN)
      IF(IERROR.NE.0)RETURN
C
C  Normalize the tangent vector.
C
      TNORM=SNRM2(NVAR,TAN,1)
      IF(TNORM.EQ.ZERO)THEN
        IERROR=6
      ELSE
        SKALE=ONE/TNORM
        CALL SSCAL(NVAR,SKALE,TAN,1)
        IERROR=0
        ENDIF
      RETURN
      END
      SUBROUTINE DENSLV(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,
     1 IWORK,LIW,JOB,NVAR,RWORK,LRW,X,Y)
C
C***********************************************************************
C
C  DENSLV SOLVES THE NVAR*NVAR LINEAR SYSTEM
C
C               (  DF(X)  )
C  (1)          (         ) * Y   = Y
C               (  E(IPC) )
C
C  WHERE DF(X) IS THE JACOBIAN OF THE GIVEN FUNCTION FX AT X, AND
C  E(IPC) IS THE IP-TH NATURAL BASIS VECTOR; THAT IS, THE LAST ROW
C  CONSISTS OF ALL ZEROS EXCEPT FOR A 1 IN COLUMN IPC.
C
C  THE MATRIX IS STORED IN DENSE FORM IN THE ONE-DIMENSIONAL ARRAY
C  RWORK(LRBEG)-RWORK(LRLST) WHERE LRLST=LRBEG+NVAR*NVAR - 1 AND
C  LRBEG IS STORED IN IWORK(15). A DATA ERROR IS CAUSED IF LRLST.GT.LRW
C
C  THE ROUTINE STORES THE PIVOT ARRAY IN IWORK(IWORK(13))) THROUGH 
C  IWORK(LILST) WHERE LILST=IWORK(13)+NVAR - 1.
C  THE ROUTINE RETURNS WITH A DATA ERROR IF LILST.GT.LIW.
C
C  SINCE THE MATRIX OF THE SYSTEM (1) IS STORED IN DENSE FORM THE USER
C  PROGRAM DF MAY HAVE THE FORM
C
C         SUBROUTINE DF(NVAR,FPAR,IPAR,X,A,IERR)
C         DIMENSION A(NVAR,NVAR), X(NVAR), FPAR(*),IPAR(*)
C         .....
C         A(I,J)=DF(I)/DX(J) , I=1,...,NVAR-1,J=1,...,NVAR
C         .....
C
C  WHERE DF(I)/DX(J) DENOTES THE DERIVATIVE OF THE I-TH COMPONENT
C  OF THE FUNCTION FX WITH RESPECT TO THE J-TH VARIABLE.  THE USER
C  NEED NOT SUPPLY THE LAST ROW OF THE MATRIX OF THE SYSTEM (1).
C
C  THE LINPACK SUBROUTINES SGEFA, SGEDI, AND SGESL ARE USED.
C
C  LIST OF VARIABLES
C
C  DETS    THE SIGN OF THE DETERMINANT OF THE MATRIX IN (1).
C  FX      AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED
C          SUBROUTINE WHICH COMPUTES THE NVAR-1 DIMENSIONAL
C          FUNCTION.
C
C  DF      AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED
C          SUBROUTINE WHICH COMPUTES THE JACOBIAN OF THE FUNCTION.
C  FPAR    A REAL ARRAY FOR THE TRANSMISSION OF PARAMETERS TO
C          DENSLV AND DF.
C  IERROR  ERROR FLAG
C         = 0   NORMAL RETURN
C         = 1   DATA OR STORAGE ERROR
C         = 2   ERROR IN DERIVATIVE ROUTINE DF
C         = 3   NUMERICALLY SINGULAR MATRIX DETECTED.
C  IPC     THE LOCATION OF THE 1 IN THE LAST ROW OF THE MATRIX
C  IPAR    AN INTEGER ARRAY FOR TRANSMISSION OF PARAMETERS TO DENSLV
C          AND DF.
C  IWORK   INTEGER WORK ARRAY USED TO STORE THE PIVOT ARRAY
C          (SEE ABOVE)
C  LIW     DIMENSION OF IWORK IN THE CALLING PROGRAM
C  JOB     WORK REQUEST SWITCH
C          = 0   DECOMPOSE THE MATRIX, SOLVE THE SYSTEM, AND
C          COMPUTE THE SIGN OF THE DETERMINANT.
C          = 1   SOLVE THE SYSTEM, ASSUMING THAT THE MATRIX HAS
C          BEEN PREVIOUSLY DECOMPOSED.
C          2 = Factor system, compute determinant.
C          3 = Check jacobian matrix.  Call user jacobian routine,
C          multiply by -1.0, add finite difference jacobian,
C          print largest entry.
C
C  NVAR    THE DIMENSION OF THE SYSTEM.
C  RWORK   REAL WORK ARRAY USED TO STORE THE MATRIX (SEE ABOVE).
C  LRW     THE DIMENSION OF RWORK IN THE CALLING PROGRAM.
C  X       A REAL VECTOR OF LENGTH NVAR, THE POINT AT WHICH THE
C          JACOBIAN IS TO BE EVALUATED.
C  Y       A REAL VECTOR OF LENGTH NVAR. ON INPUT Y CONTAINS THE
C          RIGHT HAND SIDE OF THE LINEAR SYSTEM (1), ON RETURN
C          WITH IERROR = 0 THE SOLUTION OF THE SYSTEM IS GIVEN IN Y.
C
      REAL      ONE
      REAL      ZERO
C
      PARAMETER (ONE=1.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  DENJAC
      EXTERNAL  DF
      EXTERNAL  FX
      EXTERNAL  ISAMAX
      EXTERNAL  SGEDI
      EXTERNAL  SGEFA
      EXTERNAL  SGESL
      EXTERNAL  SSCAL
C
      INTRINSIC MOD
      INTRINSIC SQRT
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      REAL      DET(2)
      REAL      DETS
      REAL      EPS
      REAL      FPAR(*)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(*)
      INTEGER   IPC
      INTEGER   ISAMAX
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   JAC
      INTEGER   JACK
      INTEGER   JOB
      INTEGER   JOB2 
      INTEGER   K
      INTEGER   LDF
      INTEGER   LFXM
      INTEGER   LFXP
      INTEGER   LILST
      INTEGER   LOUNIT
      INTEGER   LPIV
      INTEGER   LRLST
      INTEGER   NDIM
      REAL      RWORK(LRW)
      REAL      SKALE
      REAL      X(NVAR)
      REAL      Y(NVAR)
C
C  CHECK INPUT VALUES OF NVAR, AND IPC.
C  DEPENDING ON VALUE OF JOB, EITHER SET UP
C  AUGMENTED JACOBIAN, DECOMPOSE INTO L-U FACTORS, AND GET DETERMINANT,
C  OR USE CURRENT FACTORED JACOBIAN WITH NEW RIGHT HAND SIDE.
C
      LOUNIT=IWORK(8)
      IERROR=0
      LPIV=IWORK(13)
      LDF=IWORK(15)
      LILST=LPIV+NVAR-1
      JAC=IWORK(9)
      IF(JAC.NE.0.AND.JOB.EQ.3)THEN
        IERROR=4
        WRITE(LOUNIT,*)'DENSLV - Error!  Jacobian check requested'
        WRITE(LOUNIT,*)'         but no user Jacobian routine.'
        RETURN
        ENDIF
      IF(JAC.EQ.0.AND.JOB.NE.3)THEN
        LRLST=LDF-1+NVAR*NVAR
      ELSE
        LFXP=LDF+NVAR*NVAR
        LFXM=LFXP+NVAR
        LRLST=LDF-1+NVAR*NVAR+2*NVAR
        ENDIF
      IF((LILST.GT.LIW).OR.(LRLST.GT.LRW))THEN
        IERROR=1
        WRITE(LOUNIT,1040)LILST,LIW
        WRITE(LOUNIT,1050)LRLST,LRW
        RETURN
        ENDIF
      IF(JOB.EQ.1)GO TO 10
C
C  IF JOB IS 0 OR 2, THEN WE MUST EVALUATE THE JACOBIAN, USING THE USER
C  ROUTINE OR FINITE DIFFERENCES, AND THEN FACTOR IT AND GET THE
C  SIGN OF THE DETERMINANT.
C
      NDIM=NVAR*NVAR
      DO 5 I=1,NDIM
        RWORK(LDF+I-1)=ZERO
5       CONTINUE
      IF(JAC.EQ.0)THEN
        CALL DF(NVAR,FPAR,IPAR,X,RWORK(LDF),IERROR)
        IWORK(19)=IWORK(19)+1
        RWORK(LDF+IPC*NVAR-1)=ONE
        ENDIF
      IF(JOB.EQ.3)THEN
        SKALE=-ONE
        CALL SSCAL(NDIM,SKALE,RWORK(LDF),1)
        ENDIF
      IF(JAC.EQ.1.OR.JAC.EQ.2.OR.JOB.EQ.3)THEN
        JACK=JAC
        IF(JOB.EQ.3)JACK=2
        EPS=SQRT(SQRT(RWORK(8)))
        CALL DENJAC(EPS,FPAR,RWORK(LDF),FX,IERROR,IPAR,IPC,IWORK,
     *  JACK,LIW,LOUNIT,NVAR,X,RWORK(LFXP),RWORK(LFXM))
        ENDIF
      IF(IERROR.NE.0)RETURN
C
C  If JOB=3, print out largest entry of matrix.
C
      IF(JOB.EQ.3)THEN
        K=ISAMAX(NDIM,RWORK(LDF),1)
        I=MOD(K-1,NVAR)+1
        J=(K-I)/NVAR+1
        WRITE(LOUNIT,1070)RWORK(LDF+K-1),I,J
        IF(IWORK(1).EQ.(-2))THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'DENSLV - Entire difference matrix:'
          WRITE(LOUNIT,*)' '
          DO 99 I=1,NVAR
            DO 98 J=1,NVAR
              K=LDF+(J-1)*NVAR+I-1
              WRITE(LOUNIT,1080)RWORK(K),I,J
98            CONTINUE
            WRITE(LOUNIT,*)' '
99          CONTINUE
          ENDIF
        RETURN
        ENDIF
C
C  DECOMPOSE MATRIX
C
      CALL SGEFA(RWORK(LDF),NVAR,NVAR,IWORK(LPIV),IERROR)
      IWORK(20)=IWORK(20)+1
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,1492)IERROR
        IERROR=3
        RETURN
        ENDIF
C
C  COMPUTE DETERMINANT
C
      JOB2=10
      CALL SGEDI(RWORK(LDF),NVAR,NVAR,IWORK(LPIV),DET,Y,JOB2)
C
C  RECORD THE SIGN OF THE DETERMINANT
C
      DETS=ZERO
      IF(DET(1).GT.ZERO)THEN
        DETS=ONE
      ELSEIF(DET(1).LT.ZERO)THEN
        DETS=-ONE
        ENDIF
      IF(JOB.EQ.2)RETURN
C
C  USING L-U FACTORED MATRIX, SOLVE SYSTEM USING FORWARD-BACKWARD
C  ELIMINATION, AND OVERWRITE RIGHT HAND SIDE WITH SOLUTION
C
10    CONTINUE
      CALL SGESL(RWORK(LDF),NVAR,NVAR,IWORK(LPIV),Y,0)
      IWORK(21)=IWORK(21)+1
      RETURN
1040  FORMAT(' DENSLV - Need LIW=',I6,', have LIW=',I6)
1050  FORMAT(' DENSLV - Need LRW=',I6,', have LRW=',I6)
1070  FORMAT(' DENSLV - Maximum difference between user and estimated'/
     *       '          jacobian is ',G14.6,' row ',I6,' column ',I6)
1080  FORMAT(1X,G14.6,' =FP(I,J)-DELF(I,J), I, J=',2I6)
1492  FORMAT(' DENSLV - Zero pivot, SGEFA returns INFO=',I6)
      END
      SUBROUTINE DENJAC(EPS,FPAR,FPRIME,FX,IERROR,IPAR,IPC,IWORK,
     * JAC,LIW,LOUNIT,NVAR,X,WORK1,WORK2)
C
C************************************************************************
C
C  DENJAC ESTIMATES THE JACOBIAN FPRIME OF THE FUNCTION FX ASSUMING
C  THAT FULL (DENSE) STORAGE IS USED.  IT IS CALLED AUTOMATICALLY
C  BY DENSLV WHEN JAC=1 OR 2.
C
C  EPS    - A TOLERANCE TO USE FOR SHIFTING THE X VALUES.  
C           EPS=SQRT(MACHINE EPSILON) FOR THE VAX, OR
C           EPS=SQRT(SQRT(MACHINE EPSILON)) FOR THE PC 
C           HAVE BEEN FOUND SUITABLE.
C
C  FPAR   - USER REAL PARAMETER VECTOR, PASSED TO FX BUT NOT OTHERWISE
C           REFERENCED.
C
C  FPRIME - REAL(NVAR,NVAR), THE ARRAY INTO WHICH THE JACOBIAN WILL
C           BE STORED, INCLUDING THE LAST AUXILLIARY ROW.
C
C  FX     - EXTERNAL NAME OF USER ROUTINE DEFINING FUNCTION, OF THE
C           FORM SUBROUTINE FX(NVAR,FPAR,IPAR,X,F,IERROR).
C
C  IERROR - ERROR RETURN FLAG.  IERROR NONZERO MEANS THERE WAS A PROBLEM
C           IN THE USER ROUTINE FX OR IN DENJAC.
C
C  IPAR   - USER INTEGER PARAMETER VECTOR, PASSED TO FX BUT NOT OTHERWISE
C           REFERENCED.
C
C  IPC    - INDEX OF CURRENT CONTINUATION PARAMETER, DETERMINES FORM OF
C           LAST ROW OF JACOBIAN.
C
C  IWORK  - INTEGER WORK VECTOR OF DIMENSION LIW.  ONLY REQUIRED IN ORDER
C           TO UPDATE IWORK(22), WHICH COUNTS FUNCTION EVALUATIONS.
C
C  JAC    - USER DEFINED JACOBIAN OPTION.  0 MEANS USER SUPPLIES JACOBIAN.
C           1 MEANS USE FORWARD, 2 MEANS USE CENTRAL DIFFERENCES TO
C           ESTIMATE JACOBIAN FROM FX.
C
C  LIW    - DIMENSION OF IWORK.
C
C  LOUNIT - OUTPUT UNIT.
C
C  NVAR   - DIMENSION OF X, FPRIME, WORK1, WORK2.  THE NUMBER OF VARIABLES.
C
C  X      - POINT AT WHICH JACOBIAN IS TO BE ESTIMATED.
C
C  WORK1  - WORK VECTOR OF DIMENSION NVAR.
C
C  WORK2  - WORK VECTOR OF DIMENSION NVAR.
C
      REAL      ONE
      REAL      ZERO
C
      PARAMETER (ONE=1.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  FX
      EXTERNAL  SAXPY
      EXTERNAL  SSCAL
C
      INTRINSIC ABS
C
      INTEGER   LIW
      INTEGER   NVAR
C
      REAL      DELM
      REAL      DELP
      REAL      EPS
      REAL      FPAR(*)
      REAL      FPRIME(NVAR,NVAR)
      INTEGER   IPAR(*)
      INTEGER   IERROR
      INTEGER   IPC
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   JAC
      INTEGER   LOUNIT
      REAL      SKALE
      REAL      X(NVAR)
      REAL      XSAVE
      REAL      WORK1(NVAR)
      REAL      WORK2(NVAR)
C
      IF(JAC.EQ.1)THEN
        CALL FX(NVAR,FPAR,IPAR,X,WORK2,IERROR)
        DELM=ZERO
        IWORK(22)=IWORK(22)+1
        IF(IERROR.NE.0)RETURN
        ENDIF
C
      DO 20 J=1,NVAR
        XSAVE=X(J)
        DELP=EPS*(ONE+ABS(X(J)))
        X(J)=X(J)+DELP
        CALL FX(NVAR,FPAR,IPAR,X,WORK1,IERROR)
        IWORK(22)=IWORK(22)+1
        IF(IERROR.NE.0)RETURN
        IF(JAC.EQ.2)THEN
          DELM=-DELP
          X(J)=XSAVE+DELM
          CALL FX(NVAR,FPAR,IPAR,X,WORK2,IERROR)
          IWORK(22)=IWORK(22)+1
          IF(IERROR.NE.0)RETURN
          ENDIF
        X(J)=XSAVE
C
C  Compute DFDX(*,J)= (F(X+)-F(X-))/DELX.
C  Note that we contrive to ADD this quantity to DFDX,
C  This makes it possible to use this routine also to check
C  the accuracy of the user's jacobian routine.
C
        SKALE=-ONE
        CALL SAXPY(NVAR-1,SKALE,WORK2,1,WORK1,1)
        SKALE=ONE/(DELP-DELM)
        CALL SSCAL(NVAR-1,SKALE,WORK1,1)
        SKALE=ONE
        CALL SAXPY(NVAR-1,SKALE,WORK1,1,FPRIME(1,J),1)
20      CONTINUE
C
      FPRIME(NVAR,IPC)=FPRIME(NVAR,IPC)+ONE
      RETURN
      END
      SUBROUTINE BANSLV(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,LIW,
     1 JOB,NVAR,RWORK,LRW,X,Y)
C
C***********************************************************************
C
C  BANSLV SOLVES A LINEAR SYSTEM OF THE FORM
C
C  ( DFDY  DFDZ )
C  (            ) * X=B
C  (   E(IPC)   )
C
C  WHERE DFDY IS FORMALLY AN (NVAR-1) BY (NVAR-1) MATRIX WHICH IS
C  TO BE STORED IN LINPACK GENERAL BAND STORAGE, DFDZ IS A COLUMN
C  VECTOR, AND E(IPC) IS AN NVAR-DIMENSIONAL ROW VECTOR WHICH IS 0
C  EXCEPT FOR A 1 IN THE IPC-TH POSITION.  IT IS ASSUMED THAT THE
C  FULL MATRIX ABOVE REPRESENTS THE JACOBIAN OF A SET OF NONLINEAR
C  EQUATIONS, WITH THE LAST EQUATION BEING X(IPC)=CONSTANT.
C
C  THE PURPOSE OF THIS ROUTINE IS TO SOLVE THE FULL SYSTEM, WHILE
C  TAKING ADVANTAGE OF SPECIAL PROPERTIES OF THE (NVAR-1) SUBSYSTEM
C  (IN THIS CASE BANDEDNESS), WHICH WOULD BE LOST IF THE FULL
C  SYSTEM WAS SOLVED DIRECTLY.
C
C  SEE THE PAPER BY CHAN LISTED ABOVE FOR DETAILS.
C
C  SUBROUTINE ARGUMENTS
C
C  DETS   - THE SIGN OF THE DETERMINANT OF THE FULL MATRIX.
C
C  FX     - AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED
C           SUBROUTINE WHICH COMPUTES THE NVAR-1 DIMENSIONAL
C           FUNCTION.
C
C  DF     - AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED
C           SUBROUTINE WHICH COMPUTES THE JACOBIAN OF THE NVAR-1
C           NONLINEAR FUNCTIONS.  THE LAST, AUGMENTING FUNCTION,
C           IS HANDLED BY THIS ROUTINE.  BECAUSE OF THE SPECIAL
C           STORAGE ARRANGEMENTS, GREAT CARE MUST BE TAKEN TO
C           STORE INFORMATION PROPERLY.
C
C           LET MU AND ML BE THE UPPER AND LOWER BANDWIDTHS.
C           THESE VALUES MUST BE STORED IN IPAR(*) AND IPAR(2),
C           RESPECTIVELY. SET NBAND=2*ML+MU+1.
C           THE INDEX K OF RWORK IN WHICH WE ARE TO STORE THE
C           (I,J) ENTRY OF THE JACOBIAN IS DETERMINED AS FOLLOWS:
C
C           IF(J.EQ.NVAR) K=(NVAR-1)*NBAND+I
C
C           OTHERWISE
C
C           IF((J-I.LE.ML).AND.(I-J.LE.MU)) K=I+J*(NBAND-1)-ML
C
C           FOR ALL OTHER VALUES OF I AND J THE ENTRY IS ZERO AND
C           NOT STORED.
C
C           THE FORM OF THE CALLING SEQUENCE OF DF IS
C           SUBROUTINE DF(NVAR,FPAR,IPAR,X,RWORK,IERR)
C           DIMENSION FPAR(*),IPAR(*),X(NVAR),RWORK(1)
C
C           SET IERR NONZERO IF EXECUTION IS TO BE HALTED FOR ANY REASON
C
C  FPAR   - A REAL ARRAY FOR TRANSMISSION OF PARAMETERS FROM THE
C           USER CALLING PROGRAM THROUGH THE CONTINUATION CODE TO
C           USER SUBROUTINES.  THE CONTINUATION CODE ASSUMES A
C           DIMENSION OF 1 FOR FPAR, AND NEVER ACCESSES ANY ENTRIES.
C
C  IERROR - ERROR RETURN FLAG.
C           IERROR= 0, NO ERRORS DETECTED.
C           IERROR= 1, DATA OR STORAGE ERROR.
C                      ILLEGAL VALUES OF NVAR, IPC, MU, ML, OR
C                      INSUFFICIENT STORAGE IN RWORK OR IWORK.
C           IERROR= 2, ERROR RETURN FROM DF.
C                      THE USER ROUTINE DF HAS SET AN ERROR FLAG.
C           IERROR= 3, NUMERICALLY SINGULAR MATRIX DETECTED.
C                      EITHER THE SUBMATRIX AS MODIFIED BY THIS
C                      ROUTINE IS SINGULAR, AND THE DECOMPOSITION
C                      OR THE CHOICE OF IPC IS INAPPROPRIATE,
C                      OR THE FULL SYSTEM MAY ACTUALLY BE SINGULAR.
C
C  IPC    - AN INDEX DETERMINED BY THE CONTINUATION CODE, WHICH
C           DEFINES THE AUGMENTING EQUATION, AND HENCE THE FORM
C           OF THE AUGMENTING ROW OF THE JACOBIAN.  THE FULL
C           ALGORITHM EMPLOYED HERE IS REALLY ONLY REQUIRED
C           WHEN IPC IS NOT EQUAL TO NVAR.  OTHERWISE, THE SOLUTION
C           IS QUITE EASY.
C
C  IPAR   - AN INTEGER ARRAY SIMILAR TO RPAR, EXCEPT THAT LOCATIONS
C           1 AND 2 OF IPAR ARE REFERENCED BY BANSLV.  IPAR(*)=ML,
C           THE LOWER BANDWIDTH OF THE JACOBIAN, AND IPAR(2)=MU,
C           THE UPPER BANDWIDTH.
C
C  IWORK  - AN INTEGER WORK ARRAY, USED BY THE CONTINUATION CODE
C           TO STORE STATISTICS, POINTERS, AND THE PIVOT VECTOR.
C           NOTE THAT IWORK(13) STORES THE FIRST ENTRY IN IWORK
C           DEVOTED TO THE PIVOT VECTOR, IWORK(15) STORES THE
C           FIRST ELEMENT OF RWORK DEVOTED TO THE JACOBIAN SUBSYSTEM.
C           IWORK(20) COUNTS MATRIX FACTORIZATIONS, AND IWORK(21)
C           COUNTS BACK-SOLVES.
C
C  LIW    - DIMENSION OF IWORK IN THE CALLING PROGRAM.
C
C  JOB    - WORK CONTROL SWITCH.
C           JOB=0, DECOMPOSE MATRIX, GET DETERMINANT, SOLVE SYSTEM.
C           JOB=1, SOLVE A NEW SYSTEM WITH SAME OLD MATRIX.
C           JOB=2, DECOMPOSE, COMPUTE DETERMINANT.
C           3, Check jacobian matrix.  Call user jacobian routine,
C           multiply by -1.0, add finite difference jacobian,
C           print largest entry.
C
C  NVAR   - THE NUMBER OF VARIABLES X, THE FORMAL DIMENSION OF
C           THE FULL LINEAR SYSTEM.
C
C  RWORK  - A REAL ARRAY USED FOR STORING THE MATRIX AND THE TWO
C           AUXILIARY VECTORS NEEDED.
C
C  LRW    - THE DIMENSION OF RWORK IN THE CALLING PROGRAM.
C
C  X      - A REAL VECTOR OF LENGTH NVAR, THE POINT AT WHICH
C           THE FULL LINEAR SYSTEM IS TO BE EVALUATED.
C
C  Y      - A REAL VECTOR OF LENGTH NVAR, WHICH ON INPUT IS
C           THE RIGHT HAND SIDE OF THE LINEAR SYSTEM TO BE SOLVED,
C           AND ON SUCCESSFUL RETURN WITH IERROR=0, IS THE
C           SOLUTION OF THAT LINEAR SYSTEM.
C
C
      REAL      ONE
      REAL      ZERO
C
      PARAMETER (ONE=1.0)
      PARAMETER (ZERO=0.0)
C
      EXTERNAL  BANJAC
      EXTERNAL  DF
      EXTERNAL  FX
      EXTERNAL  ISAMAX
      EXTERNAL  SAXPY
      EXTERNAL  SDOT
      EXTERNAL  SGBDI
      EXTERNAL  SGBFA
      EXTERNAL  SGBSL
      EXTERNAL  SSCAL
      EXTERNAL  SSWAP
C
      INTRINSIC MAX
      INTRINSIC MIN
      INTRINSIC SQRT
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      REAL      AK
      REAL      DET(2)
      REAL      DETS
      REAL      EPS
      REAL      FPAR(*)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(*)
      INTEGER   IPC
      INTEGER   IRL
      INTEGER   IRU
      INTEGER   ISAMAX
      INTEGER   ITEMP
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   JAC
      INTEGER   JACK
      INTEGER   JOB
      INTEGER   JTEMP
      INTEGER   K
      INTEGER   LDFL
      INTEGER   LDFX
      INTEGER   LDX
      INTEGER   LFXM
      INTEGER   LFXP
      INTEGER   LILST
      INTEGER   LOUNIT
      INTEGER   LPIV
      INTEGER   LRLST
      INTEGER   LROWIP
      INTEGER   MBAND
      INTEGER   ML
      INTEGER   MU
      INTEGER   NBAND
      INTEGER   NDIM
      INTEGER   NEQN
      INTEGER   NSWAP
      REAL      RWORK(LRW)
      REAL      SDOT
      REAL      SKALE
      REAL      TEMP
      REAL      X(NVAR)
      REAL      Y(NVAR)
C
C  CHECK INPUT DATA
C
      IERROR=0
      LOUNIT=IWORK(8)
      ML=IPAR(1)
      MU=IPAR(2)
      NEQN=NVAR-1
      IF(ML.LT.0.OR.ML.GT.NEQN)THEN
        IERROR=1
        WRITE(LOUNIT,1020)ML
        RETURN
        ENDIF
      IF(MU.LT.0.OR.MU.GT.NEQN)THEN
        IERROR=1
        WRITE(LOUNIT,1030)MU
        RETURN
        ENDIF
      MBAND=ML+MU+1
      NBAND=MBAND+ML
      LPIV=IWORK(13)
      LILST=LPIV+NEQN-1
      LDFX=IWORK(15)
      LDFL=LDFX+NBAND*NEQN
      LROWIP=LDFL+NVAR
      JAC=IWORK(9)
      IF(JAC.EQ.0.AND.JOB.NE.3)THEN
        LRLST=LROWIP+NVAR-1
      ELSE
        LFXP=LROWIP+NVAR
        LFXM=LFXP+NVAR
        LDX=LFXM+NVAR
        LRLST=LDX+NVAR-1
        ENDIF
      NDIM=NEQN*NBAND+NVAR+NVAR-1
      IF(LILST.GT.LIW.OR.LRLST.GT.LRW)THEN
        WRITE(LOUNIT,1040)LILST,LIW
        WRITE(LOUNIT,1050)LRLST,LRW
        IERROR=1
        RETURN
        ENDIF
      IF(JOB.EQ.1)THEN
        IF(IPC.EQ.NVAR)THEN
          GO TO 70
        ELSE
          GO TO 40
        ENDIF
      ENDIF
C
C  JOB=0 OR 2 MEANS WE MUST EVALUATE THE JACOBIAN, EITHER WITH THE USER'S
C  ROUTINE OR BY FINITE DIFFERENCES, FACTOR IT AND GET THE SIGN
C  OF THE DETERMINANT
C
      DO 10 I=1,NDIM
        RWORK(LDFX+I-1)=ZERO
10      CONTINUE
      IF(JAC.EQ.0)THEN
        CALL DF(NVAR,FPAR,IPAR,X,RWORK(LDFX),IERROR)
        IWORK(19)=IWORK(19)+1
        IF(IERROR.NE.0)RETURN
        RWORK(LROWIP-1+IPC)=ONE
        ENDIF
      IF(JOB.EQ.3)THEN
        SKALE=-ONE
        CALL SSCAL(NDIM,SKALE,RWORK(LDFX),1)
        ENDIF
      IF(JAC.EQ.1.OR.JAC.EQ.2.OR.JOB.EQ.3)THEN
        JACK=JAC
        IF(JOB.EQ.3)JACK=2
        EPS=SQRT(SQRT(RWORK(8)))
        CALL BANJAC(EPS,RWORK(LDFL),FPAR,RWORK(LDFX),RWORK(LROWIP),
     *  FX,IERROR,IPAR,IPC,IWORK,JACK,LIW,LOUNIT,NBAND,NEQN,NVAR,X,
     *  RWORK(LDX),RWORK(LFXP),RWORK(LFXM))
        ENDIF
C
C  Figure out formulas for I and J
C
      IF(JOB.EQ.3)THEN
        K=ISAMAX(NDIM,RWORK(LDFX),1)
        AK=RWORK(LDFX+K-1)
        IF(K.LE.NEQN*NBAND)THEN
          J=((K-1)/NBAND)+1
          I=K-(J-1)*NBAND+J-ML-MU-1
        ELSEIF(K.LE.NEQN*NBAND+NEQN)THEN
          I=K-NEQN*NBAND
          J=NVAR
        ELSE
          I=NVAR
          J=K-NEQN*NBAND-NEQN
          ENDIF
        WRITE(LOUNIT,1070)AK,I,J
        IF(IWORK(1).EQ.-2)THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'BANSLV - Entire difference matrix:'
          WRITE(LOUNIT,*)' '
          DO 99 I=1,NVAR
            DO 98 J=1,NVAR
              IF(J.EQ.NVAR)THEN
                K=LDFX-1+(NVAR-1)*NBAND+I
                WRITE(LOUNIT,1080)RWORK(K),I,J
              ELSEIF((J-I.LE.ML).AND.(I-J.LE.MU))THEN
                K=LDFX-1+I+J*(NBAND-1)-ML
                WRITE(LOUNIT,1080)RWORK(K),I,J
                ENDIF
98            CONTINUE
            WRITE(LOUNIT,*)' '
99          CONTINUE
          ENDIF
        RETURN
        ENDIF
      IF(IPC.EQ.NVAR)GO TO 60
C
C  SWITCH THE NVAR-TH AND IPC-TH ROWS.
C
      IRL=MAX(1,IPC-ML)
      IRU=MIN(NEQN,IPC+MU)
      NSWAP=IRU+1-IRL
      ITEMP=LDFX-ML-1+IPC+IRL*(NBAND-1)
      JTEMP=NBAND-1
      CALL SSWAP(NSWAP,RWORK(LROWIP-1+IRL),1,RWORK(ITEMP),JTEMP)
C
C  DECOMPOSE THE SUBMATRIX AND OBTAIN DETERMINANT SIGN
C
      CALL SGBFA(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),IERROR)
      IWORK(20)=IWORK(20)+1
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,1000)IERROR,IPC
        IERROR=3
        RETURN
        ENDIF
      CALL SGBDI(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),DET)
      DETS=ZERO
      IF(DET(1).GT.ZERO)THEN
        DETS=-ONE
      ELSEIF(DET(1).LT.ZERO)THEN
        DETS=ONE
        ENDIF
C
C  SET THE RIGHT HAND SIDE OF THE AUXILLIARY SYSTEM TO
C  THE LAST COLUMN OF THE JACOBIAN, MINUS E(IPC).  SHUFFLE THE
C  IPC-TH AND NVAR-TH ENTRIES OF THIS RIGHT HAND SIDE
C  TO REFLECT THE PIVOTING OF THE EQUATIONS, THEN SOLVE.
C
      TEMP=RWORK(LDFL+IPC-1)-ONE
      RWORK(LDFL+IPC-1)=RWORK(LDFL+NVAR-1)
      RWORK(LDFL+NVAR-1)=TEMP
      CALL SGBSL(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),RWORK(LDFL)
     * ,0)
      IWORK(21)=IWORK(21)+1
C
C  SOLVE FOR LAST ENTRY OF AUXILLIARY SOLUTION
C
      RWORK(LDFL+NVAR-1)=RWORK(LDFL+NVAR-1)
     * -SDOT(NEQN,RWORK(LROWIP),1,RWORK(LDFL),1)
C
C  ADJUST DETERMINANT SIGN
C
      IF(ONE+RWORK(LDFL+NVAR-1).EQ.ZERO)THEN
        IERROR=3
        WRITE(LOUNIT,*)'BANSLV - Algorithm fails, DENOM=0.0'
        RETURN
        ENDIF
      IF(ONE+RWORK(LDFL+NVAR-1).LT.ZERO)DETS=-DETS
      IF(JOB.EQ.2)RETURN
C
C  SOLVE SYSTEM
C
40    CONTINUE
      IF(IPC.EQ.NVAR)GO TO 70
C
C  MODIFY RIGHT HAND SIDE OF MAIN SYSTEM
C
      TEMP=Y(IPC)
      Y(IPC)=Y(NVAR)
      Y(NVAR)=TEMP
C
C  SOLVE SUBSYSTEM
C
      CALL SGBSL(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),Y,0)
      IWORK(21)=IWORK(21)+1
C
C  SOLVE FOR LAST ENTRY OF MAIN SOLUTION
C
      Y(NVAR)=Y(NVAR)-SDOT(NEQN,RWORK(LROWIP),1,Y,1)
C
C  CORRECT MAIN SOLUTION WITH MULTIPLE OF SUBSOLUTION
C
      SKALE=-Y(NVAR)/(ONE+RWORK(LDFL+NVAR-1))
      CALL SAXPY(NVAR,SKALE,RWORK(LDFL),1,Y,1)
      RETURN
C
C  FACTOR MATRIX FOR THE SPECIAL CASE IPC=NVAR
C
   60 CONTINUE
      CALL SGBFA(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),IERROR)
      IWORK(20)=IWORK(20)+1
      IF(IERROR.NE.0)THEN
        WRITE(LOUNIT,1000)IERROR,IPC
        IERROR=3
        RETURN
        ENDIF
      CALL SGBDI(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),DET)
      DETS=ZERO
      IF(DET(1).LT.ZERO)THEN
        DETS=-ONE
      ELSEIF(DET(1).GT.ZERO)THEN
        DETS=ONE
        ENDIF
      IF(JOB.EQ.2)RETURN
C
C  SOLVE SYSTEM, FOR SPECIAL CASE IPC=NVAR
C
   70 CONTINUE
      SKALE=-Y(NVAR)
      CALL SAXPY(NEQN,SKALE,RWORK(LDFL),1,Y,1)
      CALL SGBSL(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),Y,0)
      IWORK(21)=IWORK(21)+1
      RETURN
1000  FORMAT(' BANSLV - Zero pivot, SGBFA returns INFO=',I6,' IPC=',I6)
1020  FORMAT(' BANSLV - Illegal ML=',I6)
1030  FORMAT(' BANSLV - Illegal MU=',I6)
1040  FORMAT(' BANSLV - Need LIW=',I6,', have LIW=',I6)
1050  FORMAT(' BANSLV - Need LRW=',I6,', have LRW=',I6)
1070  FORMAT(' BANSLV - Maximum difference between user and estimated'/
     *       '          jacobian is ',G14.6,' row ',I6,' column ',I6)
1080  FORMAT(1X,G14.6,' =FP(I,J)-DELF(I,J), I, J=',2I6)
      END
      SUBROUTINE BANJAC(EPS,FCOL,FPAR,FPRIME,FROW,FX,IERROR,IPAR,IPC,
     * IWORK,JAC,LIW,LOUNIT,NBAND,NEQN,NVAR,X,XTEMP,WORK1,WORK2)
C
C**********************************************************************
C
C  BANJAC estimates the jacobian matrix FPRIME of the function FX, 
C  using forward or central finite differences.  BANJAC is called by 
C  BANSLV when the user has specified the jacobian option as 1 or 2.
C
C  EPS    Input, a tolerance used for shifting the X values.  A value
C         of the square root of the machine precision is usually
C         appropriate.
C
C  FCOL   Output, REAL FCOL(NEQN), the last column of the approximate 
C         jacobian, which is allowed to be "full".  This comprises
C         matrix entries FPRIME(1,NVAR) through FPRIME(NEQN,NVAR).
C
C  FPAR   Input/output, REAL FPAR(*), user parameter vector, not
C         touched by this routine, but passed on to user routines.
C
C  FPRIME Output, REAL FPRIME(NBAND,NEQN), is the array into which the
C         the banded portion of the computed jacobian will be stored.
C         The LINPACK general band format is used, assigning entry (I,J)
C         to FPRIME(I-J+ML+MU+1,J), where ML and MU are the lower and
C         upper half bandwidths respectively.
C
C  FROW   Output, REAL FROW(NVAR), storage for the last (augmenting) row
C         of the jacobian, which will be all zero except for a 1 in
C         location IPC.
C
C  FX     Input, EXTERNAL FX, the name of the routine which evaluates the 
C         function.
C
C         FX computes the value of the nonlinear function.  This name must be 
C         declared EXTERNAL in the calling program.  FX should evaluate the 
C         NVAR-1 function components at the input point X, and store the result 
C         in the vector FVEC.  An augmenting equation will be stored in entry 
C         NVAR of FVEC by the PITCON program.
C
C         FX should have the following form:
C
C           SUBROUTINE FX(NVAR,FPAR,IPAR,X,FVEC,IERROR)
C
C           NVAR   Input, INTEGER NVAR, number of variables.
C
C           FPAR   Input/output, REAL FPAR(*), array of user parameters.
C
C           IPAR   Input/output, INTEGER IPAR(*), array of user parameters.
C
C           X      Input, REAL X(NVAR), the point at which function evaluation
C                  is required.
C
C           FVEC   Output, REAL FVEC(NVAR), the value of the function at point 
C                  X.  Only the first NVAR-1 entries of FVEC are to be set by 
C                  the routine.  PITCON sets the final value itself.
C 
C           IERROR Output, INTEGER IERROR, error flag.  FX should set this to 0 
C                  if there are no problems, or to 2 if there is a problem.
C
C  IERROR Output, INTEGER IERROR, error flag.  A nonzero value means that
C         there was an error in the user routine FX, or in BANJAC itself.
C         In either case, the jacobian has not been computed.
C
C  IPAR   Input, INTEGER IPAR(*), a user parameter vector passed to FX.  
C         However, because this is a problem with a banded jacobian, entries
C         IPAR(1) and IPAR(2) are read by this routine.  IPAR(1) contains
C         ML, the lower half bandwidth of the jacobian, and IPAR(2) contains
C         MU, the upper half bandwidth of the jacobian.
C
C  IPC    Input, INTEGER IPC, the index of the current continuation parameter,
C         which is needed to determine the form of FROW.
C
C  IWORK  Input, INTEGER IWORK(LIW), work and statistics vector.  Only
C         required here so that we can count the number of function
C         evaluations.
C
C  JAC    Input, INTEGER JAC, the user requested jacobian option.  For
C         our purposes, the only two values of interest are:
C
C           1 = estimate jacobian with forward differences,
C           2 = estimate jacobian with central differences (twice the work)
C
C  LIW    Input, INTEGER LIW, the dimension of IWORK.
C
C  LOUNIT Input, INTEGER LOUNIT, the FORTRAN output unit for messages.
C
C  NBAND  Input, INTEGER NBAND, the first dimension of the jacobian matrix
C         FPRIME, NBAND=ML+MU+1.
C
C  NEQN   Input, INTEGER NEQN, the number of equations, equal to NVAR-1.
C
C  NVAR   Input, INTEGER NVAR, the number of variables.
C
C  X      Input, REAL X(NVAR), the point at which the jacobian is desired.
C
C  XTEMP,
C  WORK1,
C  WORK2  Work arrays, REAL XTEMP(NVAR), WORK1(NVAR), WORK2(NVAR).
C
      REAL      ONE
      REAL      TWO
C
      PARAMETER (ONE=1.0)
      PARAMETER (TWO=2.0)
C
      EXTERNAL  FX
      EXTERNAL  SAXPY
      EXTERNAL  SCOPY
      EXTERNAL  SSCAL
C
      INTRINSIC ABS
      INTRINSIC MAX
      INTRINSIC MIN
C
      INTEGER   LIW
      INTEGER   NBAND
      INTEGER   NEQN
      INTEGER   NVAR
C
      REAL      EPS
      REAL      FCOL(NEQN)
      REAL      FPAR(*)
      REAL      FPRIME(NBAND,NEQN)
      REAL      FROW(NVAR)
      INTEGER   IBAND
      INTEGER   IERROR
      INTEGER   IHI
      INTEGER   ILO
      INTEGER   IPAR(2)
      INTEGER   IPC
      INTEGER   IROW
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   JAC
      INTEGER   KCALL
      INTEGER   LOUNIT
      INTEGER   MBAND
      INTEGER   ML
      INTEGER   MU
      REAL      SKALE
      REAL      X(NVAR)
      REAL      XJAC
      REAL      XTEMP(NVAR)
      REAL      WORK1(NVAR)
      REAL      WORK2(NVAR)
C
      ML=IPAR(1)
      MU=IPAR(2)
      MBAND=ML+MU+1
      IF(JAC.EQ.1)THEN
        CALL FX(NVAR,FPAR,IPAR,X,WORK2,IERROR)
        IWORK(22)=IWORK(22)+1
        IF(IERROR.NE.0)RETURN
        ENDIF
      XJAC=ONE
      IF(JAC.EQ.2)XJAC=TWO
      DO 40 KCALL=1,MBAND
        CALL SCOPY(NVAR,X,1,XTEMP,1)
        DO 10 J=KCALL,NEQN,MBAND
          XTEMP(J)=X(J)+EPS*(ONE+ABS(X(J)))
10        CONTINUE
        CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK1,IERROR)
        IWORK(22)=IWORK(22)+1
        IF(IERROR.NE.0)RETURN
        IF(JAC.EQ.2)THEN
          CALL SCOPY(NVAR,X,1,XTEMP,1)
          DO 20 J=KCALL,NEQN,MBAND
            XTEMP(J)=X(J)-EPS*(ONE+ABS(X(J)))
20          CONTINUE
          CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK2,IERROR)
          IWORK(22)=IWORK(22)+1
          IF(IERROR.NE.0)RETURN
          ENDIF
        DO 30 J=KCALL,NEQN,MBAND
          ILO=MAX(1,J-MU)
          IHI=MIN(NEQN,J+ML)
          IROW=ILO-J+ML+MU+1
          IBAND=IHI-ILO+1
          SKALE=-ONE
          CALL SAXPY(IBAND,SKALE,WORK2(ILO),1,WORK1(ILO),1)
          SKALE=ONE/(XJAC*EPS*(ONE+ABS(X(J))))
          CALL SSCAL(IBAND,SKALE,WORK1(ILO),1)
          SKALE=ONE
          CALL SAXPY(IBAND,SKALE,WORK1(ILO),1,FPRIME(IROW,J),1)
30        CONTINUE
40      CONTINUE
C
C  Compute last column of jacobian, rows 1 to NEQN
C
      CALL SCOPY(NVAR,X,1,XTEMP,1)
      XTEMP(NVAR)=X(NVAR)+EPS*(ONE+ABS(X(NVAR)))
      CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK1,IERROR)
      IWORK(22)=IWORK(22)+1
      IF(IERROR.NE.0)RETURN
      IF(JAC.EQ.2)THEN
        XTEMP(NVAR)=X(NVAR)-EPS*(ONE+ABS(X(NVAR)))
        CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK2,IERROR)
        IWORK(22)=IWORK(22)+1
        IF(IERROR.NE.0)RETURN  
        ENDIF
      SKALE=-ONE
      CALL SAXPY(NEQN,SKALE,WORK2,1,WORK1,1)
      SKALE=ONE/(XJAC*EPS*(ONE+ABS(X(NVAR))))
      CALL SSCAL(NEQN,SKALE,WORK1,1)
      SKALE=ONE
      CALL SAXPY(NEQN,SKALE,WORK1,1,FCOL,1)
C
C  Do last row, J=1,NVAR
C
      FROW(IPC)=FROW(IPC)+ONE
      RETURN
      END
      FUNCTION REPS()
C
C********************************************************************
C
C  Compute the machine relative precision.
C
C  If you would prefer a simpler computation, you can replace this routine 
C  by a line that reads "REPS=1.0E-7" for example, or "REPS=R1MACH(4)"
C  if you have the PORT/SLATEC machine constant routines.  
C
      REAL      HALF
      REAL      ONE
      REAL      TWO
      REAL      ZERO
C
      PARAMETER (HALF=0.5)
      PARAMETER (ONE=1.0)
      PARAMETER (TWO=2.0)
      PARAMETER (ZERO=0.0)
C
      REAL      EPS
      INTEGER   I
      REAL      RADD
      REAL      RADDOLD
      REAL      REPS
      REAL      RMANT
      REAL      RMANTOLD
      REAL      TEMP1
      REAL      TEMP2
C
      RMANT=ONE
      RADD=ONE
      DO 10 I=1,100
        RADDOLD=RADD
        RADD=HALF*RADD
        RMANTOLD=RMANT
        RMANT=RMANT+RADD
        IF((RMANT.EQ.RMANTOLD).OR.
     *     (RMANT-RADD.NE.RMANTOLD))GO TO 20
        IF(RMANT.EQ.TWO)THEN
          RADDOLD=RADD
          GO TO 20
          ENDIF
10      CONTINUE
C
C  Loop doesn't terminate after 100 steps!
C  We'll just use that value.
C
20    CONTINUE
C
      RMANT=RMANTOLD
      EPS=RADDOLD
C
C  Additional requirement on the value of EPS is that
C  (1.0+EPS)-EPS.EQ.(1.0-EPS)+EPS.EQ.1.0)
C
30    CONTINUE
      TEMP1=ONE+EPS
      TEMP2=ONE-EPS
      IF(TEMP1-EPS.NE.ONE.OR.TEMP2+EPS.NE.ONE)THEN
        EPS=TWO*EPS
        GO TO 30
        ENDIF
C
      REPS=EPS
      RETURN
      END 
