share/hompack/fortran/fixpnf.f

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