SUBROUTINE FIXPNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,NFE, & ARCLEN,YP,YOLD,YPOLD,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR, & PAR,IPAR) C C SUBROUTINE FIXPNF FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE C OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X). FOR THE FIXED C POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL C INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (LAMBDA(S), X(S)) USING A HERMITE CUBIC PREDICTOR AND A C CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL C TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF C D(HOMOTOPY MAP)/DS ). C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE C OF THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO C BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR C ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET C OF E**M SATISFIES C C RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N C C FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,0,X0)/DX ] = N . C C WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING C FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY C SOLVING THE ORDINARY DIFFERENTIAL EQUATION C D RHO(A,LAMBDA(S),X(S))/DS = 0 FOR Y(S) = (LAMBDA(S), X(S)), C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY C MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0 . C C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY C A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X AND RETURNS THE C VECTOR F(X) IN V, AND A SUBROUTINE FJAC(X,V,K) WHICH RETURNS IN V C THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X. FOR C THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE C RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO AT C (A,LAMBDA,X) AND RETURNS THE VECTOR RHO(A,LAMBDA,X) IN V, AND A C SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR) WHICH RETURNS IN V THE KTH C COLUMN OF THE N X (N+1) JACOBIAN MATRIX [D RHO/D LAMBDA, D RHO/DX] C EVALUATED AT (A,LAMBDA,X). FIXPNF DIRECTLY OR INDIRECTLY USES C THE SUBROUTINES STEPNF , TANGNF , ROOTNF , ROOT , F (OR RHO ), C FJAC (OR RHOJAC ), D1MACH , AND THE BLAS FUNCTIONS DDOT AND C DNRM2 . ONLY D1MACH CONTAINS MACHINE DEPENDENT CONSTANTS. C NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(2),...,Y(N+1)) = A IS THE C STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND C ZERO FINDING PROBLEMS. (Y(2),...,Y(N+1)) = X0 FOR THE CURVE C TRACKING PROBLEM. C C IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE C FIRST CALL TO FIXPNF FOR THE PROBLEM X=F(X), -1 FOR THE C PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,LAMBDA,X)=0. C IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPNF, C AND FIXPNF CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCRE , ARCAE ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, C ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE. IF C ARC?E .LE. 0.0 ON INPUT IT IS RESET TO .5*SQRT(ANS?E) . C NORMALLY ARC?E SHOULD BE CONSIDERABLY LARGER THAN ANS?E . C C ANSRE , ANSAE ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR C THE ANSWER AT LAMBDA = 1. THE ACCEPTED ANSWER Y = (LAMBDA, X) C SATISFIES C C |Y(1) - 1| .LE. ANSRE + ANSAE .AND. C C ||Z|| .LE. ANSRE*||X|| + ANSAE WHERE C C (.,Z) IS THE NEWTON STEP TO Y. C C TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR C INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON C THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . C C A(1:*) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT C AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE C USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE C TRACKING PROBLEM, A MUST BE INITIALIZED BY THE USER. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT THE CURRENT POINT Y . C C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND C ON THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO C THE ZERO CURVE AT YOLD . C C QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) , W(1:N+1) , C WP(1:N+1) , Z0(1:N+1) , Z1(1:N+1) ARE ALL WORK ARRAYS USED BY C STEPNF TO CALCULATE THE TANGENT VECTORS AND NEWTON STEPS. 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 IF SSPAR(J) .LE. 0.0 ON INPUT, IT IS RESET TO A DEFAULT VALUE C BY FIXPNF . OTHERWISE THE INPUT VALUE OF SSPAR(J) IS USED. C SEE THE COMMENTS BELOW AND IN STEPNF FOR MORE INFORMATION ABOUT C THESE CONSTANTS. 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, RHOJAC. C C C ON OUTPUT: C C N , TRACE , A ARE UNCHANGED. C C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA C MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). C C IFLAG = C -2 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM C X = F(X) (USE ON FIRST CALL). C C 1 NORMAL RETURN. C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. SOME OR ALL OF C ARCRE , ARCAE , ANSRE , ANSAE HAVE BEEN INCREASED TO C SUITABLE VALUES. TO CONTINUE, JUST CALL FIXPNF AGAIN C WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPNF HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPNF AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM C HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE C FOLLOWED ANY FURTHER). C C 5 THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE C HOMOTOPY MAP AND IS NOT MAKING PROGRESS. THE ERROR TOLERANCES C ARC?E AND ANS?E WERE TOO LENIENT. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPNF WITH SMALLER ERROR TOLERANCES C AND IFLAG = 0 (-1, -2). C C 6 THE NORMAL FLOW NEWTON ITERATION IN STEPNF OR ROOTNF C FAILED TO CONVERGE. THE ERROR TOLERANCES ANS?E MAY BE TOO C STRINGENT. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCRE , ARCAE , ANSRE , ANSAE ARE UNCHANGED AFTER A NORMAL RETURN C (IFLAG = 1). THEY ARE INCREASED TO APPROPRIATE VALUES ON THE C RETURN IFLAG = 2 . C C NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF C JACOBIAN EVALUATIONS). C C ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. C C C DOUBLE PRECISION ABSERR,ANSAE,ANSRE,ARCAE,ARCLEN,ARCRE, 1 CURSW,CURTOL,D1MACH,DNRM2,H,HOLD,RELERR,S INTEGER IFLAG,IFLAGC,ITER,JW,LIMIT,LIMITD,N,NC,NFE,NFEC,NP1, 1 TRACE LOGICAL CRASH,POLSYS,START C C ***** ARRAY DECLARATIONS. ***** C DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N), & QR(N,N+2),ALPHA(N),TZ(N+1),W(N+1),WP(N+1),Z0(N+1), & Z1(N+1),SSPAR(8),PAR(1) INTEGER PIVOT(N+1),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C SAVE C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMITD=1000) C C SWITCH FROM THE TOLERANCE ARC?E TO THE (FINER) TOLERANCE ANS?E IF C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. PARAMETER (CURSW=10.0) C C C C : : : : : : : : : : : : : : : : : : : : : : : : C SET LOGICAL SWITCH TO REFLECT ENTRY POINT. POLSYS=.FALSE. GO TO 11 ENTRY POLYNF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,NFE, & ARCLEN,YP,YOLD,YPOLD,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR, & PAR,IPAR) POLSYS=.TRUE. 11 CONTINUE C IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0) & IFLAG=7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 20 IF (IFLAG .EQ. 2) GO TO 120 IF (IFLAG .EQ. 3) GO TO 90 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 20 ARCLEN=0.0 IF (ARCRE .LE. 0.0) ARCRE=.5*SQRT(ANSRE) IF (ARCAE .LE. 0.0) ARCAE=.5*SQRT(ANSAE) NC=N NFEC=0 IFLAGC=IFLAG NP1=N+1 C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNF . START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 S=0.0 YPOLD(1)=1.0 YP(1)=1.0 Y(1)=0.0 DO 40 JW=2,NP1 YPOLD(JW)=0.0 YP(JW)=0.0 40 CONTINUE C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. C LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE C DAVIDENKO FLOW AND Y THEIR LIMIT. C IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]|| IF (SSPAR(1) .LE. 0.0) SSPAR(1)= .5 C IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])|| IF (SSPAR(2) .LE. 0.0) SSPAR(2)= .01 C IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y|| IF (SSPAR(3) .LE. 0.0) SSPAR(3)= .5 C MINIMUM STEP SIZE HMIN . IF (SSPAR(4) .LE. 0.0) SSPAR(4)= (SQRT(N+1.0)+4.0)*D1MACH(4) C MAXIMUM STEP SIZE HMAX . IF (SSPAR(5) .LE. 0.0) SSPAR(5)= 1.0 C MINIMUM STEP SIZE REDUCTION FACTOR BMIN . IF (SSPAR(6) .LE. 0.0) SSPAR(6)= .1 C MAXIMUM STEP SIZE EXPANSION FACTOR BMAX . IF (SSPAR(7) .LE. 0.0) SSPAR(7)= 3.0 C ASSUMED OPERATING ORDER P . IF (SSPAR(8) .LE. 0.0) SSPAR(8)= 2.0 C C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. IF (IFLAGC .GE. -1) THEN DO 60 JW=2,NP1 A(JW-1)=Y(JW) 60 CONTINUE ENDIF 90 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 120 DO 400 ITER=1,LIMIT IF (Y(1) .LT. 0.0) THEN ARCLEN=S IFLAG=5 RETURN ENDIF C C SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH C CURVATURE COMPONENTS. 140 CURTOL=CURSW*HOLD RELERR=ARCRE ABSERR=ARCAE DO 160 JW=1,NP1 IF (ABS(YP(JW)-YPOLD(JW)) .GT. CURTOL) THEN RELERR=ANSRE ABSERR=ANSAE GO TO 200 ENDIF 160 CONTINUE C C TAKE A STEP ALONG THE CURVE. 200 CALL STEPNF(NC,NFEC,IFLAGC,START,CRASH,HOLD,H,RELERR,ABSERR, + S,Y,YP,YOLD,YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP,Z0,Z1,SSPAR, + PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,217) ITER,NFEC,S,Y(1),(Y(JW),JW=2,NP1) 217 FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X, & 'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4)) ENDIF NFE=NFEC C CHECK IF THE STEP WAS SUCCESSFUL. IF (IFLAGC .GT. 0) THEN ARCLEN=S IFLAG=IFLAGC RETURN ENDIF IF (CRASH) THEN C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. IF (ARCRE .LT. RELERR) ARCRE=RELERR IF (ANSRE .LT. RELERR) ANSRE=RELERR IF (ARCAE .LT. ABSERR) ARCAE=ABSERR IF (ANSAE .LT. ABSERR) ANSAE=ABSERR C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMIT=LIMIT-ITER RETURN ENDIF C IF (Y(1) .GE. 1.0) THEN C C USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE C ANSWER AT LAMBDA = 1.0 . C C SAVE YOLD FOR ARC LENGTH CALCULATION LATER. DO 260 JW=1,NP1 Z0(JW)=YOLD(JW) 260 CONTINUE CALL ROOTNF(NC,NFEC,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,YPOLD, & A,QR,ALPHA,TZ,PIVOT,W,WP,PAR,IPAR) C NFE=NFEC IFLAG=1 C SET ERROR FLAG IF ROOTNF COULD NOT GET THE POINT ON THE ZERO C CURVE AT LAMBDA = 1.0 . IF (IFLAGC .GT. 0) IFLAG=IFLAGC C CALCULATE FINAL ARC LENGTH. DO 290 JW=1,NP1 W(JW)=Y(JW) - Z0(JW) 290 CONTINUE ARCLEN=S - HOLD + DNRM2(NP1,W,1) RETURN ENDIF C C FOR POLYNOMIAL SYSTEMS AND THE POLSYS HOMOTOPY MAP, C D LAMBDA/DS .GE. 0 NECESSARILY. THIS CONDITION IS FORCED HERE IF C THE ENTRY POINT WAS POLYNF . C IF (POLSYS) THEN IF (YP(1) .LT. 0.0) THEN C REVERSE TANGENT DIRECTION SO D LAMBDA/DS = YP(1) > 0 . DO 310 JW=1,NP1 YP(JW)=-YP(JW) YPOLD(JW)=YP(JW) 310 CONTINUE C FORCE STEPNF TO USE THE LINEAR PREDICTOR FOR THE NEXT STEP ONLY. START=.TRUE. ENDIF ENDIF C 400 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 ARCLEN=S RETURN C END