SUBROUTINE ROOTNS(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,YPOLD,
     $   A,QR,LENQR,PIVOT,WORK,PAR,IPAR)
C
C ROOTNS  FINDS THE POINT  YBAR = (XBAR, 1)  ON THE ZERO CURVE OF THE
C HOMOTOPY MAP.  IT STARTS WITH TWO POINTS YOLD=(XOLD,LAMBDAOLD) AND
C Y=(X,LAMBDA) SUCH THAT  LAMBDAOLD < 1 <= LAMBDA , AND ALTERNATES
C BETWEEN HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION UNTIL
C CONVERGENCE.
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 RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS
C    CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(X,LAMBDA) IS FOUND
C    SUCH THAT
C
C    |Y(NP1) - 1| <= RELERR + ABSERR              AND
C
C    ||Z|| <= RELERR*||X|| + ABSERR  ,          WHERE
C
C    (Z,?) IS THE NEWTON STEP TO Y=(X,LAMBDA).
C
C Y(1:N+1) = POINT (X(S), LAMBDA(S)) ON ZERO CURVE OF 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 DIFFERENT FROM  Y  ON THE ZERO CURVE.
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    AND INTERPOLATION.
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 , RELERR , ABSERR , A  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 PRECONDITIONED CONJUGATE GRADIENT ITERATION FAILED TO
C        CONVERGE (MOST LIKELY DUE TO A JACOBIAN MATRIX WITH RANK < N).
C        THE ITERATION WAS NOT COMPLETED.
C
C    = 6 IF THE INTERPOLATION/NEWTON ITERATION FAILED TO CONVERGE.
C        Y  AND  YOLD  CONTAIN THE LAST TWO POINTS FOUND ON THE
C        ZERO CURVE.
C
C Y  IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT  LAMBDA = 1 .
C
C
C CALLS  D1MACH , DAXPY , DCOPY , DNRM2 , ROOT , TANGNS .
C
      DOUBLE PRECISION ABSERR,AERR,D1MACH,DD001,DD0011,DD01,DD011,
     $   DELS,DNRM2,F0,F1,FP0,FP1,QOFS,QSOUT,RELERR,RERR,S,SA,SB,
     $   SOUT,U
      INTEGER IFLAG,IPP,IRHO,ITANGW,ITZ,IW,IWP,IZ0,IZ1,JUDY,JW,
     $   LCODE,LENQR,LIMIT,N,NFE,NP1
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),PAR(1)
      INTEGER PIVOT(N+2),IPAR(1)
C
C ***** END OF DIMENSIONAL INFORMATION. *****
C
C THE LIMIT ON THE NUMBER OF ITERATIONS ALLOWED MAY BE CHANGED BY
C CHANGING THE FOLLOWING PARAMETER STATEMENT:
      PARAMETER (LIMIT=20)
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
      U=D1MACH(4)
      RERR=MAX(RELERR,U)
      AERR=MAX(ABSERR,0.0D0)
      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
C
C *****  MAIN LOOP.  *****
C
100   DO 300 JUDY=1,LIMIT
      DO 110 JW=1,NP1
        WORK(ITZ+JW-1)=Y(JW)-YOLD(JW)
110   CONTINUE
      DELS=DNRM2(NP1,WORK(ITZ),1)
C
C USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT
C THE HERMITE CUBIC INTERPOLANT Q(S).  THEN USE  ROOT  TO FIND THE S
C CORRESPONDING TO  LAMBDA = 1 .  THE TWO POINTS ON THE ZERO CURVE ARE
C ALWAYS CHOSEN TO BRACKET LAMBDA=1, WITH THE BRACKETING INTERVAL
C ALWAYS BEING [0, DELS].
C
      SA=0.0
      SB=DELS
      LCODE=1
130   CALL ROOT(SOUT,QSOUT,SA,SB,RERR,AERR,LCODE)
      IF (LCODE .GT. 0) GO TO 140
      QSOUT=QOFS(YOLD(NP1),YPOLD(NP1),Y(NP1),YP(NP1),DELS,SOUT) - 1.0
      GO TO 130
C IF LAMBDA = 1 WERE BRACKETED,  ROOT  CANNOT FAIL.
140   IF (LCODE .GT. 2) THEN
        IFLAG=6
        RETURN
      ENDIF
C
C CALCULATE Q(SA) AS THE INITIAL POINT FOR A NEWTON ITERATION.
      DO 150 JW=1,NP1
        WORK(IW+JW-1)=QOFS(YOLD(JW),YPOLD(JW),Y(JW),YP(JW),DELS,SA)
150   CONTINUE
C CALCULATE NEWTON STEP AT Q(SA).
      CALL TANGNS(SA,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 NEXT POINT = CURRENT POINT + NEWTON STEP.
      CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1)
C GET THE TANGENT  WP  AT  W  AND THE NEXT NEWTON STEP IN  TZ .
      CALL TANGNS(SA,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 TAKE NEWTON STEP AND CHECK CONVERGENCE.
      CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1)
      IF ((ABS(WORK(IW+N)-1.0) .LE. RERR+AERR) .AND.
     $   (DNRM2(N,WORK(ITZ),1) .LE. RERR*DNRM2(N,WORK(IW),1)+AERR))
     $   THEN
        CALL DCOPY(NP1,WORK(IW),1,Y,1)
        RETURN
      ENDIF
C IF THE ITERATION HAS NOT CONVERGED, DISCARD ONE OF THE OLD POINTS
C SUCH THAT  LAMBDA = 1  IS STILL BRACKETED.
      IF ((YOLD(NP1)-1.0)*(WORK(IW+N)-1.0) .GT. 0.0) THEN
        CALL DCOPY(NP1,WORK(IW),1,YOLD,1)
        CALL DCOPY(NP1,WORK(IWP),1,YPOLD,1)
      ELSE
        CALL DCOPY(NP1,WORK(IW),1,Y,1)
        CALL DCOPY(NP1,WORK(IWP),1,YP,1)
      ENDIF
300   CONTINUE
C
C ***** END OF MAIN LOOP. *****
C
C THE ALTERNATING OSCULATORY CUBIC INTERPOLATION AND NEWTON ITERATION
C HAS NOT CONVERGED IN  LIMIT  STEPS.  ERROR RETURN.
      IFLAG=6
      RETURN
      END

.