share/hompack/fortran/fixpns.f

   1       SUBROUTINE FIXPNS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,
   2      $   NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,WORK,SSPAR,
   3      $   PAR,IPAR)
   4 C
   5 C SUBROUTINE  FIXPNS  FINDS A FIXED POINT OR ZERO OF THE
   6 C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE
   7 C OF A GENERAL HOMOTOPY MAP RHO(A,X,LAMBDA).  FOR THE FIXED
   8 C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL
   9 C INTO ITSELF.  THE EQUATION  X = F(X)  IS SOLVED BY
  10 C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP
  11 C
  12 C  LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A)  ,
  13 C
  14 C STARTING FROM LAMBDA = 0, X = A.  THE CURVE IS PARAMETERIZED
  15 C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY
  16 C DIFFERENTIAL EQUATION  D(HOMOTOPY MAP)/DS = 0  FOR
  17 C Y(S) = (X(S), LAMBDA(S)) USING A HERMITE CUBIC PREDICTOR AND A
  18 C CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL
  19 C TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF
  20 C D(HOMOTOPY MAP)/DS ).
  21 C
  22 C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP
  23 C SUCH THAT FOR SOME R > 0,  X*F(X) >= 0  WHENEVER NORM(X) = R.
  24 C THE EQUATION  F(X) = 0  IS SOLVED BY FOLLOWING THE ZERO CURVE
  25 C OF THE HOMOTOPY MAP
  26 C
  27 C   LAMBDA*F(X) + (1 - LAMBDA)*(X - A)
  28 C
  29 C EMANATING FROM LAMBDA = 0, X = A.
  30 C
  31 C  A  MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS.
  32 C
  33 C FOR THE CURVE TRACKING PROBLEM RHO(A,X,LAMBDA) IS ASSUMED TO
  34 C BE A C2 MAP FROM E**M X E**N X [0,1) INTO E**N, WHICH FOR
  35 C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET
  36 C OF E**M SATISFIES
  37 C
  38 C  RANK [D RHO(A,X,LAMBDA)/DX , D RHO(A,X,LAMBDA)/D LAMBDA] = N
  39 C
  40 C FOR ALL POINTS (X,LAMBDA) SUCH THAT RHO(A,X,LAMBDA)=0.  IT IS
  41 C FURTHER ASSUMED THAT
  42 C
  43 C           RANK [ D RHO(A,X0,0)/DX ] = N  .
  44 C
  45 C WITH A FIXED, THE ZERO CURVE OF RHO(A,X,LAMBDA) EMANATING
  46 C FROM  LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY
  47 C SOLVING THE ORDINARY DIFFERENTIAL EQUATION
  48 C D RHO(A,X(S),LAMBDA(S))/DS = 0  FOR  Y(S) = (X(S), LAMBDA(S)),
  49 C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE.  ALSO THE HOMOTOPY
  50 C MAP RHO(A,X,LAMBDA) IS ASSUMED TO BE CONSTRUCTED SUCH THAT
  51 C
  52 C              D LAMBDA(0)/DS > 0  .
  53 C
  54 C
  55 C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER
  56 C MUST SUPPLY A SUBROUTINE  F(X,V)  WHICH EVALUATES F(X) AT X
  57 C AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE
  58 C  FJACS(X,QR,LENQR,PIVOT)  WHICH EVALUATES THE (SYMMETRIC)
  59 C JACOBIAN MATRIX OF F(X) AT X, AND RETURNS THE SYMMETRIC
  60 C JACOBIAN MATRIX IN PACKED SKYLINE STORAGE FORMAT IN QR.  LENQR
  61 C AND PIVOT DESCRIBE THE DATA STRUCTURE IN QR.  FOR THE CURVE
  62 C TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE
  63 C  RHO(A,LAMBDA,X,V,PAR,IPAR)  WHICH EVALUATES THE HOMOTOPY MAP RHO
  64 C AT (A,X,LAMBDA) AND RETURNS THE VECTOR RHO(A,X,LAMBDA) IN V, AND
  65 C A SUBROUTINE  RHOJS(A,LAMBDA,X,QR,LENQR,PIVOT,PP,PAR,IPAR)  WHICH
  66 C RETURNS IN QR THE SYMMETRIC N X N JACOBIAN MATRIX [D RHO/DX]
  67 C EVALUATED AT (A,X,LAMBDA) AND STORED IN PACKED SKYLINE FORMAT, AND
  68 C RETURNS IN PP THE VECTOR -(D RHO/D LAMBDA) EVALUATED AT
  69 C (A,X,LAMBDA).  LENQR AND PIVOT DESCRIBE THE DATA STRUCTURE
  70 C IN QR.
  71 C *** NOTE THE MINUS SIGN IN THE DEFINITION OF PP. ***
  72 C
  73 C
  74 C FUNCTIONS AND SUBROUTINES DIRECTLY OR INDIRECTLY CALLED BY FIXPDS:
  75 C  D1MACH , F (OR  RHO ), FJACS (OR  RHOJS ), GMFADS , MFACDS ,
  76 C  MULTDS , PCGDS , PCGNS , QIMUDS , ROOT , ROOTNS , SOLVDS ,
  77 C  STEPNS , TANGNS , AND THE BLAS FUNCTIONS  DAXPY , DCOPY , DDOT ,
  78 C  DNRM2 , DSCAL , IDAMAX .  ONLY  D1MACH  CONTAINS MACHINE DEPENDENT
  79 C CONSTANTS.  NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED.
  80 C
  81 C
  82 C ON INPUT:
  83 C
  84 C N  IS THE DIMENSION OF X, F(X), AND RHO(A,X,LAMBDA).
  85 C
  86 C Y  IS AN ARRRAY OF LENGTH  N + 1.  (Y(1),...,Y(N)) = A  IS THE
  87 C    STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND
  88 C    ZERO FINDING PROBLEMS.  (Y(1),...,Y(N)) = X0  FOR THE CURVE
  89 C    TRACKING PROBLEM.
  90 C
  91 C IFLAG  CAN BE -2, -1, 0, 2, OR 3.  IFLAG  SHOULD BE 0 ON THE
  92 C    FIRST CALL TO  FIXPNS  FOR THE PROBLEM  X=F(X), -1 FOR THE
  93 C    PROBLEM  F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,X,LAMBDA)=0.
  94 C    IN CERTAIN SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPNS,
  95 C    AND  FIXPNS  CAN BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
  96 C
  97 C ARCRE , ARCAE  ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY,
  98 C    ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE.  IF
  99 C    ARC?E .LE. 0.0  ON INPUT IT IS RESET TO  .5*SQRT(ANS?E) .
 100 C    NORMALLY  ARC?E SHOULD BE CONSIDERABLY LARGER THAN  ANS?E .
 101 C
 102 C ANSRE , ANSAE  ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR
 103 C    THE ANSWER AT LAMBDA = 1.  THE ACCEPTED ANSWER  Y = (X, LAMBDA)
 104 C    SATISFIES
 105 C
 106 C       |Y(NP1) - 1|  .LE.  ANSRE + ANSAE           .AND.
 107 C
 108 C       ||Z||  .LE.  ANSRE*||X|| + ANSAE          WHERE
 109 C
 110 C    (Z,.) IS THE NEWTON STEP TO Y.
 111 C
 112 C TRACE  IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR
 113 C    INTERMEDIATE OUTPUT.  IF  TRACE .GT. 0  THE POINTS COMPUTED ON
 114 C    THE ZERO CURVE ARE WRITTEN TO I/O UNIT  TRACE .
 115 C
 116 C A(1:*)  CONTAINS THE PARAMETER VECTOR  A .  FOR THE FIXED POINT
 117 C    AND ZERO FINDING PROBLEMS, A  NEED NOT BE INITIALIZED BY THE
 118 C    USER, AND IS ASSUMED TO HAVE LENGTH  N.  FOR THE CURVE
 119 C    TRACKING PROBLEM, A  MUST BE INITIALIZED BY THE USER.
 120 C
 121 C YP(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO
 122 C    THE ZERO CURVE AT THE CURRENT POINT  Y .
 123 C
 124 C YOLD(1:N+1)  IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND
 125 C    ON THE ZERO CURVE.
 126 C
 127 C YPOLD(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO
 128 C    THE ZERO CURVE AT  YOLD .
 129 C
 130 C QR(1:LENQR)  IS A WORK ARRAY CONTAINING THE N X N SYMMETRIC
 131 C    JACOBIAN MATRIX WITH RESPECT TO X STORED IN PACKED SKYLINE
 132 C    STORAGE FORMAT.  LENQR  AND  PIVOT  DESCRIBE THE DATA
 133 C    STRUCTURE IN  QR .
 134 C
 135 C LENQR  IS THE LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR  .
 136 C
 137 C PIVOT(1:N+2)  IS A WORK ARRAY CONTAINING THE INDICES OF THE
 138 C    DIAGONAL ELEMENTS OF THE N X N SYMMETRIC JACOBIAN MATRIX
 139 C    (WITH RESPECT TO X) WITHIN  QR .
 140 C
 141 C WORK(1:13*(N+1)+2*N+LENQR)  IS A WORK ARRAY SPLIT UP AND USED
 142 C    FOR THE CALCULATION OF THE JACOBIAN MATRIX KERNEL, THE
 143 C    NEWTON STEP, INTERPOLATION, AND THE ESTIMATION OF THE OPTIMAL
 144 C    STEP SIZE  H .
 145 C
 146 C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P)  IS
 147 C    A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION.
 148 C    IF  SSPAR(J) .LE. 0.0  ON INPUT, IT IS RESET TO A DEFAULT VALUE
 149 C    BY  FIXPNS .  OTHERWISE THE INPUT VALUE OF  SSPAR(J)  IS USED.
 150 C    SEE THE COMMENTS BELOW AND IN  STEPNS  FOR MORE INFORMATION ABOUT
 151 C    THESE CONSTANTS.
 152 C
 153 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
 154 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
 155 C    RHO, RHOJS.
 156 C
 157 C
 158 C ON OUTPUT:
 159 C
 160 C N , TRACE , A  ARE UNCHANGED.
 161 C
 162 C (Y(1),...,Y(N)) = X, Y(NP1) = LAMBDA, AND Y IS AN APPROXIMATE
 163 C    ZERO OF THE HOMOTOPY MAP.  NORMALLY LAMBDA = 1 AND X IS A
 164 C    FIXED POINT(ZERO) OF F(X).  IN ABNORMAL SITUATIONS LAMBDA
 165 C    MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO).
 166 C
 167 C IFLAG =
 168 C  -2   CAUSES  FIXPNS  TO INITIALIZE EVERYTHING FOR THE PROBLEM
 169 C       RHO(A,X,LAMBDA) = 0 (USE ON FIRST CALL).
 170 C
 171 C  -1   CAUSES  FIXPNS  TO INITIALIZE EVERYTHING FOR THE PROBLEM
 172 C       F(X) = 0 (USE ON FIRST CALL).
 173 C
 174 C   0   CAUSES  FIXPNS  TO INITIALIZE EVERYTHING FOR THE PROBLEM
 175 C       X = F(X) (USE ON FIRST CALL).
 176 C
 177 C   1   NORMAL RETURN.
 178 C
 179 C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  SOME OR ALL OF
 180 C       ARCRE , ARCAE , ANSRE , ANSAE  HAVE BEEN INCREASED TO
 181 C       SUITABLE VALUES.  TO CONTINUE, JUST CALL  FIXPNS  AGAIN
 182 C       WITHOUT CHANGING ANY PARAMETERS.
 183 C
 184 C   3   STEPNS  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
 185 C       FIXPNS  AGAIN WITHOUT CHANGING ANY PARAMETERS.
 186 C
 187 C   4   THE PRECONDITIONED CONJUGATE GRADIENT ITERATION FAILED TO
 188 C       CONVERGE (PROBABLY BECAUSE THE JACOBIAN MATRIX DID NOT HAVE
 189 C       FULL RANK).  THE ALGORITHM HAS FAILED (THE ZERO CURVE OF
 190 C       THE HOMOTOPY MAP CANNOT BE FOLLOWED ANY FURTHER).
 191 C
 192 C   5   THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE
 193 C       HOMOTOPY MAP AND IS NOT MAKING PROGRESS.  THE ERROR TOLERANCES
 194 C       ARC?E  AND  ANS?E  WERE TOO LENIENT.  THE PROBLEM SHOULD BE
 195 C       RESTARTED BY CALLING  FIXPNS  WITH SMALLER ERROR TOLERANCES
 196 C       AND  IFLAG = 0 (-1, -2).
 197 C
 198 C   6   THE NORMAL FLOW NEWTON ITERATION IN  STEPNS  OR  ROOTNS
 199 C       FAILED TO CONVERGE.  THE ERROR TOLERANCES  ANS?E  MAY BE TOO
 200 C       STRINGENT.
 201 C
 202 C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
 203 C
 204 C ARCRE , ARCAE , ANSRE , ANSAE  ARE UNCHANGED AFTER A NORMAL RETURN
 205 C    (IFLAG = 1).  THEY ARE INCREASED TO APPROPRIATE VALUES ON THE
 206 C    RETURN  IFLAG = 2 .
 207 C
 208 C NFE  IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF
 209 C    JACOBIAN EVALUATIONS).
 210 C
 211 C ARCLEN  IS THE LENGTH OF THE PATH FOLLOWED.
 212 C
 213 C
 214 C
 215       DOUBLE PRECISION ABSERR,ANSAE,ANSRE,ARCAE,ARCLEN,ARCRE,
 216      1   CURSW,CURTOL,D1MACH,DNRM2,H,HOLD,RELERR,S
 217       INTEGER IFLAG,IFLAGC,IPP,IRHO,ITANGW,ITER,ITZ,IW,IWP,
 218      1   IZ0,IZ1,JW,LENQR,LIMIT,LIMITD,N,NC,NFE,NFEC,NP1,TRACE
 219       LOGICAL START,CRASH
 220 C
 221 C ***** ARRAY DECLARATIONS. *****
 222 C
 223       DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N),
 224      $   QR(LENQR),WORK(13*(N+1)+2*N+LENQR),SSPAR(8),PAR(1)
 225       INTEGER PIVOT(N+2),IPAR(1)
 226 C
 227 C ***** END OF DIMENSIONAL INFORMATION. *****
 228 C
 229       SAVE
 230 C
 231 C LIMITD  IS AN UPPER BOUND ON THE NUMBER OF STEPS.  IT MAY BE
 232 C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT:
 233       PARAMETER (LIMITD=1000)
 234 C
 235 C SWITCH FROM THE TOLERANCE  ARC?E  TO THE (FINER) TOLERANCE  ANS?E  IF
 236 C THE CURVATURE OF ANY COMPONENT OF  Y  EXCEEDS  CURSW.
 237       PARAMETER (CURSW=10.0)
 238 C
 239 C
 240 C
 241 C :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :
 242       IF (N .LE. 0  .OR.  ANSRE .LE. 0.0  .OR.  ANSAE .LT. 0.0)
 243      $                                                     IFLAG=7
 244       IF (IFLAG .GE. -2  .AND.  IFLAG .LE. 0) GO TO 20
 245       IF (IFLAG .EQ. 2) GO TO 120
 246       IF (IFLAG .EQ. 3) GO TO 90
 247 C ONLY VALID INPUT FOR  IFLAG  IS -2, -1, 0, 2, 3.
 248       IFLAG=7
 249       RETURN
 250 C
 251 C *****  INITIALIZATION BLOCK.  *****
 252 C
 253 20    ARCLEN=0.0
 254       IF (ARCRE .LE. 0.0) ARCRE=.5*SQRT(ANSRE)
 255       IF (ARCAE .LE. 0.0) ARCAE=.5*SQRT(ANSAE)
 256       NC=N
 257       NFEC=0
 258       IFLAGC=IFLAG
 259       NP1=N+1
 260 C SET INDICES FOR SPLITTING UP WORK ARRAY.
 261       IPP=1
 262       IRHO=N+1
 263       IW=IRHO+N
 264       IWP=IW+NP1
 265       ITZ=IWP+NP1
 266       IZ0=ITZ+NP1
 267       IZ1=IZ0+NP1
 268       ITANGW=IZ1+NP1
 269 C SET INITIAL CONDITIONS FOR FIRST CALL TO  STEPNS .
 270       START=.TRUE.
 271       CRASH=.FALSE.
 272       HOLD=1.0
 273       H=.1
 274       S=0.0
 275       YPOLD(NP1)=1.0
 276       YP(NP1)=1.0
 277       Y(NP1)=0.0
 278       DO 40 JW=1,N
 279         YPOLD(JW)=0.0
 280         YP(JW)=0.0
 281 40    CONTINUE
 282       DO 50 JW=ITANGW,ITANGW+NP1+N
 283         WORK(JW)=0.0
 284 50    CONTINUE
 285 C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS.
 286 C LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE
 287 C DAVIDENKO FLOW AND Y THEIR LIMIT.
 288 C IDEAL CONTRACTION FACTOR:  ||Z[2] - Z[1]|| / ||Z[1] - Z[0]||
 289       IF (SSPAR(1) .LE. 0.0) SSPAR(1)= .5
 290 C IDEAL RESIDUAL FACTOR:  ||RHO(A, Z[1])|| / ||RHO(A, Z[0])||
 291       IF (SSPAR(2) .LE. 0.0) SSPAR(2)= .01
 292 C IDEAL DISTANCE FACTOR:  ||Z[1] - Y|| / ||Z[0] - Y||
 293       IF (SSPAR(3) .LE. 0.0) SSPAR(3)= .5
 294 C MINIMUM STEP SIZE  HMIN .
 295       IF (SSPAR(4) .LE. 0.0) SSPAR(4)= (SQRT(N+1.0)+4.0)*D1MACH(4)
 296 C MAXIMUM STEP SIZE  HMAX .
 297       IF (SSPAR(5) .LE. 0.0) SSPAR(5)= 1.0
 298 C MINIMUM STEP SIZE REDUCTION FACTOR  BMIN .
 299       IF (SSPAR(6) .LE. 0.0) SSPAR(6)= .1
 300 C MAXIMUM STEP SIZE EXPANSION FACTOR  BMAX .
 301       IF (SSPAR(7) .LE. 0.0) SSPAR(7)= 3.0
 302 C ASSUMED OPERATING ORDER  P .
 303       IF (SSPAR(8) .LE. 0.0) SSPAR(8)= 2.0
 304 C
 305 C LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
 306       IF (IFLAGC .GE. -1) THEN
 307         CALL DCOPY(N,Y,1,A,1)
 308       ENDIF
 309 90    LIMIT=LIMITD
 310 C
 311 C *****  END OF INITIALIZATION BLOCK.  *****
 312 C
 313 C
 314 C *****  MAIN LOOP.  *****
 315 C
 316 120   DO 400 ITER=1,LIMIT
 317       IF (Y(NP1) .LT. 0.0) THEN
 318         ARCLEN=S
 319         IFLAG=5
 320         RETURN
 321       ENDIF
 322 C
 323 C SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH
 324 C CURVATURE COMPONENTS.
 325 140   CURTOL=CURSW*HOLD
 326       RELERR=ARCRE
 327       ABSERR=ARCAE
 328       DO 160 JW=1,NP1
 329         IF (ABS(YP(JW)-YPOLD(JW)) .GT. CURTOL) THEN
 330           RELERR=ANSRE
 331           ABSERR=ANSAE
 332           GO TO 200
 333         ENDIF
 334 160   CONTINUE
 335 C
 336 C TAKE A STEP ALONG THE CURVE.
 337 200   CALL STEPNS(NC,NFEC,IFLAGC,START,CRASH,HOLD,H,RELERR,ABSERR,
 338      +     S,Y,YP,YOLD,YPOLD,A,QR,LENQR,PIVOT,WORK,SSPAR,PAR,IPAR)
 339 C PRINT LATEST POINT ON CURVE IF REQUESTED.
 340       IF (TRACE .GT. 0) THEN
 341         WRITE (TRACE,217) ITER,NFEC,S,Y(NP1),(Y(JW),JW=1,NC)
 342 217     FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X,
 343      $  'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4))
 344       ENDIF
 345       NFE=NFEC
 346 C CHECK IF THE STEP WAS SUCCESSFUL.
 347       IF (IFLAGC .GT. 0) THEN
 348         ARCLEN=S
 349         IFLAG=IFLAGC
 350         RETURN
 351       ENDIF
 352       IF (CRASH) THEN
 353 C RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
 354         IFLAG=2
 355 C CHANGE ERROR TOLERANCES.
 356         IF (ARCRE .LT. RELERR) ARCRE=RELERR
 357         IF (ANSRE .LT. RELERR) ANSRE=RELERR
 358         IF (ARCAE .LT. ABSERR) ARCAE=ABSERR
 359         IF (ANSAE .LT. ABSERR) ANSAE=ABSERR
 360 C CHANGE LIMIT ON NUMBER OF ITERATIONS.
 361         LIMIT=LIMIT-ITER
 362         RETURN
 363       ENDIF
 364 C
 365       IF (Y(NP1) .GE. 1.0) THEN
 366 C
 367 C USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE
 368 C ANSWER AT LAMBDA = 1.0 .
 369 C
 370 C SAVE  YOLD  FOR ARC LENGTH CALCULATION LATER.
 371         CALL DCOPY(NP1,YOLD,1,WORK(IZ0),1)
 372 C
 373         CALL ROOTNS(NC,NFEC,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,YPOLD,
 374      $              A,QR,LENQR,PIVOT,WORK,PAR,IPAR)
 375 C
 376         NFE=NFEC
 377         IFLAG=1
 378 C SET ERROR FLAG IF  ROOTNS  COULD NOT GET THE POINT ON THE ZERO
 379 C CURVE AT  LAMBDA = 1.0  .
 380         IF (IFLAGC .GT. 0) IFLAG=IFLAGC
 381 C CALCULATE FINAL ARC LENGTH.
 382         DO 290 JW=1,NP1
 383           WORK(JW)=Y(JW) - WORK(IZ0+JW-1)
 384 290     CONTINUE
 385         ARCLEN=S - HOLD + DNRM2(NP1,WORK,1)
 386         RETURN
 387       ENDIF
 388 C
 389 400   CONTINUE
 390 C
 391 C *****  END OF MAIN LOOP.  *****
 392 C
 393 C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS.
 394       IFLAG=3
 395       ARCLEN=S
 396       RETURN
 397 C
 398 C
 399       END