      SUBROUTINE ALCON1 (FY,N,X,XW,TAU,TAUMIN,TAUMAX,EPS,INFO,
     &   RWORK,LRWORK,IWORK,LIWORK)
C*    Begin Prologue ALCON1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL FY
      DIMENSION RWORK(LRWORK),IWORK(LIWORK),X(N),XW(N),INFO(9)
C
C---------------------------------------------------------------------
C
C*  Title
C
C     (Al)gebraic system of equations (Con)tinuation method.
C
C*  Written by        P. Deuflhard, P. Kunkel
C*  Purpose           Solution of parameter dependent systems of 
C                     nonlinear equations.
C*  Method            Numerical pathfollowing with automatic steplength 
C                     control
C*  Category          F4 - Parameter Dependent Nonlinear Equation 
C                          Systems
C*  Keywords          Numerical pathfollowing, Homotopy Method
C*  Version           0.9
C*  Revision          September 1985
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C*  Copyright     (c) Konrad-Zuse-Zentrum fuer
C                     Informationstechnik Berlin (ZIB)
C                     Takustrasse 7, D-14195 Berlin-Dahlem
C                     phone : + 49/30/84185-0
C                     fax   : + 49/30/84185-125
C*  Contact           Bodo Erdmann
C                     ZIB, Division Scientific Computing, 
C                          Department Numerical Analysis and Modelling
C                     phone : + 49/30/84185-185
C                     fax   : + 49/30/84185-107
C                     e-mail: erdmann@zib.de
C
C*    Reference:
C
C     /1/ P. Deuflhard, B. Fiedler, P. Kunkel:
C         Efficient Numerical Pathfollowing Beyond Critical Points
C         Univ. Heidelberg, Sfb 123, Tech. Rep. 278 (1984)
C
C  ---------------------------------------------------------------
C
C* Licence
C    You may use or modify this code for your own non commercial
C    purposes for an unlimited time.
C    In any case you should not deliver this code without a special
C    permission of ZIB.
C    In case you intend to use the code commercially, we oblige you
C    to sign an according licence agreement with ZIB.
C
C* Warranty
C    This code has been tested up to a certain level. Defects and
C    weaknesses, which may be included in the code, do not establish
C    any warranties by ZIB. ZIB does not take over any liabilities
C    which may follow from aquisition or application of this code.
C
C* Software status
C    This code is under partial care of ZIB and belongs to ZIB
C    software class 2.
C
C     ------------------------------------------------------------
C
C*   Summary
C
C     Continuation method for systems of algebraic equations  f(x,tau)=0
C     Optional computation of turning and (simple) bifurcation points
C     ------------------------------------------------------------
C
C  INPUT PARAMETERS  (* MARKS INOUT PARAMETERS)
C
C    - FY (X,TAU,F)
C               EXTERNAL SUBROUTINE DEFINING SYSTEM OF
C               ALGEBRAIC EQUATIONS
C                 X(N)       VECTOR OF STATE VARIABLES    INPUT
C                 TAU        PARAMETER                    INPUT
C                 F(N)       VALUES OF F(X,TAU)           OUTPUT
C
C    - N        NUMBER OF ALGEBRAIC EQUATIONS
C
C    - X(N)   * ESTIMATE OF SOLUTION ON CONTINUATION PATH
C               FOR INITIAL TAU
C
C    - XW(N)    INITIAL SCALING QUANTITIES FOR X(N)
C               (ALCON1 USES ADAPTIVE INTERNAL SCALING
C               ALONG THE CONTINUATION PATH)
C
C    - TAU    * INITIAL VALUE OF PARAMETER TAU
C
C    - TAUMIN   MINIMUM ALLOWED VALUE FOR TAU
C
C    - TAUMAX   MAXIMUM ALLOWED VALUE FOR TAU
C
C    - EPS      REQUIRED RELATIVE ACCURACY
C
C    - INFO(9)  ARRAY USED FOR COMMUNICATION BETWEEN
C               THE PARTICULAR ROUTINES
C                 INFO(1)  PRINT PARAMETER
C                            0  NO PRINT
C                            1  PRINT OF X(N), TAU, AND INFORMATION
C                               ILLUSTRATING THE CONTINUATION PROCESS
C                            2  ADDITIONALLY INFORMATION ABOUT THE
C                               CONVERGENCE BEHAVIOR OF THE
C                               GAUSS-NEWTON METHOD
C                 INFO(2)  MAXIMUM PERMITTED NUMBER OF
C                          CONTINUATION STEPS
C                 INFO(3)  INITIAL DIRECTION
C                           +1  PATHFOLLOWING ALONG A SINGLE BRANCH
C                               STARTING IN POSITIVE TAU-DIRECTION
C                           -1  PATHFOLLOWING ALONG A SINGLE BRANCH
C                               STARTING IN NEGATIVE TAU-DIRECTION
C                 INFO(4)  CRITICAL POINT OPTIONS
C                            0  NO CRITICAL POINT DETERMINATION
C                            1  TURNING POINT COMPUTATION
C                            2  BIFURCATION POINT COMPUTATION
C                            3  TURNING AND BIFURCATION POINT
C                               COMPUTATION  (RECOMMENDED OPTION)
C                 INFO(5)  ERROR EXIT PARAMETER  (SEE OUTPUT PARAMETERS)
C                 INFO(6)  INTERNALLY USED
C                 INFO(7)  INTERNALLY USED
C                 INFO(8)  INTERNALLY USED
C                 INFO(9)  INTERNALLY USED
C
C    - RWORK    REAL WORKSPACE
C
C    - LRWORK   LENGTH OF REAL WORKSPACE
C               MUST BE AT LEAST   3 * N**2 + 25 * N + 19
C
C    - IWORK    INTEGER WORKSPACE
C
C    - LIWORK   LENGTH OF INTEGER WORKSPACE
C               MUST BE AT LEAST   N + 1
C
C    REMARK: FOR EASE OF IMPLEMENTATION ONLY REGULAR SOLUTIONS ARE
C            ASSUMED NEAR INITIAL TAU, TAUMIN, AND TAUMAX. COMPUTATIONS
C            ARE PERFORMED STRICTLY BETWEEN TAUMIN AND TAUMAX.
C
C
C  OUTPUT PARAMETERS
C
C    - X(N)     FINAL SOLUTION VALUES
C
C    - TAU      FINAL PARAMETER VALUE
C
C    - INFO(5)  ERROR EXIT PARAMETER
C                 0  NO ERROR OCCURRED
C                 1  MORE THAN INFO(2) CONTINUATION STEPS
C                 2  STEPLENGTH IN CONTINUATION PROCESS TOO SMALL
C                    (RELATIVE DIFFERENCE LESS THAN 10*EPS)
C                 3  RANK DEFICIENT JACOBIAN
C                 4  INITIAL GUESS OF X(N) TOO BAD
C                 5  NOT USED IN ALCON1
C                 6  WORKSPACE TOO SMALL
C
C
C  EXTERNAL UNITS
C
C      COMMON /A1UNIT/ UPR,UDIAG
C
C      - UPR     PRINT UNIT
C                (STANDARD UNIT 6)
C      - UDIAG   PLOT INFORMATION UNIT, INPUT TO PLOT ROUTINE PLTHM
C                DISC FILE WITH CARD IMAGE  TO BE DECLARED BY THE USER
C                SET UDIAG=0, IF NO PLOT IS DESIRED
C                (STANDARD UNIT 2)
C
C     ------------------------------------------------------------
C*    End Prologue
      INTEGER UPR,UDIAG
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A1UNIT/ UPR,UDIAG
C  MACHINE DEPENDENT CONSTANTS
C            EPMACH  UNIT ROUND OFF
C            SMALL   SQRT OF SMALLEST MACHINE NUMBER
      CALL ZIBCONST(EPMACH,SMALL)
C  WORKSPACE REQUIREMENTS
      N1=N+1
      NN1=N*N1
      N1Q=N1*N1
      LRWREQ=3*N**2+25*N+19
      LIWREQ=N+1
      IF (LRWREQ.GT.LRWORK .OR. LIWREQ.GT.LIWORK) GOTO 9000
C  REAL WORKSPACE DISTRIBUTION
      LY=1
      LYQ=LY+N1
      LYD=LYQ+N1
      LYDH=LYD+N1
      LYH=LYDH+N1
      LDYT=LYH+N1
      LDYTA=LDYT+N1
      LDYTH=LDYTA+N1
      LV=LDYTH+N1
      LD=LV+N1
      LYW=LD+N1
      LETA=LYW+N1
      LDFY=LETA+N1
      LA=LDFY+NN1
      LAH=LA+NN1
      LYA=LAH+N1Q
      LDY=LYA+N1
      LDYQ=LDY+N1
      LW=LDYQ+N1
      LF=LW+N1
      LFH=LF+N
      LU=LFH+N
      LW1=LU+N
      LW2=LW1+N1
C  INTEGER WORKSPACE DISTRIBUTION
      LIPIV=1
C  REORGANIZE INITIAL VALUE
      DO 1000 I=1,N
      RWORK(I)=X(I)
1000  CONTINUE
      RWORK(N1)=TAU
C  EXTERNAL INITIAL SCALING
      DO 2000 I=1,N
      RWORK(LYW+I-1)=XW(I)
2000  CONTINUE
      RWORK(LYW+N)=DMAX1(DABS(TAUMIN),DABS(TAUMAX))
C  CALL CONTINUATION ROUTINE
      CALL A1HOMQ(FY,N,N1,TAUMIN,TAUMAX,EPS,INFO,
     &   RWORK(LY),RWORK(LYQ),RWORK(LYD),RWORK(LYDH),RWORK(LYH),
     &   RWORK(LDYT),RWORK(LDYTA),RWORK(LDYTH),
     &   RWORK(LV),RWORK(LD),RWORK(LYW),RWORK(LETA),
     &   RWORK(LDFY),RWORK(LA),RWORK(LAH),
     &   RWORK(LYA),RWORK(LDY),RWORK(LDYQ),RWORK(LW),
     &   RWORK(LF),RWORK(LFH),RWORK(LU),
     &   RWORK(LW1),RWORK(LW2),IWORK(LIPIV))
C  RESTORE FINAL VALUES
      DO 3000 I=1,N
      X(I)=RWORK(I)
3000  CONTINUE
      TAU=RWORK(N1)
      RETURN
C  FAIL EXIT
9000  CONTINUE
      INFO(5)=6
      KPRINT=INFO(1)
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0 .AND. LRWREQ.GT.LRWORK) WRITE(UPR,60002) LRWREQ
      IF (KPRINT.GT.0 .AND. LIWREQ.GT.LIWORK) WRITE(UPR,60003) LIWREQ
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      RETURN
60001 FORMAT(///)
60002 FORMAT('   REAL  WORKSPACE NOT SUFFICIENT',4X,'AT LEAST',I5,
     &   ' LOCATIONS NEEDED')
60003 FORMAT(' INTEGER WORKSPACE NOT SUFFICIENT',4X,'AT LEAST',I5,
     &   ' LOCATIONS NEEDED')
      END
C
C
      SUBROUTINE A1HOMQ (FY,N,N1,TAUMIN,TAUMAX,EPS,INFO,Y,YQ,YD,YDH,YH,
     &   DYT,DYTA,DYTH,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,
     &   W1,W2,IPIV)
C
C  SUBROUTINE  A1HOMQ  TO BE USED WITH ROUTINE ALCON1
C
C  DRIVER ROUTINE FOR ALCON1
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),YQ(N1),YD(N1),YDH(N1),YH(N1)
      DOUBLE PRECISION DYT(N1),DYTA(N1),DYTH(N1)
      DOUBLE PRECISION V(N1),D(N1),YW(N1),ETA(N1)
      DOUBLE PRECISION DFY(N,N1),A(N,N1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(N1),DYQ(N1),W(N1)
      DOUBLE PRECISION F(N),FH(N),U(N)
      DOUBLE PRECISION W1(N1),W2(N1)
      INTEGER INFO(9),IPIV(N1)
      INTEGER UPR,UDIAG
      EXTERNAL FY
      COMMON /COUNT/ IFCTEV,ITER,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A1UNIT/ UPR,UDIAG
      COMMON /A1VALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      COMMON /A1METH/ FACTOR,IEXTR,IJACM
C-----------------------------------------------------------------------
C  PREPARATIONS
      EPDIFF=1.D3*EPMACH
      ETADIF=DSQRT(EPDIFF)
      ITER=0
      IHALT=0
      IRED=0
      KPRINT=INFO(1)
      DO 1010 I=1,N1
      ETA(I)=ETADIF
1010  CONTINUE
      IFCTEV=0
      IDECS=0
      ISOLS=0
      INFO(6)=IGNMAX
C  CHECK FOR CONSISTENCY OF INITIAL VALUES
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60201)
      IF (KPRINT.GT.2) WRITE(UPR,60500) Y(N1),(Y(I),I=1,N)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INFO(8)=1
      CALL A1GNHM(N,N1,Y,FY,EPS,INFO,SUMYH,SUMYQH,TAUMIN,TAUMAX,
     &     IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U)
      IERR=INFO(9)
      IF (IERR.NE.0) GOTO 9940
      ITS=INFO(7)
      CALL A1TAND(N,N1,FY,Y,F,FH,YW,ETA,DFY,A,COND,DETFA,IPIV,AH,V,
     &   D,IRANK)
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) GOTO 9940
      DETHMA=DETFA
      IF (IPIVS.NE.N1) CALL A1DETH(N,N1,A,AH,D,N1,IPIV,DETHMA)
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
      IF (KPRINT.GT.0) WRITE(UPR,60500) Y(N1),(Y(I),I=1,N)
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
C.......................................................................
C  START OF BRANCH FOLLOWING
      INIT=0
      ITER=1
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60202) ITER
      DO 1110 I=1,N1
      YQ(I)=Y(I)
1110  CONTINUE
C.......................................................................
C  INITIAL ESTIMATE FOR STEPLENGTH
      DO 1210 I=1,N
      DYT(IPIV(I))=-V(I)/D(N1)
1210  CONTINUE
      DYT(IPIVS)=1.D0/D(N1)
      SIGMA=SMALL1
      SIGNUM=DBLE(INFO(3))
      IF (DYT(N1)*SIGNUM.GT.0.D0) GOTO 1221
      DO 1220 I=1,N1
      DYT(I)=-DYT(I)
1220  CONTINUE
1221  CONTINUE
      IPIVA=IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20102) IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) Y,DYT,YW
      GOTO 2200
C-----------------------------------------------------------------------
C  STEPLENGTH PREDICTOR
2000  CONTINUE
      ITER=ITER+1
      IF (ITER.GT.INFO(2)) GOTO 9910
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60202) ITER
      IF (ITS.EQ.0) GOTO 2020
      TH=DSQRT(SUMYQH/SUMYH)
      VAL1=0.5D0*THMAX*DSQRT(SUMYH)/DABS(VAL1)
      VAL2=0.D0
      DO 2010 I=1,N1
      VAL2=VAL2+((YDH(I)-YQ(I))/YW(I))**2
2010  CONTINUE
      VAL2=TH*DSQRT(VAL2)
      IF (VAL2.LT.THR*VAL1) VAL2=THR*VAL1
      R=DSQRT(VAL1/VAL2)
      SIGMA=R*SIGMA
      GOTO 2030
C  EMPIRICAL STEPLENGTH INCREASE IN NEARLY LINEAR CASE
2020  CONTINUE
      SIGMA=SIGMA/DSQRT(THR)
C  EXTRAPOLATED STEPLENGTH BOUND
2030  CONTINUE
      IF (IHALT.EQ.0) GOTO 2040
      SIGMAQ=FACTOR*(YHALT-YQ(IPIVA))/(DYT(IPIVA)*YW(IPIVA))
      IF (SIGMAQ.LE.10.D0*EPS .OR. SIGMA.LE.SIGMAQ) GOTO 2040
      SIGMA=SIGMAQ
      IF (KPRINT.GT.1) WRITE(UPR,60206)
2040  CONTINUE
      IF (INIT.NE.0) IPIVA=IPIVS
      GOTO 2200
C.......................................................................
C  STEPLENGTH CORRECTOR
2110  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60401)
      TH=DSQRT(SUMYQH/SUMYH)
      R=DSQRT(0.5D0*THMAX/TH)
      IF (R.GT.0.7D0) R=0.7D0
      IF (R.LT.0.1D0) R=0.1D0
      GOTO 2190
2120  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60402)
      R=0.7D0
      GOTO 2190
2130  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60403)
      R=0.7D0
2190  CONTINUE
      IF (SIGMA.LT.EPS) GOTO 9920
      INIT=1
      SIGMA=R*SIGMA
C.......................................................................
C  TRIAL VALUE
2200  CONTINUE
      SIGMIN=0.D0
      SIGMAX=0.D0
      IF (DYT(N1).EQ.0.D0) GOTO 2211
      SIGMAX=(TAUMAX-YQ(N1))/(DYT(N1)*YW(N1))
      IF (SIGMAX.GT.0.D0 .AND. SIGMA.GT.SIGMAX) SIGMA=SIGMAX
      SIGMIN=(TAUMIN-YQ(N1))/(DYT(N1)*YW(N1))
      IF (SIGMIN.GT.0.D0 .AND. SIGMA.GT.SIGMIN) SIGMA=SIGMIN
2211  CONTINUE
      DO 2220 I=1,N1
      YD(I)=YQ(I)+SIGMA*DYT(I)*YW(I)
      YDH(I)=YD(I)
2220  CONTINUE
C.......................................................................
C  RESCALING
      IF (IERR.NE.0) GOTO 2490
      DO 2410 I=1,N1
      T=0.5D0*(DABS(Y(I))+DABS(YQ(I)))
      IF (T.LT.YW(I)) T=YW(I)
      DYTH(I)=DYT(I)*YW(I)/T
      YW(I)=T
2410  CONTINUE
      DO 2420 I=1,N1
      Y(I)=YQ(I)
2420  CONTINUE
2490  CONTINUE
      IF (SIGMA.EQ.SIGMIN .OR. SIGMA.EQ.SIGMAX) GOTO 3100
C-----------------------------------------------------------------------
C  ITERATION BACK TO CONTINUATION PATH BY GAUSS-NEWTON METHOD
      IF (KPRINT.GT.0) WRITE(UPR,60205) SIGMA,IPIVA
      IF (KPRINT.GT.2) WRITE(UPR,60500) YD(N1),(YD(I),I=1,N)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INFO(8)=2
      CALL A1GNHM(N,N1,YD,FY,EPS,INFO,SUMYH,SUMYQH,TAUMIN,TAUMAX,
     &     IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U)
      IERR=INFO(9)
      ITS=INFO(7)
      IF (IERR.NE.0 .OR. ITS.NE.0) INIT=1
      IF (INIT.EQ.0) GOTO 2020
      GOTO (2110,2110,9930,3110),IERR
      VAL1=0.D0
      DO 3010 I=1,N1
      VAL1=VAL1+DYT(I)*V(I)
3010  CONTINUE
      CALL A1TAND(N,N1,FY,YD,F,FH,YW,ETA,DFY,A,COND,DETF,IPIV,AH,V,
     &   D,IRANK)
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) GOTO 9930
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
      GOTO 4000
C.......................................................................
C  HIT FINAL VALUE
3100  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60205) SIGMA,IPIVA
      IF (SIGMA.EQ.SIGMIN) YD(N1)=TAUMIN
      IF (SIGMA.EQ.SIGMAX) YD(N1)=TAUMAX
      GOTO 3120
3110  CONTINUE
      IF (YD(N1).LT.TAUMIN) YD(N1)=TAUMIN
      IF (YD(N1).GT.TAUMAX) YD(N1)=TAUMAX
3120  CONTINUE
      IF (KPRINT.GT.2) WRITE(UPR,60500) YD(N1),(YD(I),I=1,N)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INIT=2
      INFO(8)=3
      CALL A1GNHM(N,N1,YD,FY,EPS,INFO,SUMYH,SUMYQH,TAUMIN,TAUMAX,
     &     IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U)
      IERR=INFO(9)
      ITS=INFO(7)
      GOTO (2110,2110,9930),IERR
      CALL A1TAND(N,N1,FY,YD,F,FH,YW,ETA,DFY,A,COND,DETF,IPIV,AH,V,
     &   D,IRANK)
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) GOTO 9930
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
C-----------------------------------------------------------------------
C  DETERMINANTS
4000  CONTINUE
      DETHM=DETF
      IF (IPIVS.NE.N1) CALL A1DETH(N,N1,A,AH,D,N1,IPIV,DETHM)
      DETFH=DETF
      IF (IPIVS.NE.IPIVA) CALL A1DETH(N,N1,A,AH,D,IPIVA,IPIV,DETFH)
      IF (KPRINT.GT.1) WRITE(UPR,60501) DETFA,DETFH,DETHMA,DETHM
C.......................................................................
C  NEW NORMALIZED TANGENT
      DO 4110 I=1,N
      W(IPIV(I))=-V(I)/D(N1)
4110  CONTINUE
      W(IPIVS)=1.D0/D(N1)
C  CHECK SIGN OF TANGENT
      SIGNUM=DSIGN(1.D0,YD(IPIVS)-Y(IPIVS))
      IF (W(IPIVS)*SIGNUM.GT.0.D0) GOTO 4121
      DO 4120 I=1,N1
      W(I)=-W(I)
4120  CONTINUE
4121  CONTINUE
C.......................................................................
C  CHECK FOR CONFLICT IN SIGNS
      IF (SIGMA.LE.10.D0*EPS) GOTO 4210
      IF (DYT(N1)*W(N1).LE.0.D0 .AND.
     &   (DETHMA*DETHM.GT.0.D0 .OR. IPIVS+IPIVA.EQ.2*N1)) GOTO 2120
      IF (DYT(N1)*W(N1).GT.0.D0 .AND. DETHMA*DETHM.LT.0.D0 .AND.
     &   DETFH*DETFA.GT.0.D0) GOTO 2120
      IF (DETFH*DETFA.LE.0.D0 .AND. DETHM*DETHMA.GT.0.D0) GOTO 2120
      IF (DYT(N1)*W(N1).LE.0.D0 .AND.
     &   DETHMA*DETHM.LE.0.D0 .AND. DETFH*DETFA.LE.0.D0) GOTO 2120
      IF (DYT(IPIVS)*W(IPIVS).LT.0.D0) GOTO 2120
      IF (DYT(IPIVA)*W(IPIVA).LT.0.D0) GOTO 2120
4210  CONTINUE
C.......................................................................
C  SAVE VALUES FOR NEXT ITERATE
      DO 4310 I=1,N1
      DYTA(I)=DYT(I)
      DYT(I)=W(I)
      YQ(I)=YD(I)
4310  CONTINUE
C-----------------------------------------------------------------------
C  CHECK FOR TURNING POINT
      IF (DYTA(N1)*DYT(N1).GT.0.D0) GOTO 5010
      CALL A1TURN(N,N1,IPIVA,EPS,INFO,Y,W1,YD,YH,DYT,DYTH,W2,V,YW,
     &   IPIV,FY,TAUMIN,TAUMAX,D,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U)
      IERRQ=INFO(9)
      IF (IERRQ.GT.0 .AND. SIGMA.GE.10.D0*EPS) GOTO 5200
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20103) IPIV(N1)
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20000) YH,W,YW
5010  CONTINUE
C.......................................................................
C  CHECK FOR BIFURCATION POINT
      IF (DETFH*DETFA.GT.0.D0) GOTO 6000
      DO 5111 I=1,N1
      YD(I)=YQ(I)
5111  CONTINUE
      CALL A1BIFU(N,N1,IPIVA,EPS,INFO,Y,W1,YD,YH,V,YW,IPIV,FY,
     &   TAUMIN,TAUMAX,D,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,DETFH,DETFA)
      IERRQ=INFO(9)
      IF (IERRQ.GT.0) GOTO 5200
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20101) IPIV(N1)
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20000) YH,W,YW
      GOTO 6000
C.......................................................................
C  RESTORE FORMER VALUES
5200  CONTINUE
      DO 5210 I=1,N1
      YQ(I)=Y(I)
      DYT(I)=DYTA(I)
5210  CONTINUE
      INIT=1
      GOTO 2130
C-----------------------------------------------------------------------
C  FINISH OUTPUT FOR CURRENT CONTINUATION STEP
6000  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60500) YQ(N1),(YQ(I),I=1,N)
      IF (UDIAG.GT.0) WRITE(UDIAG,20100) IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) YQ,DYT,YW
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INIT.EQ.2) GOTO 9900
C  EXTRAPOLATED STEPLENGTH BOUND
      IHALT=0
      IF (DETFA*DETFH.LE.0.D0 .OR. DABS(DETFA).LE.DABS(DETFH) .OR.
     &   IEXTR.NE.1) GOTO 6020
      IHALT=1
      S=DETFH/DETFA
      YHALT=(YQ(IPIVA)-S*Y(IPIVA))/(1.D0-S)
6020  CONTINUE
      DETFA=DETF
      DETHMA=DETHM
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9900  CONTINUE
      INFO(5)=0
      IF (KPRINT.GT.0) WRITE(UPR,60900) IFCTEV
      DO 9901 I=1,N1
      Y(I)=YQ(I)
9901  CONTINUE
      RETURN
C  FAIL EXIT
C  MORE THAN ITMAX STEPS
9910  CONTINUE
      INFO(5)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      RETURN
C  STEPLENGTH TOO SMALL
9920  CONTINUE
      INFO(5)=2
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      RETURN
C  RANK DEFICIENT SOLUTION OBTAINED
9930  CONTINUE
      INFO(5)=3
      IF (KPRINT.GT.0) WRITE(UPR,60500) YD(N1),(YD(I),I=1,N)
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60903)
      RETURN
C  STARTING VALUES TOO BAD
9940  CONTINUE
      INFO(5)=4
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60904)
      RETURN
C-----------------------------------------------------------------------
20000 FORMAT(4D18.10)
20100 FORMAT('0',I4)
20101 FORMAT('1',I4)
20102 FORMAT('2',I4)
20103 FORMAT('3',I4)
60001 FORMAT(/)
60002 FORMAT(///)
60100 FORMAT(' ',132('*'))
60201 FORMAT(' STEP  0       CHECK FOR CONSISTENCY OF INITIAL VALUES'/)
60202 FORMAT(' STEP',I4/)
60205 FORMAT(' SCALED STEPLENGTH',D11.4/
     &   ' CURRENT CONTINUATION PARAMETER',I4/)
60206 FORMAT(' STEPLENGTH REDUCED BY EXTRAPOLATED STEPLENGTH BOUND')
60300 FORMAT(/' GAUSS-NEWTON METHOD REQUIRED',I3,' ITERATIONS'//)
60401 FORMAT(/' GAUSS-NEWTON METHOD FAILED'/' ',132('.')/)
60402 FORMAT(/' CONFLICT IN SIGNS'/' ',132('.')/)
60403 FORMAT(/' STEP RETRIED WITH REDUCED STEPLENGTH'/' ',132('.')/)
60500 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60501 FORMAT(/' DET WRT ACTUAL CON PARM',2D17.7/
     &   ' DET WRT  REAL  CON PARM',2D17.7/)
60900 FORMAT(' ALCON1 REQUIRED',I5,' FUNCTION EVALUATIONS'//)
60901 FORMAT(' TERMINATION AFTER INFO(2) CONTINUATION STEPS'//)
60902 FORMAT(' TERMINATION SINCE STEPLENGTH TOO SMALL'//)
60903 FORMAT(' TERMINATION SINCE JACOBIAN RANK DEFICIENT'//)
60904 FORMAT(' TERMINATION SINCE INITIAL GUESS TOO FAR AWAY FROM',
     &   ' A REGULAR POINT OF THE CONTINUATION PATH'/
     &   ' USE NONLINEAR EQUATION SOLVER FOR BETTER INITIAL DATA'//)
      END
C
C
      SUBROUTINE A1GNHM (N,N1,Y,FY,EPS,INFO,SUMYH,SUMYQH,TAUMIN,TAUMAX,
     &   PIVOT,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U)
C
C  SUBROUTINE  A1GNHM  TO BE USED WITH ROUTINE ALCON1
C
C  GAUSS-NEWTON METHOD AS CORRECTOR IN PATHFOLLOWING PROCEDURE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),V(N1),D(N1),YW(N1),ETA(N1)
      DOUBLE PRECISION DFY(N,N1),A(N,N1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(N1),DYQ(N1),W(N1)
      DOUBLE PRECISION F(N),FH(N),U(N)
      INTEGER PIVOT(N1),INFO(9)
      INTEGER UPR,UDIAG
      EXTERNAL FY
      COMMON /COUNT/ IFCTEV,IDUMQ,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A1UNIT/ UPR,UDIAG
      COMMON /A1VALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      COMMON /A1METH/ FACTOR,IEXTR,IJACM
C  INTERNAL PARAMETERS
      SIGMA=0.5D0/THMAX
C  INITIAL PREPARATIONS TO START A1GNHM
      THMAX2=1.D0
      IF (INFO(8).NE.1) THMAX2=THMAX**2
      KPRINT=INFO(1)
      ITMAX=INFO(6)
      ITER=0
      LEVEL=0
      NEW=0
      NQ=N
      IF (INFO(8).EQ.2) NQ=N1
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y
C-----------------------------------------------------------------------
C  COMPUTATION OF RESIDUAL VECTOR
1     CALL FY(Y,Y(N1),F)
      IFCTEV=IFCTEV+1
      IF (LEVEL.EQ.1 .AND. IJACM.EQ.3 .AND. NEW.GT.0) GOTO 4
      IF (LEVEL.EQ.1) GOTO 43
C-----------------------------------------------------------------------
C  DIFFERENCE APPROXIMATION OF THE SCALED JACOBIAN MATRIX A
C  INCLUDING FEED-BACK DEVICE
2     NEW=0
      CALL A1DERI(FY,N,N1,Y,F,FH,YW,ETA,DFY)
      GOTO 4
C-----------------------------------------------------------------------
C  RANK-1 APPROXIMATION OF THE JACOBIAN DUE TO BROYDEN
C  (SCALED VERSION)
3     NEW=NEW+1
301   ST=1.D0/SUMY
      DO 31 I=1,N
      T=F(I)*ST
      DO 311 K=1,NQ
311   IF (DFY(I,K).NE.0.D0) DFY(I,K)=DFY(I,K)+T*DY(K)
31    CONTINUE
      IF (IJACM.EQ.3) GOTO 434
C-----------------------------------------------------------------------
C  SOLUTION OF THE LINEAR SYSTEM
4     DO 41 K=1,NQ
      DO 411 I=1,N
411   A(I,K)=-DFY(I,K)
41    CONTINUE
C  HOUSEHOLDER TRIANGULARIZATION
      COND=1.D0/EPMACH
      IRANK=N
      CALL A1DECC(A,N,N1,0,N,NQ,IRANK,COND,D,PIVOT,0,AH,W)
      IDECS=IDECS+1
      D1=DABS(D(1))
      IF (IRANK.LT.N) GOTO 93
      IF (ITER.NE.0 .OR. LEVEL.EQ.1 .OR. INFO(8).NE.2) GOTO 42
      DO 421 I=1,N
421   V(PIVOT(I))=-W(I)/D(N1)
      V(PIVOT(N1))=1.D0/D(N1)
42    CONTINUE
C  LINEAR LEAST SQUARES SOLUTION
43    DO 431 K=1,N
431   U(K)=F(K)
      CALL A1SOLC(A,N,N1,0,N,NQ,DYQ,U,IRANK,D,PIVOT,0,AH,W)
      ISOLS=ISOLS+1
      IF (LEVEL.EQ.1) GOTO 44
432   CONV=0.D0
      DO 433 L=1,NQ
      T=DYQ(L)
      CONV=DMAX1(CONV,DABS(T))
      YA(L)=Y(L)
433   Y(L)=Y(L)+T*YW(L)
C  TEST OF ACCURACY
      IF (CONV.LE.EPS) GOTO 9
      IF (IJACM.EQ.3 .AND. NEW.GT.0) GOTO 301
434   SUMY=0.D0
      DO 435 L=1,NQ
      T=DYQ(L)
      SUMY=SUMY+T*T
435   DY(L)=T
      LEVEL=1
      GOTO 1
44    SUMYQ=0.D0
      DO 441 L=1,NQ
      T=DYQ(L)
441   SUMYQ=SUMYQ+T*T
C-----------------------------------------------------------------------
C  RESTRICTED NATURAL MONOTONICITY TEST
      SUMYQ=DMAX1(SMALL*SUMY,SUMYQ)
      IF (ITER.NE.0) GOTO 51
      SUMYH=SUMY
      SUMYQH=SUMYQ
51    CONTINUE
      IF (SUMYQ.GT.THMAX2*SUMY) GOTO 69
C  A-POSTERIORI ESTIMATE OF RELAXATION FACTOR
      FCH=0.5D0*DSQRT(SUMY/SUMYQ)
C-----------------------------------------------------------------------
C  PREPARATIONS TO START THE FOLLOWING ITERATION STEP
      LEVEL=0
      ITER=ITER+1
      IF (KPRINT.GT.1 .AND. ITER.EQ.1) WRITE(UPR,60003)
      IF (KPRINT.GT.1) WRITE(UPR,60002)
     &   ITER,SUMY,SUMYQ,NEW,IRANK,COND,D1,(PIVOT(I),I=1,NQ)
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y
      IF (Y(N1).LT.TAUMIN .OR. Y(N1).GT.TAUMAX) GOTO 94
      IF (ITER.GE.ITMAX) GOTO 91
      IF (FCH.LT.SIGMA .OR. IRANK.LT.N .OR. IJACM.EQ.2) GOTO 2
      IF (IJACM.EQ.1) GOTO 3
      IF (IJACM.NE.3) GOTO 63
      ALPHA=0.D0
      DO 61 L=1,NQ
61    ALPHA=ALPHA+DY(L)*DYQ(L)
      ALPHA=ALPHA/SUMY
      T=1.D0/(1.D0-ALPHA)
      DO 62 L=1,NQ
62    DYQ(L)=T*DYQ(L)
      NEW=NEW+1
      GOTO 432
63    NEW=NEW-1
      GOTO 43
69    CONTINUE
      IF (INFO(8).NE.1 .OR. NEW.EQ.0) GOTO 92
      DO 691 I=1,NQ
691   Y(I)=YA(I)
      LEVEL=0
      GOTO 1
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9     INFO(7)=ITER
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  TERMINATION AFTER ITMAX ITERATIONS
91    INFO(9)=1
      IF (KPRINT.GT.1) WRITE(UPR,60901)
      RETURN
C  RESTRICTED MONOTONICITY TEST VIOLATED
92    INFO(9)=2
      IF (KPRINT.GT.1) WRITE(UPR,60902)
      SUMYH=SUMY
      SUMYQH=SUMYQ
      RETURN
C  RANK DEFICIENCY OF JACOBIAN
93    INFO(9)=3
      IF (KPRINT.GT.1) WRITE(UPR,60903)
      RETURN
C  TAUMIN OR TAUMAX CROSSED
94    INFO(9)=4
      IF (KPRINT.GT.1) WRITE(UPR,60904)
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60002 FORMAT(5X,I2,4X,D10.4,3X,D10.4,4X,I2,4X,I2,
     &   4X,D8.2,3X,D8.2,3X,20I3)
60003 FORMAT(4X,'ITER',5X,'LEVELX',6X,'LEVELXQ',6X,'NEW',2X,'RANK',
     &   5X,'COND',7X,'SENS',9X,'PIVOTS')
60901 FORMAT(' TERMINATION AFTER ITMAX ITERATIONS')
60902 FORMAT(' TERMINATION SINCE RESTRICTED MONOTONICITY TEST WAS ',
     &   'VIOLATED')
60903 FORMAT(' RANK DEFICIENT JACOBIAN')
60904 FORMAT(' TERMINATION SINCE TAUMIN OR TAUMAX WAS CROSSED')
      END
C
C
      SUBROUTINE A1TURN (N,N1,IPIVA,EPS,INFO,Y,Y0,Y1,YH,DYT,DYT0,DYT1,V,
     &   YW,IPIV,FY,TAUMIN,TAUMAX,D,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U)
C
C  SUBROUTINE  A1TURN  TO BE USED WITH ROUTINE ALCON1
C
C  DETERMINATION OF TURNING POINTS
C  INTERVAL METHOD WITH CUBIC HERMITE INTERPOLATION AS PREDICTOR
C  AND GAUSS-NEWTON METHOD AS CORRECTOR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),Y0(N1),Y1(N1),YH(N1),DYT(N1),DYT0(N1),
     $                 DYT1(N1)
      DOUBLE PRECISION V(N1),YW(N1)
      DOUBLE PRECISION D(N1),ETA(N1),DFY(N,N1),AQ(N,N1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(N1),DYQ(N1),W(N1)
      DOUBLE PRECISION F(N),FH(N),U(N)
      INTEGER INFO(9),IPIV(N1)
      INTEGER UPR,UDIAG
      EXTERNAL FY
      COMMON /A1UNIT/ UPR,UDIAG
      COMMON /A1VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C  PREPARATIONS
      KPRINT=INFO(1)
      INFO(9)=-1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INFO(4).EQ.0 .OR. INFO(4).EQ.2) GOTO 9990
      IPIVH=IPIVA
      IF (IPIVH.EQ.N1) IPIVH=IPIV(N1)
      DO 1000 I=1,N1
      Y0(I)=Y(I)
      DYT1(I)=DYT(I)
1000  CONTINUE
      F0P=DYT0(N1)/DYT0(IPIVH)
      F1P=DYT1(N1)/DYT1(IPIVH)
      INFO(1)=0
      ITER=0
      X0=Y0(IPIVH)/YW(IPIVH)
      X1=Y1(IPIVH)/YW(IPIVH)
      H=X1-X0
      DEL=DABS(H)
C-----------------------------------------------------------------------
C  ITERATION LOOP
2000  CONTINUE
C  HERMITE INTERPOLATION
      F0=Y0(N1)/YW(N1)
      F1=Y1(N1)/YW(N1)
      A=6.D0*(F0-F1)+3.D0*H*(F0P+F1P)
      B=-6.D0*(F0-F1)-2.D0*H*(2.D0*F0P+F1P)
      C=H*F0P
      Z=B**2-4.D0*A*C
      Z=DSIGN(DSQRT(Z),-B)
      Z=(-B+Z)/(2.D0*A)
      IF (Z.LT.0.D0 .OR. Z.GT.1.D0) Z=(C/A)/Z
      YH(IPIVH)=Y0(IPIVH)+H*Z*YW(IPIVH)
      DO 2020 I=1,N1
      IF (I.EQ.IPIVH) GOTO 2010
      F0H=Y0(I)/YW(I)
      F1H=Y1(I)/YW(I)
      F0PH=DYT0(I)/DYT0(IPIVH)
      F1PH=DYT1(I)/DYT1(IPIVH)
      A=2.D0*(F0H-F1H)+H*(F0PH+F1PH)
      B=-3.D0*(F0H-F1H)-H*(2.D0*F0PH+F1PH)
      C=H*F0PH
      R=F0H+Z*(C+Z*(B+Z*A))
      YH(I)=R*YW(I)
2010  CONTINUE
2020  CONTINUE
C  ITERATION BACK TO CONTINUATION PATH
      ITER=ITER+1
      IF (ITER.GT.ITMAX) GOTO 9910
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      INFO(8)=2
      CALL A1GNHM(N,N1,YH,FY,EPS,INFO,SUMYH,SUMYQH,TAUMIN,TAUMAX,
     &     IPIV,V,D,YW,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U)
      IF (INFO(9).NE.0) GOTO 9920
      CALL A1TAND(N,N1,FY,YH,F,FH,YW,ETA,DFY,AQ,COND,DETF,IPIV,AH,V,
     &   D,IRANK)
      IF (IRANK.LT.N) GOTO 9920
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      IF (DABS(YH(IPIVH)-Y0(IPIVH)).GT.DEL*YW(IPIVH) .OR.
     &    DABS(YH(IPIVH)-Y1(IPIVH)).GT.DEL*YW(IPIVH)) GOTO 9930
      DO 2030 I=1,N
      W(IPIV(I))=-V(I)/D(N1)
2030  CONTINUE
      W(IPIV(N1))=1.D0/D(N1)
      FQP=W(N1)/W(IPIVH)
C  NEW INCLUSION INTERVAL
      IF (F0P*FQP.GT.0.D0) GOTO 2200
      DQ=DABS(YH(IPIVH)-Y1(IPIVH))
      F1P=FQP
      DO 2110 I=1,N1
      Y1(I)=YH(I)
      DYT1(I)=W(I)
2110  CONTINUE
      GOTO 2300
2200  CONTINUE
      DQ=DABS(YH(IPIVH)-Y0(IPIVH))
      F0P=FQP
      DO 2210 I=1,N1
      Y0(I)=YH(I)
      DYT0(I)=W(I)
2210  CONTINUE
2300  CONTINUE
      X0=Y0(IPIVH)/YW(IPIVH)
      X1=Y1(IPIVH)/YW(IPIVH)
      H=X1-X0
      DEL=DABS(H)
      IF (KPRINT.GT.0 .AND. ITER.EQ.1) WRITE(UPR,60101)
      IF (KPRINT.GT.0) WRITE(UPR,60102) ITER,YH(IPIVH),DEL,FQP,INFO(7)
      IF (FQP.EQ.0.D0) GOTO 9000
      IF (DEL.LT.EPS) GOTO 9000
      IF (DQ/YW(IPIVH).LT.EPS) GOTO 9000
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9000  CONTINUE
      INFO(9)=0
      IF (KPRINT.GT.0) WRITE(UPR,60103) ITER
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(N1),(YH(I),I=1,N)
      GOTO 9990
C  FAIL EXIT
C  MORE THAN ITMAX ITERATIONS
9910  CONTINUE
      INFO(9)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      GOTO 9990
C  GAUSS-NEWTON METHOD FAILED
9920  CONTINUE
      INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      GOTO 9990
C  ITERATION STRATEGY FAILED
9930  CONTINUE
      INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60903)
C  RETURN TO CALLING ROUTINE
9990  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      INFO(1)=KPRINT
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(' ',132('-'))
60002 FORMAT(' POSSIBLY TURNING POINT DETECTED')
60003 FORMAT(/)
60101 FORMAT(/4X,'IGNC',7X,'YHIT',8X,'DELTAX',8X,'DERIV',6X,'ITER')
60102 FORMAT(5X,I2,4X,D11.4,3X,D10.4,3X,D11.4,4X,I2)
60103 FORMAT(/' TURN REQUIRED',I3,' GAUSS-NEWTON-CALLS'/)
60104 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60201 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,4(D14.6,6X)/))
60901 FORMAT(/' TURN TERMINATED AFTER ITMAX ITERATIONS')
60902 FORMAT(/' TURN TERMINATED SINCE GAUSS-NEWTON METHOD FAILED')
60903 FORMAT(/' TURN TERMINATED SINCE ITERATION STRATEGY FAILED')
      END
C
C
      SUBROUTINE A1BIFU (N,N1,IPIVA,EPS,INFO,Y,Y0,Y1,YH,V,YW,IPIV,FY,
     &   TAUMIN,TAUMAX,D,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U,DETFS,DETFAH)
C
C  SUBROUTINE  A1BIFU  TO BE USED WITH ROUTINE ALCON1
C
C  DETERMINATION OF BIFURCATION POINTS
C  INTERVAL METHOD WITH BRENT ALGORITHM IN CONNECTION WITH
C  LINEAR INTERPOLATION AS PREDICTOR
C  AND GAUSS-NEWTON METHOD AS CORRECTOR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),Y0(N1),Y1(N1),YH(N1)
      DOUBLE PRECISION V(N1),YW(N1)
      DOUBLE PRECISION D(N1),ETA(N1),DFY(N,N1),AQ(N,N1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(N1),DYQ(N1),W(N1)
      DOUBLE PRECISION F(N),FH(N),U(N)
      INTEGER INFO(9),IPIV(N1)
      INTEGER UPR,UDIAG
      EXTERNAL FY
      COMMON /A1UNIT/ UPR,UDIAG
      COMMON /A1VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C  PREPARATIONS
      KPRINT=INFO(1)
      INFO(9)=-1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INFO(4).LT.2) GOTO 9990
      DETF=DETFS
      DETFA=DETFAH
      DO 1000 I=1,N1
      Y0(I)=Y(I)
1000  CONTINUE
      INFO(1)=0
      ITER=0
      X0=Y0(IPIVA)/YW(IPIVA)
      X1=Y1(IPIVA)/YW(IPIVA)
      H=X1-X0
      DEL=DABS(H)
      AZ=X0
      FA=DETFA
      BZ=X1
      FB=DETF
C-----------------------------------------------------------------------
C  ITERATION LOOP
C  BRENT ALGORITHM
2000  CONTINUE
      CZ=AZ
      FC=FA
      DZ=BZ-CZ
      EZ=DZ
2010  IF (DABS(FC).GE.DABS(FB)) GOTO 2015
      AZ=BZ
      BZ=CZ
      CZ=AZ
      FA=FB
      FB=FC
      FC=FA
2015  CONTINUE
      RM=(CZ-BZ)*0.5D0
      IF (DABS(CZ-BZ).LE.EPS .AND. ITER.GT.0) GOTO 9000
      IF (DABS(EZ).LT.EPS) GOTO 2030
      IF (DABS(FA).LE.DABS(FB)) GOTO 2030
      SZ=FB/FA
      IF (AZ.NE.CZ) GOTO 2020
      PZ=(CZ-BZ)*SZ
      QZ=1.D0-SZ
      GOTO 2025
2020  QZ=FA/FC
      RZ=FB/FC
      RZ1=RZ-1.D0
      PZ=SZ*((CZ-BZ)*QZ*(QZ-RZ)-(BZ-AZ)*RZ1)
      QZ=(QZ-1.D0)*RZ1*(SZ-1.D0)
2025  IF (PZ.GT.0.D0) QZ=-QZ
      IF (PZ.LT.0.D0) PZ=-PZ
      SZ=EZ
      EZ=DZ
      IF (PZ+PZ.GE.3.D0*RM*QZ) GOTO 2030
      IF (PZ+PZ.GE.DABS(SZ*QZ)) GOTO 2030
      DZ=PZ/QZ
      GOTO 2035
2030  EZ=RM
      DZ=EZ
2035  AZ=BZ
      FA=FB
      TEMP=DZ
      IF (DABS(TEMP).LE.0.5D0*EPS) TEMP=DSIGN(0.5D0*EPS,RM)
      BZ=BZ+TEMP
      SZ=BZ
      Z=(SZ-AZ)/(CZ-AZ)
      IF (FA.EQ.DETF) Z=1.D0-Z
C  LINEAR INTERPOLATION
      DO 2080 I=1,N1
      F0H=Y0(I)/YW(I)
      F1H=Y1(I)/YW(I)
      R=F0H+Z*(F1H-F0H)
      YH(I)=R*YW(I)
2080  CONTINUE
C  ITERATION BACK TO CONTINUATION PATH
      ITER=ITER+1
      IF (ITER.GT.ITMAX) GOTO 9910
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      INFO(8)=2
      CALL A1GNHM(N,N1,YH,FY,EPS,INFO,SUMYH,SUMYQH,TAUMIN,TAUMAX,
     &     IPIV,V,D,YW,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U)
      IERR=INFO(9)
      INFO(9)=-1
      GOTO (9920,9920,9010,9920),IERR
      CALL A1TAND(N,N1,FY,YH,F,FH,YW,ETA,DFY,AQ,COND,DETFH,IPIV,AH,V,
     &   D,IRANK)
      IF (IRANK.LT.N) GOTO 9010
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      IF (DABS(YH(IPIVA)-Y0(IPIVA)).GT.DEL*YW(IPIVA) .OR.
     &    DABS(YH(IPIVA)-Y1(IPIVA)).GT.DEL*YW(IPIVA)) GOTO 9930
      DO 2090 I=1,N
      W(IPIV(I))=-V(I)/D(N1)
2090  CONTINUE
      W(IPIV(N1))=1.D0/D(N1)
      IF (IPIV(N1).NE.IPIVA) CALL A1DETH(N,N1,AQ,AH,D,IPIVA,IPIV,DETFH)
C  NEW INCLUSION INTERVAL
      IF (DETFA*DETFH.GT.0.D0) GOTO 2200
      DETF=DETFH
      DO 2110 I=1,N1
      Y1(I)=YH(I)
2110  CONTINUE
      GOTO 2300
2200  CONTINUE
      DETFA=DETFH
      DO 2210 I=1,N1
      Y0(I)=YH(I)
2210  CONTINUE
2300  CONTINUE
      X0=Y0(IPIVA)/YW(IPIVA)
      X1=Y1(IPIVA)/YW(IPIVA)
      H=X1-X0
      DEL=DABS(H)
      IF (KPRINT.GT.0 .AND. ITER.EQ.1) WRITE(UPR,60101)
      IF (KPRINT.GT.0) WRITE(UPR,60102) ITER,YH(IPIVA),DEL,DETFH,INFO(7)
      FB=DETFH
      IF (FB*FC.LE.0.D0) GOTO 2010
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9000  CONTINUE
      INFO(9)=0
9010  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60103) ITER
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(N1),(YH(I),I=1,N)
      IF (KPRINT.GT.0 .AND. IRANK.LT.N) WRITE(UPR,60109)
      GOTO 9990
C  FAIL EXIT
C  MORE THAN ITMAX ITERATIONS
9910  CONTINUE
      INFO(9)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      GOTO 9990
C  GAUSS-NEWTON METHOD FAILED
9920  CONTINUE
      INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      GOTO 9990
C  ITERATION STRATEGY FAILED
9930  CONTINUE
      INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60903)
C  RETURN TO CALLING ROUTINE
9990  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      INFO(1)=KPRINT
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(' ',132('-'))
60002 FORMAT(' POSSIBLY BIFURCATION POINT DETECTED')
60003 FORMAT(/)
60101 FORMAT(/4X,'IGNC',7X,'YHIT',8X,'DELTAX',8X,'DETER',6X,'ITER')
60102 FORMAT(5X,I2,4X,D11.4,3X,D10.4,3X,D11.4,4X,I2)
60103 FORMAT(/' BIFURC REQUIRED',I3,' GAUSS-NEWTON-CALLS'/)
60104 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60109 FORMAT(/' NO TANGENT AT SOLUTION AVAILABLE')
60201 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60901 FORMAT(/' BIFURC TERMINATED AFTER ITMAX ITERATIONS')
60902 FORMAT(/' BIFURC TERMINATED SINCE GAUSS-NEWTON METHOD FAILED')
60903 FORMAT(/' BIFURC TERMINATED SINCE ITERATION STRATEGY FAILED')
      END
C
C
      SUBROUTINE A1DERI (FY,N,N1,Y,F,FH,YW,ETA,DFY)
C
C  SUBROUTINE  A1DERI  TO BE USED WITH ROUTINE ALCON1
C
C  COMPUTATION OF JACOBIAN F' BY FINITE DIFFERENCES
C  INCLUDING FEED BACK DEVICE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),F(N),FH(N),YW(N1),ETA(N1),DFY(N,N1)
      EXTERNAL FY
      COMMON /COUNT/ IFCTEV,IDUMQ,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A1VALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      EPDIFF=1.D3*EPMACH
      ETAD=DSQRT(1.D1*EPMACH)
      DO 1 K=1,N1
      K1=0
11    T=Y(K)
      ETAH=ETA(K)
      IF (T.LT.0.D0) ETAH=-ETAH
      S=ETAH*YW(K)
      Y(K)=T+S
      CALL FY(Y,Y(N1),FH)
      IFCTEV=IFCTEV+1
      Y(K)=T
      SU=0.D0
      DO 13 I=1,N
      HG=DMAX1(DABS(F(I)),DABS(FH(I)))
      T=FH(I)-F(I)
      IF (HG.NE.0.D0) SU=SU+(T/HG)**2
13    DFY(I,K)=T/ETAH
      SU=DSQRT(SU/DBLE(N1))
      IF ((SU.EQ.0.D0).OR.(K1.GT.0)) GOTO 1
      ETA(K)=DMAX1(EPDIFF,DSQRT(ETAD/SU)*ETA(K))
      ETA(K)=DMIN1(SMALL1,ETA(K))
      K1=1
      IF (SU.LT.EPDIFF) GOTO 11
1     CONTINUE
      RETURN
      END
C
C
      SUBROUTINE A1DETH (N,N1,A,AH,D,IPIVA,IPIV,DET)
C
C  SUBROUTINE  A1DETH  TO BE USED WITH ROUTINE ALCON1
C
C  COMPUTATION OF DETERMINANT OF A HESSENBERG MATRIX
C  BY GIVENS ROTATIONS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(N,N1),AH(N,N1),D(N1),IPIV(N1)
C  SEARCH OLD PIVOT
      DO 1010 I=1,N1
      IQ=I
      IF (IPIV(I).EQ.IPIVA) GOTO 1020
1010  CONTINUE
1020  CONTINUE
C  PRODUCT OF FIRST IQ-1 DIAGONAL ELEMENTS
      DET=1.D0
      IF (IQ.EQ.1) GOTO 2011
      IEND=IQ-1
      DO 2010 I=1,IEND
      DET=DET*D(I)
2010  CONTINUE
2011  CONTINUE
C  N-IQ GIVENS TRANSFORMATIONS
      IF (IQ.EQ.N1) GOTO 3050
      IBEG=IQ+1
      X1=A(IBEG-1,IBEG)
      IF (IBEG.GT.N) GOTO 3040
      DO 3010 I=IBEG,N1
      AH(IQ,I)=A(IQ,I)
3010  CONTINUE
      DO 3030 I=IBEG,N
      Y1=D(I)
      R=DSQRT(X1**2+Y1**2)
      C=X1/R
      S=Y1/R
      DET=DET*(C*X1+S*Y1)
      I1=I+1
      DO 3020 J=I1,N1
      X2=AH(I-1,J)
      Y2=A(I,J)
      AH(I,J)=-S*X2+C*Y2
3020  CONTINUE
      X1=AH(I,I1)
3030  CONTINUE
3040  CONTINUE
      DET=DET*X1
3050  CONTINUE
C  SIGN OF PIVOT VECTOR
      DO 4020 I=1,N
      IF (I.EQ.IQ) GOTO 4020
      I1=I+1
      DO 4010 J=I1,N1
      IF (J.EQ.IQ) GOTO 4010
      IF (IPIV(J).LT.IPIV(I)) DET=-DET
4010  CONTINUE
4020  CONTINUE
      RETURN
      END
C
C
      SUBROUTINE A1TAND (N,N1,FY,Y,F,FH,YW,ETA,DFY,A,COND,DET,PIVOT,AH,
     &   V,D,IRANK)
C
C  SUBROUTINE  A1TAND  TO BE USED WITH ROUTINE ALCON1
C
C  COMPUTATION OF TANGENT AND DETERMINANT AT GIVEN POINT OF
C  CONTINUATION PATH BY QR DECOMPOSITION
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),F(N),FH(N),YW(N1),ETA(N1),DFY(N,N1),A(N,N1)
      DOUBLE PRECISION AH(N1,N1),V(N1),D(N1)
      INTEGER PIVOT(N1)
      EXTERNAL FY
      COMMON /COUNT/ IFCTEV,IDUMQ,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
C  COMPUTE NEW TANGENT AND DETERMINANT
      CALL FY(Y,Y(N1),F)
      IFCTEV=IFCTEV+1
      CALL A1DERI(FY,N,N1,Y,F,FH,YW,ETA,DFY)
      DO 1010 I=1,N
      DO 1010 J=1,N1
      A(I,J)=DFY(I,J)
1010  CONTINUE
      IRANK=N
      COND=1.D0/EPMACH
      CALL A1DECC(A,N,N1,0,N,N1,IRANK,COND,D,PIVOT,0,AH,V)
      IDECS=IDECS+1
      IF (IRANK.LT.N) RETURN
      DET=1.D0
      DO 1020 I=1,N
      DET=DET*D(I)
1020  CONTINUE
      NM1=N-1
      IF (NM1.EQ.0) GOTO 1031
      DO 1030 I=1,NM1
      I1=I+1
      DO 1030 J=I1,N
      IF (PIVOT(J).LT.PIVOT(I)) DET=-DET
1030  CONTINUE
1031  CONTINUE
      RETURN
      END
C
C
      BLOCK DATA A1BLOK
C
C  BLOCK DATA  TO BE USED WITH ROUTINE ALCON1
C
C    COMMON /A1UNIT/ UPR,UDIAG,UBIF
C      INPUT/OUTPUT UNITS
C      - UPR     PRINT UNIT                            OUTPUT
C      - UDIAG   UNIT FOR PLOT INFORMATION             OUTPUT
C                (IN CONNECTION WITH ROUTINE PLTHM)
C
C    COMMON /A1VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C      LIMIT VALUES IN CONTINUATION PROCESS
C      - SMALL1  USED AS INITIAL STEPLENGTH ESTIMATE AND AS MAXIMUM
C                PERMITTED STEPSIZE IN NUMERICAL DIFFERENTIATION
C      - THMAX   MAXIMUM PERMITTED RATIO IN RESTRICTED MONOTONICITY TEST
C      - THR     THRESHOLD IN NEARLY LINEAR CASE
C      - ITMAX   MAXIMUM PERMITTED NUMBER OF ITERATIONS
C                IN INTERVAL METHODS
C      - IGNMAX  MAXIMUM PERMITTED NUMBER OF GAUSS-NEWTON ITERATIONS
C                PER CORRECTOR CALL
C
C    COMMON /A1METH/ FACTOR,IEXTR,IJACM
C      METHOD DESCRIBING PARAMETERS
C      - FACTOR  FACTOR FOR EXTRAPOLATED STEPLENGTH BOUND
C      - IEXTR   METHOD FOR EXTRAPOLATED STEPLENGTH BOUND
C                  0  NO EXTRAPOLATED STEPLENGTH BOUND
C                  1  WITH EXTRAPOLATED STEPLENGTH BOUND
C      - IJACM   METHOD FOR UPDATING JACOBIAN IN GAUSS-NEWTON METHOD
C                  0  KEEP FIRST JACOBIAN
C                  1  BROYDEN UPDATES OF JACOBIAN
C                  2  NEW JACOBIAN IN EACH ITERATION
C                  3  AS 1  BUT IN A COMPUTING TIME SAVING WAY
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER UPR,UDIAG
      COMMON /A1UNIT/ UPR,UDIAG
      COMMON /A1VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
      COMMON /A1METH/ FACTOR,IEXTR,IJACM
      DATA UPR/6/,UDIAG/2/
      DATA SMALL1/0.01D0/,THMAX/0.25D0/,THR/0.1D0/,ITMAX/20/,IGNMAX/10/
      DATA FACTOR/1.1D0/,IEXTR/1/,IJACM/3/
      END
