share/hompack/fortran/fixpdf.f

   1       SUBROUTINE FIXPDF(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,NFE,
   2      $     ARCLEN,YP,YPOLD,QR,ALPHA,TZ,PIVOT,WT,PHI,P,PAR,IPAR)
   3 C
   4 C SUBROUTINE  FIXPDF  FINDS A FIXED POINT OR ZERO OF THE
   5 C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE
   6 C OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X).  FOR THE FIXED
   7 C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL
   8 C INTO ITSELF.  THE EQUATION  X = F(X)  IS SOLVED BY
   9 C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP
  10 C
  11 C  LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A)  ,
  12 C
  13 C STARTING FROM LAMBDA = 0, X = A.  THE CURVE IS PARAMETERIZED
  14 C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY
  15 C DIFFERENTIAL EQUATION  D(HOMOTOPY MAP)/DS = 0  FOR
  16 C Y(S) = (LAMBDA(S), X(S)).
  17 C
  18 C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP
  19 C SUCH THAT FOR SOME R > 0,  X*F(X) >= 0  WHENEVER NORM(X) = R.
  20 C THE EQUATION  F(X) = 0  IS SOLVED BY FOLLOWING THE ZERO CURVE
  21 C OF THE HOMOTOPY MAP
  22 C
  23 C   LAMBDA*F(X) + (1 - LAMBDA)*(X - A)
  24 C
  25 C EMANATING FROM LAMBDA = 0, X = A.
  26 C
  27 C  A  MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS.
  28 C
  29 C FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO
  30 C BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR
  31 C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET
  32 C OF E**M SATISFIES
  33 C
  34 C  RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N
  35 C
  36 C FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0.  IT IS
  37 C FURTHER ASSUMED THAT
  38 C
  39 C           RANK [ D RHO(A,0,X0)/DX ] = N  .
  40 C
  41 C WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING
  42 C FROM  LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY
  43 C SOLVING THE ORDINARY DIFFERENTIAL EQUATION
  44 C D RHO(A,LAMBDA(S),X(S))/DS = 0  FOR  Y(S) = (LAMBDA(S), X(S)),
  45 C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE.  ALSO THE HOMOTOPY
  46 C MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT
  47 C
  48 C              D LAMBDA(0)/DS > 0  .
  49 C
  50 C THIS CODE IS BASED ON THE ALGORITHM IN L. T. WATSON, A
  51 C GLOBALLY CONVERGENT ALGORITHM FOR COMPUTING FIXED POINTS OF
  52 C C2 MAPS, APPL. MATH. COMPUT., 5 (1979) 297-311.
  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  FJAC(X,V,K)
  58 C WHICH RETURNS IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF
  59 C F(X) EVALUATED AT X.  FOR THE CURVE TRACKING PROBLEM, THE USER MUST
  60 C SUPPLY A SUBROUTINE  RHOA(V,LAMBDA,X,PAR,IPAR)  WHICH GIVEN
  61 C (LAMBDA,X) RETURNS A PARAMETER VECTOR A IN V SUCH THAT
  62 C RHO(A,LAMBDA,X)=0, AND A SUBROUTINE  RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR)
  63 C WHICH RETURNS IN V THE KTH COLUMN OF THE N X (N+1) JACOBIAN
  64 C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT (A,LAMBDA,X).
  65 C  FIXPDF  DIRECTLY OR INDIRECTLY USES THE SUBROUTINES
  66 C  STEPS , SINTRP , ROOT , FODE , F (OR  RHOA ),
  67 C FJAC (OR  RHOJAC ), DCPOSE , D1MACH , AND THE BLAS FUNCTIONS
  68 C DDOT  AND  DNRM2 .  ONLY  D1MACH  CONTAINS MACHINE
  69 C DEPENDENT CONSTANTS.  NO OTHER MODIFICATIONS BY THE USER ARE
  70 C REQUIRED.
  71 C
  72 C ***WARNING:  THIS SUBROUTINE IS GENERALLY MORE ROBUST THAN  FIXPNF
  73 C AND  FIXPQF , BUT MAY BE SLOWER THAN THOSE SUBROUTINES BY A
  74 C FACTOR OF TWO.
  75 C
  76 C
  77 C ON INPUT:
  78 C
  79 C N  IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X).
  80 C
  81 C Y  IS AN ARRRAY OF LENGTH  N + 1.  (Y(2),...,Y(N+1)) = A  IS THE
  82 C    STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND
  83 C    ZERO FINDING PROBLEMS.  (Y(2),...,Y(N+1)) = X0  FOR THE CURVE
  84 C    TRACKING PROBLEM.
  85 C
  86 C IFLAG  CAN BE -2, -1, 0, 2, OR 3.  IFLAG  SHOULD BE 0 ON THE
  87 C    FIRST CALL TO  FIXPDF  FOR THE PROBLEM  X=F(X), -1 FOR THE
  88 C    PROBLEM  F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,LAMBDA,X)=0.
  89 C    IN CERTAIN SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPDF,
  90 C    AND  FIXPDF  CAN BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
  91 C
  92 C ARCTOL  IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN
  93 C    FOLLOWING THE ZERO CURVE.  IF  ARCTOL .LE. 0.0  ON INPUT
  94 C    IT IS RESET TO  .5*DSQRT(EPS).  NORMALLY  ARCTOL  SHOULD
  95 C    BE CONSIDERABLY LARGER THAN  EPS.
  96 C
  97 C EPS  IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN VERY
  98 C    NEAR THE FIXED POINT(ZERO).  EPS  IS APPROXIMATELY THE
  99 C    MIXED ABSOLUTE AND RELATIVE ERROR IN THE COMPUTED FIXED
 100 C    POINT(ZERO).
 101 C
 102 C TRACE  IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR
 103 C    INTERMEDIATE OUTPUT.  IF  TRACE .GT. 0  THE POINTS COMPUTED ON
 104 C    THE ZERO CURVE ARE WRITTEN TO I/O UNIT  TRACE .
 105 C
 106 C A(1:NDIMA) CONTAINS THE PARAMETER VECTOR  A .  FOR THE FIXED POINT
 107 C    AND ZERO FINDING PROBLEMS, A  NEED NOT BE INITIALIZED BY THE
 108 C    USER, AND IS ASSUMED TO HAVE LENGTH  N.  FOR THE CURVE
 109 C    TRACKING PROBLEM, A  HAS LENGTH  NDIMA  AND MUST BE INITIALIZED
 110 C    BY THE USER.
 111 C
 112 C NDIMA  IS THE DIMENSION OF  A, AND IS USED ONLY FOR THE CURVE
 113 C    TRACKING PROBLEM.
 114 C
 115 C YP(1:N+1) IS A WORK ARRAY CONTAINING THE CURRENT TANGENT
 116 C    VECTOR TO THE ZERO CURVE.
 117 C
 118 C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS TANGENT
 119 C    VECTOR TO THE ZERO CURVE.
 120 C
 121 C QR(1:N,1:N+1), ALPHA(1:N), TZ(1:N+1), AND PIVOT(1:N+1) ARE
 122 C    ALL WORK ARRAYS USED BY  FODE  TO CALCULATE THE TANGENT
 123 C    VECTOR YP.
 124 C
 125 C WT(1:N+1), PHI(1:N+1,1:16), AND P(1:N+1) ARE ALL WORK ARRAYS
 126 C    USED BY THE ODE SUBROUTINE  STEPS  .
 127 C
 128 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
 129 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
 130 C    RHOA, RHOJAC.
 131 C
 132 C Y, ARCTOL, EPS, ARCLEN, NFE, AND IFLAG SHOULD ALL BE
 133 C VARIABLES IN THE CALLING PROGRAM.
 134 C
 135 C
 136 C ON OUTPUT:
 137 C
 138 C N  AND  TRACE  ARE UNCHANGED.
 139 C
 140 C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE
 141 C    ZERO OF THE HOMOTOPY MAP.  NORMALLY LAMBDA = 1 AND X IS A
 142 C    FIXED POINT(ZERO) OF F(X).  IN ABNORMAL SITUATIONS LAMBDA
 143 C    MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO).
 144 C
 145 C IFLAG =
 146 C  -2   CAUSES  FIXPDF  TO INITIALIZE EVERYTHING FOR THE PROBLEM
 147 C       RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL).
 148 C
 149 C  -1   CAUSES  FIXPDF  TO INITIALIZE EVERYTHING FOR THE PROBLEM
 150 C       F(X) = 0 (USE ON FIRST CALL).
 151 C
 152 C   0   CAUSES  FIXPDF  TO INITIALIZE EVERYTHING FOR THE PROBLEM
 153 C       X = F(X) (USE ON FIRST CALL).
 154 C
 155 C   1   NORMAL RETURN.
 156 C
 157 C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  EPS HAS BEEN
 158 C       INCREASED TO A SUITABLE VALUE.  TO CONTINUE, JUST CALL
 159 C       FIXPDF  AGAIN WITHOUT CHANGING ANY PARAMETERS.
 160 C
 161 C   3   STEPS  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
 162 C       FIXPDF  AGAIN WITHOUT CHANGING ANY PARAMETERS.
 163 C
 164 C   4   JACOBIAN MATRIX DOES NOT HAVE FULL RANK.  THE ALGORITHM
 165 C       HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE
 166 C       FOLLOWED ANY FURTHER).
 167 C
 168 C   5   EPS  (OR  ARCTOL ) IS TOO LARGE.  THE PROBLEM SHOULD BE
 169 C       RESTARTED BY CALLING  FIXPDF  WITH A SMALLER  EPS  (OR
 170 C       ARCTOL ) AND  IFLAG = 0 (-1, -2).
 171 C
 172 C   6   I - DF(X)  IS NEARLY SINGULAR AT THE FIXED POINT (DF(X) IS
 173 C       NEARLY SINGULAR AT THE ZERO, OR  D RHO(A,LAMBDA,X)/DX  IS
 174 C       NEARLY SINGULAR AT  LAMBDA = 1 ).  ANSWER MAY NOT BE
 175 C       ACCURATE.
 176 C
 177 C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
 178 C
 179 C ARCTOL = EPS AFTER A NORMAL RETURN (IFLAG = 1).
 180 C
 181 C EPS  IS UNCHANGED AFTER A NORMAL RETURN (IFLAG = 1).  IT IS
 182 C    INCREASED TO AN APPROPRIATE VALUE ON THE RETURN IFLAG = 2.
 183 C
 184 C A  WILL (NORMALLY) HAVE BEEN MODIFIED.
 185 C
 186 C NFE  IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF
 187 C    JACOBIAN EVALUATIONS).
 188 C
 189 C ARCLEN  IS THE LENGTH OF THE PATH FOLLOWED.
 190 C
 191 C
 192       DOUBLE PRECISION AOLD,ARCLEN,ARCTOL,CURSW,CURTOL,EPS,
 193      1  EPSSTP,EPST,H,HOLD,S,S99,SA,SB,SOUT,SQNP1,XOLD,Y1SOUT
 194       INTEGER IFLAG,IFLAGC,ITER,IVC,J,JUDY,JW,K,KGI,KOLD,
 195      1  KSTEPS,LCODE,LIMIT,LIMITD,N,NDIMA,NFE,NFEC,NP1,TRACE
 196       LOGICAL START,CRASH,ST99
 197 C
 198 C *****  ARRAY DECLARATIONS.  *****
 199 C
 200 C ARRAYS NEEDED BY THE ODE SUBROUTINE  STEPS .
 201       DOUBLE PRECISION Y(N+1),WT(N+1),PHI(N+1,16),P(N+1),YP(N+1),
 202      1     ALPHAS(12),W(12),G(13),GI(11)
 203       INTEGER IV(10)
 204 C
 205 C ARRAYS NEEDED BY  FIXPDF , FODE , AND  DCPOSE .
 206       DOUBLE PRECISION YPOLD(N+1),A(N),QR(N,N+1),ALPHA(N),TZ(N+1),
 207      $     PAR(1)
 208       INTEGER PIVOT(N+1),IPAR(1)
 209 C
 210 C *****  END OF DIMENSIONAL INFORMATION.  *****
 211 C
 212       SAVE
 213       EXTERNAL FODE
 214 C
 215 C LIMITD  IS AN UPPER BOUND ON THE NUMBER OF STEPS.  IT MAY BE
 216 C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT:
 217       PARAMETER (LIMITD=1000)
 218 C
 219 C
 220 C :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :
 221       IF (N .LE. 0  .OR.  EPS .LE. 0.0 ) IFLAG=7
 222       IF (IFLAG .GE. -2  .AND.  IFLAG .LE. 0) GO TO 10
 223       IF (IFLAG .EQ. 2) GO TO 35
 224       IF (IFLAG .EQ. 3) GO TO 30
 225 C ONLY VALID INPUT FOR  IFLAG  IS -2, -1, 0, 2, 3.
 226       IFLAG=7
 227       RETURN
 228 C
 229 C *****  INITIALIZATION BLOCK.  *****
 230 C
 231 10    ARCLEN=0.0
 232       S=0.0
 233       IF (ARCTOL .LE. 0.0) ARCTOL=.5*SQRT(EPS)
 234       NFEC=0
 235       IFLAGC=IFLAG
 236       NP1=N+1
 237       SQNP1=SQRT(DBLE(NP1))
 238 C
 239 C SWITCH FROM THE TOLERANCE  ARCTOL  TO THE (FINER) TOLERANCE  EPS  IF
 240 C THE CURVATURE OF ANY COMPONENT OF  Y  EXCEEDS  CURSW.
 241 C
 242       CURSW=10.0
 243 C
 244       ST99=.FALSE.
 245       START=.TRUE.
 246       CRASH=.FALSE.
 247       HOLD=1.0
 248       H=.1
 249       EPSSTP=ARCTOL
 250       KSTEPS=0
 251 C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION.
 252       YPOLD(1)=1.0
 253       YP(1)=1.0
 254       Y(1)=0.0
 255       DO 20 J=2,NP1
 256         YPOLD(J)=0.0
 257         YP(J)=0.0
 258 20    CONTINUE
 259 C LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
 260       IF (IFLAGC .GE. -1) THEN
 261         DO 23 J=2,NP1
 262           A(J-1)=Y(J)
 263 23      CONTINUE
 264       ENDIF
 265 30    LIMIT=LIMITD
 266 C
 267 C *****  END OF INITIALIZATION BLOCK.  *****
 268 C
 269 C
 270 C *****  MAIN LOOP.  *****
 271 C
 272 35    DO 150 ITER=1,LIMIT
 273       IF (Y(1) .LT. 0.0) THEN
 274 40      ARCLEN=ARCLEN+S
 275         IFLAG=5
 276         RETURN
 277       ENDIF
 278 50    IF (S .LE. 7.0*SQNP1) GO TO 80
 279 C ARC LENGTH IS GETTING TOO LONG, THE PROBLEM WILL BE
 280 C RESTARTED WITH A DIFFERENT  A  VECTOR.
 281       ARCLEN=ARCLEN+S
 282       S=0.0
 283 60    START=.TRUE.
 284       CRASH=.FALSE.
 285 C COMPUTE A NEW  A  VECTOR.
 286       IF (IFLAGC .EQ. -2) THEN
 287         DO 63 JW=1,NDIMA
 288           QR(JW,1)=A(JW)
 289 63      CONTINUE
 290         CALL RHOA(A,Y(1),Y(2),PAR,IPAR)
 291         DO 65 JW=1,NDIMA
 292           AOLD=QR(JW,1)
 293 C IF NEW AND OLD  A  DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE.
 294           IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN
 295             ARCLEN=ARCLEN+S
 296             IFLAG=5
 297             RETURN
 298           ENDIF
 299 65      CONTINUE
 300       ELSE
 301         CALL F(Y(2),YP)
 302         DO 70 JW=1,N
 303           AOLD=A(JW)
 304           IF (IFLAGC .EQ. -1) THEN
 305             A(JW)=Y(1)*YP(JW)/(1.0 - Y(1)) + Y(JW+1)
 306           ELSE
 307             A(JW)=(Y(JW+1) - Y(1)*YP(JW))/(1.0 - Y(1))
 308           ENDIF
 309 C IF NEW AND OLD  A  DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE.
 310           IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN
 311             ARCLEN=ARCLEN+S
 312             IFLAG=5
 313             RETURN
 314           ENDIF
 315 70      CONTINUE
 316       ENDIF
 317       GO TO 100
 318 80    IF (Y(1) .LE. .99  .OR. ST99) GO TO 100
 319 C WHEN LAMBDA REACHES .99, THE PROBLEM WILL BE RESTARTED WITH
 320 C A NEW  A  VECTOR.
 321 90    ST99=.TRUE.
 322       EPSSTP=EPS
 323       ARCTOL=EPS
 324       GO TO 60
 325 C
 326 C SET DIFFERENT ERROR TOLERANCE FOR HIGH CURVATURE COMPONENTS OF THE
 327 C TRAJECTORY Y(S).
 328 100   CURTOL=CURSW*HOLD
 329       EPST=EPS/EPSSTP
 330       DO 110 JW=1,NP1
 331         IF (ABS(YP(JW)-YPOLD(JW)) .LE. CURTOL) THEN
 332           WT(JW)=(ABS(Y(JW))+1.0)
 333         ELSE
 334           WT(JW)=(ABS(Y(JW))+1.0)*EPST
 335         ENDIF
 336 110   CONTINUE
 337 C
 338 C TAKE A STEP ALONG THE CURVE.
 339       CALL STEPS(FODE,NP1,Y,S,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH,
 340      +     PHI,P,YP,ALPHAS,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI,
 341      +     YPOLD,A,QR,ALPHA,TZ,PIVOT,NFEC,IFLAGC,PAR,IPAR)
 342 C PRINT LATEST POINT ON CURVE IF REQUESTED.
 343       IF (TRACE .GT. 0) THEN
 344         WRITE (TRACE,117) ITER,NFEC,S,Y(1),(Y(JW),JW=2,NP1)
 345 117     FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X,
 346      $  'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4))
 347       ENDIF
 348       NFE=NFEC
 349 C CHECK IF THE STEP WAS SUCCESSFUL.
 350       IF (IFLAGC .EQ. 4) THEN
 351         ARCLEN=ARCLEN+S
 352         IFLAG=4
 353         RETURN
 354       ENDIF
 355 120   IF (CRASH) THEN
 356 C RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
 357         IFLAG=2
 358 C CHANGE ERROR TOLERANCES.
 359         EPS=EPSSTP
 360         IF (ARCTOL .LT. EPS) ARCTOL=EPS
 361 C CHANGE LIMIT ON NUMBER OF ITERATIONS.
 362         LIMIT=LIMIT-ITER
 363         RETURN
 364       ENDIF
 365 C
 366 130   IF (Y(1) .GE. 1.0) THEN
 367         IF (ST99) GO TO 160
 368 C
 369 C IF LAMBDA .GE. 1.0 BUT THE PROBLEM HAS NOT BEEN RESTARTED
 370 C WITH A NEW  A  VECTOR, BACK UP AND RESTART.
 371 C
 372         S99=S-.5*HOLD
 373 C GET AN APPROXIMATE ZERO Y(S) WITH  Y(1)=LAMBDA .LT. 1.0  .
 374 135     CALL SINTRP(S,Y,S99,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI,
 375      $     ALPHAS,G,W,XOLD,P)
 376         IF (WT(1) .LT. 1.0) GO TO 140
 377         S99=.5*(S-HOLD+S99)
 378         GO TO 135
 379 C
 380 140     DO 144 JUDY=1,NP1
 381           Y(JUDY)=WT(JUDY)
 382           YPOLD(JUDY)=YP(JUDY)
 383 144     CONTINUE
 384         S=S99
 385         GO TO 90
 386       ENDIF
 387 C
 388 150   CONTINUE
 389 C
 390 C *****  END OF MAIN LOOP.  *****
 391 C
 392 C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS.
 393       IFLAG=3
 394       RETURN
 395 C
 396 C
 397 C USE INVERSE INTERPOLATION TO GET THE ANSWER AT LAMBDA = 1.0 .
 398 C
 399 160   SA=S-HOLD
 400       SB=S
 401       LCODE=1
 402 170   CALL ROOT(SOUT,Y1SOUT,SA,SB,EPS,EPS,LCODE)
 403 C ROOT  FINDS S SUCH THAT Y(1)(S) = LAMBDA = 1 .
 404       IF (LCODE .GT. 0) GO TO 190
 405       CALL SINTRP(S,Y,SOUT,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI,
 406      $     ALPHAS,G,W,XOLD,P)
 407       Y1SOUT=WT(1)-1.0
 408       GO TO 170
 409 190   IFLAG=1
 410 C SET IFLAG = 6 IF  ROOT  COULD NOT GET  LAMBDA = 1.0  .
 411       IF (LCODE .GT. 2) IFLAG=6
 412       ARCLEN=ARCLEN+SA
 413 C LAMBDA(SA) = 1.0 .
 414       CALL SINTRP(S,Y,SA,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI,
 415      $     ALPHAS,G,W,XOLD,P)
 416 C
 417       DO 210 J=1,NP1
 418 210   Y(J)=WT(J)
 419       RETURN
 420       END