SUBROUTINE STEPNS(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR,
     &   ABSERR,S,Y,YP,YOLD,YPOLD,A,QR,LENQR,PIVOT,WORK,SSPAR,
     &   PAR,IPAR)
C
C  STEPNS  TAKES ONE STEP ALONG THE ZERO CURVE OF THE HOMOTOPY MAP
C USING A PREDICTOR-CORRECTOR ALGORITHM.  THE PREDICTOR USES A HERMITE
C CUBIC INTERPOLANT, AND THE CORRECTOR RETURNS TO THE ZERO CURVE ALONG
C THE FLOW NORMAL TO THE DAVIDENKO FLOW.  STEPNS  ALSO ESTIMATES A
C STEP SIZE H FOR THE NEXT STEP ALONG THE ZERO CURVE.  NORMALLY
C  STEPNS  IS USED INDIRECTLY THROUGH  FIXPNS , AND SHOULD BE CALLED
C DIRECTLY ONLY IF IT IS NECESSARY TO MODIFY THE STEPPING ALGORITHM'S
C PARAMETERS.
C
C ON INPUT:
C
C N = DIMENSION OF X AND THE HOMOTOPY MAP.
C
C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
C
C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
C
C START = .TRUE. ON FIRST CALL TO  STEPNS , .FALSE. OTHERWISE.
C
C HOLD = ||Y - YOLD||; SHOULD NOT BE MODIFIED BY THE USER.
C
C H = UPPER LIMIT ON LENGTH OF STEP THAT WILL BE ATTEMPTED.  H  MUST BE
C    SET TO A POSITIVE NUMBER ON THE FIRST CALL TO  STEPNS .
C    THEREAFTER  STEPNS  CALCULATES AN OPTIMAL VALUE FOR  H , AND  H
C    SHOULD NOT BE MODIFIED BY THE USER.
C
C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS
C    CONSIDERED TO HAVE CONVERGED WHEN A POINT W=(X,LAMBDA) IS FOUND 
C    SUCH THAT
C
C    ||Z|| <= RELERR*||W|| + ABSERR  ,          WHERE
C
C    Z IS THE NEWTON STEP TO W=(X,LAMBDA).
C
C S = (APPROXIMATE) ARC LENGTH ALONG THE HOMOTOPY ZERO CURVE UP TO
C    Y(S) = (X(S), LAMBDA(S)).
C
C Y(1:N+1) = PREVIOUS POINT (X(S), LAMBDA(S)) FOUND ON THE ZERO CURVE OF 
C    THE HOMOTOPY MAP.
C
C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
C    AT  Y .
C
C YOLD(1:N+1) = A POINT BEFORE  Y  ON THE ZERO CURVE OF THE HOMOTOPY MAP.
C
C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY
C    MAP AT  YOLD .
C
C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
C
C QR(1:LENQR) = THE N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X
C    STORED IN PACKED SKYLINE STORAGE FORMAT.  LENQR  AND  PIVOT  
C    DESCRIBE THE DATA STRUCTURE IN  QR .
C
C LENQR = LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR  USED TO CONTAIN THE
C    N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X IN PACKED
C    SKYLINE STORAGE FORMAT.
C
C PIVOT(1:N+2) = INDICES OF THE DIAGONAL ELEMENTS OF THE N X N SYMMETRIC
C    JACOBIAN MATRIX (WITH RESPECT TO X) WITHIN  QR .
C
C WORK(1:13*(N+1)+2*N+LENQR) = WORK ARRAY SPLIT UP AND USED FOR THE
C    CALCULATION OF THE JACOBIAN MATRIX KERNEL, THE NEWTON STEP,
C    INTERPOLATION, AND THE ESTIMATION OF THE NEXT STEP SIZE  H .
C
C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P)  IS
C    A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION.
C
C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHO, RHOJS.
C
C ON OUTPUT:
C
C N , A , SSPAR  ARE UNCHANGED.
C
C NFE  HAS BEEN UPDATED.
C
C IFLAG  
C    = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
C
C    = 4 IF THE CONJUGATE GRADIENT ITERATION FAILED TO CONVERGE
C        (MOST LIKELY DUE TO A JACOBIAN MATRIX WITH RANK < N).  THE
C        ITERATION WAS NOT COMPLETED.
C
C    = 6 IF THE NEWTON ITERATION FAILED TO CONVERGE.  W  CONTAINS 
C        THE LAST NEWTON ITERATE.
C
C START = .FALSE. ON A NORMAL RETURN.
C
C CRASH 
C    = .FALSE. ON A NORMAL RETURN.
C
C    = .TRUE. IF THE STEP SIZE  H  WAS TOO SMALL.  H  HAS BEEN
C      INCREASED TO AN ACCEPTABLE VALUE, WITH WHICH  STEPNS  MAY BE
C      CALLED AGAIN.
C
C    = .TRUE. IF  RELERR  AND/OR  ABSERR  WERE TOO SMALL.  THEY HAVE
C      BEEN INCREASED TO ACCEPTABLE VALUES, WITH WHICH  STEPNS  MAY
C      BE CALLED AGAIN.
C
C HOLD = ||Y - YOLD||.
C
C H = OPTIMAL VALUE FOR NEXT STEP TO BE ATTEMPTED.  NORMALLY  H  SHOULD
C    NOT BE MODIFIED BY THE USER.
C
C RELERR, ABSERR  ARE UNCHANGED ON A NORMAL RETURN.
C
C S = (APPROXIMATE) ARC LENGTH ALONG THE ZERO CURVE OF THE HOMOTOPY MAP 
C    UP TO THE LATEST POINT FOUND, WHICH IS RETURNED IN  Y .
C
C Y, YP, YOLD, YPOLD  CONTAIN THE TWO MOST RECENT POINTS AND TANGENT
C    VECTORS FOUND ON THE ZERO CURVE OF THE HOMOTOPY MAP.
C
C
C CALLS  D1MACH , DAXPY , DCOPY , DNRM2 , TANGNS .
C
      DOUBLE PRECISION ABSERR,D1MACH,DCALC,DD001,DD0011,DD01,
     &   DD011,DELS,DNRM2,F0,F1,FOURU,FP0,FP1,H,HFAIL,HOLD,HT,
     &   LCALC,QOFS,RCALC,RELERR,RHOLEN,S,TEMP,TWOU
      INTEGER IFLAG,IPP,IRHO,ITANGW,ITNUM,ITZ,IW,IWP,IZ0,IZ1,
     &   J,JUDY,LENQR,LITFH,N,NFE,NP1
      LOGICAL CRASH,FAIL,START
C
C ***** ARRAY DECLARATIONS. *****
C
      DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N),
     &   QR(LENQR),WORK(13*(N+1)+2*N+LENQR),SSPAR(8),PAR(1)
      INTEGER PIVOT(N+2),IPAR(1)
C
C ***** END OF DIMENSIONAL INFORMATION. *****
C
C THE LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED BEFORE REDUCING
C THE STEP SIZE  H  MAY BE CHANGED BY CHANGING THE FOLLOWING PARAMETER 
C STATEMENT:
      PARAMETER (LITFH=4)
C
C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
C
      DD01(F0,F1,DELS)=(F1-F0)/DELS
      DD001(F0,FP0,F1,DELS)=(DD01(F0,F1,DELS)-FP0)/DELS
      DD011(F0,F1,FP1,DELS)=(FP1-DD01(F0,F1,DELS))/DELS
      DD0011(F0,FP0,F1,FP1,DELS)=(DD011(F0,F1,FP1,DELS) - 
     &                            DD001(F0,FP0,F1,DELS))/DELS
      QOFS(F0,FP0,F1,FP1,DELS,S)=((DD0011(F0,FP0,F1,FP1,DELS)*(S-DELS) +
     &   DD001(F0,FP0,F1,DELS))*S + FP0)*S + F0
C
C
      TWOU=2.0*D1MACH(4)
      FOURU=TWOU+TWOU
      NP1=N+1
      IPP=1
      IRHO=N+1
      IW=IRHO+N
      IWP=IW+NP1
      ITZ=IWP+NP1
      IZ0=ITZ+NP1
      IZ1=IZ0+NP1
      ITANGW=IZ1+NP1
      CRASH=.TRUE.
C THE ARCLENGTH  S  MUST BE NONNEGATIVE.
      IF (S .LT. 0.0) RETURN
C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE.
      IF (H .LT. FOURU*(1.0+S)) THEN
        H=FOURU*(1.0+S)
        RETURN
      ENDIF
C IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE VALUES.
      TEMP=DNRM2(NP1,Y,1)
      IF (.5*(RELERR*TEMP+ABSERR) .GE. TWOU*TEMP) GO TO 40
      IF (RELERR .NE. 0.0) THEN
        RELERR=FOURU*(1.0+FOURU)
        ABSERR=MAX(ABSERR,0.0D0)
      ELSE
        ABSERR=FOURU*TEMP
      ENDIF
      RETURN
 40   CRASH=.FALSE.
      IF (.NOT. START) GO TO 300
C
C *****  STARTUP SECTION(FIRST STEP ALONG ZERO CURVE.  *****
C
      FAIL=.FALSE.
      START=.FALSE.
C DETERMINE SUITABLE INITIAL STEP SIZE.
      H=MIN(H, .10D0, SQRT(SQRT(RELERR*TEMP+ABSERR)))
C USE LINEAR PREDICTOR ALONG TANGENT DIRECTION TO START NEWTON ITERATION.
      YPOLD(NP1)=1.0
      DO 50 J=1,N
        YPOLD(J)=0.0
 50   CONTINUE
      CALL TANGNS(S,Y,YP,WORK(ITZ),YPOLD,A,QR,LENQR,PIVOT,
     &   WORK(IPP),WORK(IRHO),WORK(ITANGW),NFE,N,IFLAG,PAR,IPAR)
      IF (IFLAG .GT. 0) RETURN
 70   DO 80 J=1,NP1
        TEMP=Y(J) + H * YP(J)
        WORK(IW+J-1)=TEMP
        WORK(IZ0+J-1)=TEMP
 80   CONTINUE
      DO 200 JUDY=1,LITFH
        RHOLEN=-1.0
C CALCULATE THE NEWTON STEP  TZ  AT THE CURRENT POINT  W .
        CALL TANGNS(RHOLEN,WORK(IW),WORK(IWP),WORK(ITZ),YPOLD,A,
     &     QR,LENQR,PIVOT,WORK(IPP),WORK(IRHO),WORK(ITANGW),
     &     NFE,N,IFLAG,PAR,IPAR)
        IF (IFLAG .GT. 0) RETURN
C
C TAKE NEWTON STEP AND CHECK CONVERGENCE.
        CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1)
        ITNUM=JUDY
C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
        IF (JUDY .EQ. 1) THEN
          LCALC=DNRM2(NP1,WORK(ITZ),1)
          RCALC=RHOLEN
          CALL DCOPY(NP1,WORK(IW),1,WORK(IZ1),1)
        ELSE IF (JUDY .EQ. 2) THEN
          LCALC=DNRM2(NP1,WORK(ITZ),1)/LCALC
          RCALC=RHOLEN/RCALC
        ENDIF
C GO TO MOP-UP SECTION AFTER CONVERGENCE.
        IF ( DNRM2(NP1,WORK(ITZ),1) .LE. 
     &     RELERR*DNRM2(NP1,WORK(IW),1)+ABSERR )  GO TO 600
C
200   CONTINUE
C
C NO CONVERGENCE IN  LITFH  ITERATIONS.  REDUCE  H  AND TRY AGAIN.
      IF (H .LE. FOURU*(1.0 + S)) THEN
        IFLAG=6
        RETURN
      ENDIF
      H=.5 * H
      GO TO 70
C
C ***** END OF STARTUP SECTION. *****
C
C ***** PREDICTOR SECTION. *****
C
300   FAIL=.FALSE.
C COMPUTE POINT PREDICTED BY HERMITE INTERPOLANT.  USE STEP SIZE  H
C COMPUTED ON LAST CALL TO  STEPNS .
320   DO 330 J=1,NP1
        TEMP=QOFS(YOLD(J),YPOLD(J),Y(J),YP(J),HOLD,HOLD+H)
        WORK(IW+J-1)=TEMP
        WORK(IZ0+J-1)=TEMP
330   CONTINUE
C
C ***** END OF PREDICTOR SECTION. *****
C
C ***** CORRECTOR SECTION. *****
C
      DO 500 JUDY=1,LITFH
        RHOLEN=-1.0
C CALCULATE THE NEWTON STEP  TZ  AT THE CURRENT POINT  W .
        CALL TANGNS(RHOLEN,WORK(IW),WORK(IWP),WORK(ITZ),YP,A,
     &     QR,LENQR,PIVOT,WORK(IPP),WORK(IRHO),WORK(ITANGW),
     &     NFE,N,IFLAG,PAR,IPAR)
        IF (IFLAG .GT. 0) RETURN
C
C TAKE NEWTON STEP AND CHECK CONVERGENCE.
        CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1)
        ITNUM=JUDY
C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
        IF (JUDY .EQ. 1) THEN
          LCALC=DNRM2(NP1,WORK(ITZ),1)
          RCALC=RHOLEN
          CALL DCOPY(NP1,WORK(IW),1,WORK(IZ1),1)
        ELSE IF (JUDY .EQ. 2) THEN
          LCALC=DNRM2(NP1,WORK(ITZ),1)/LCALC
          RCALC=RHOLEN/RCALC
        ENDIF
C GO TO MOP-UP SECTION AFTER CONVERGENCE.
        IF ( DNRM2(NP1,WORK(ITZ),1) .LE. 
     &     RELERR*DNRM2(NP1,WORK(IW),1)+ABSERR )  GO TO 600
C
500   CONTINUE
C
C NO CONVERGENCE IN  LITFH  ITERATIONS.  RECORD FAILURE AT CALCULATED  H , 
C SAVE THIS STEP SIZE, REDUCE  H  AND TRY AGAIN.
      FAIL=.TRUE.
      HFAIL=H
      IF (H .LE. FOURU*(1.0 + S)) THEN
        IFLAG=6
        RETURN
      ENDIF
      H=.5 * H
      GO TO 320
C
C ***** END OF CORRECTOR SECTION. *****
C
C ***** MOP-UP SECTION. *****
C
C YOLD  AND  Y  ALWAYS CONTAIN THE LAST TWO POINTS FOUND ON THE ZERO
C CURVE OF THE HOMOTOPY MAP.  YPOLD  AND  YP  CONTAIN THE TANGENT
C VECTORS TO THE ZERO CURVE AT  YOLD  AND  Y , RESPECTIVELY.
C
600   CALL DCOPY(NP1,Y,1,YOLD,1)
      CALL DCOPY(NP1,YP,1,YPOLD,1)
      CALL DCOPY(NP1,WORK(IW),1,Y,1)
      CALL DCOPY(NP1,WORK(IWP),1,YP,1)
      CALL DAXPY(NP1,-1.0D0,YOLD,1,WORK(IW),1)
C UPDATE ARC LENGTH.
      HOLD=DNRM2(NP1,WORK(IW),1)
      S=S+HOLD
C
C ***** END OF MOP-UP SECTION. *****
C
C ***** OPTIMAL STEP SIZE ESTIMATION SECTION. *****
C
C CALCULATE THE DISTANCE FACTOR  DCALC .
700   CALL DAXPY(NP1,-1.0D0,Y,1,WORK(IZ0),1)
      CALL DAXPY(NP1,-1.0D0,Y,1,WORK(IZ1),1)
      DCALC=DNRM2(NP1,WORK(IZ0),1)
      IF (DCALC .NE. 0.0) DCALC=DNRM2(NP1,WORK(IZ1),1)/DCALC
C
C THE OPTIMAL STEP SIZE HBAR IS DEFINED BY
C
C   HT=HOLD * [MIN(LIDEAL/LCALC, RIDEAL/RCALC, DIDEAL/DCALC)]**(1/P)
C
C     HBAR = MIN [ MAX(HT, BMIN*HOLD, HMIN), BMAX*HOLD, HMAX ]
C
C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, SET THE CONTRACTION
C FACTOR  LCALC  TO ZERO.
      IF (ITNUM .EQ. 1) LCALC = 0.0
C FORMULA FOR OPTIMAL STEP SIZE.
      IF (LCALC+RCALC+DCALC .EQ. 0.0) THEN
        HT = SSPAR(7) * HOLD
      ELSE 
        HT = (1.0/MAX(LCALC/SSPAR(1), RCALC/SSPAR(2), DCALC/SSPAR(3)))
     &       **(1.0/SSPAR(8)) * HOLD
      ENDIF
C  HT  CONTAINS THE ESTIMATED OPTIMAL STEP SIZE.  NOW PUT IT WITHIN
C REASONABLE BOUNDS.
      H=MIN(MAX(HT,SSPAR(6)*HOLD,SSPAR(4)), SSPAR(7)*HOLD, SSPAR(5))
      IF (ITNUM .EQ. 1) THEN
C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, DON'T DECREASE  H .
        H=MAX(H,HOLD)
      ELSE IF (ITNUM .EQ. LITFH) THEN
C IF CONVERGENCE REQUIRED THE MAXIMUM  LITFH  ITERATIONS, DON'T
C INCREASE  H .
        H=MIN(H,HOLD)
      ENDIF
C IF CONVERGENCE DID NOT OCCUR IN  LITFH  ITERATIONS FOR A PARTICULAR
C H = HFAIL , DON'T CHOOSE THE NEW STEP SIZE LARGER THAN  HFAIL .
      IF (FAIL) H=MIN(H,HFAIL)
C
C
      RETURN
      END