SUBROUTINE FIXPDS(N,Y,IFLAG,ARCTOL,EPS,TRACE,A,NDIMA,NFE, $ ARCLEN,YP,YPOLD,QR,LENQR,PIVOT,PP,WORK,WT,PHI,P, $ PAR,IPAR) C C SUBROUTINE FIXPDS 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,X,LAMBDA). 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) = (X(S), LAMBDA(S)). 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,X,LAMBDA) IS ASSUMED TO C BE A C2 MAP FROM E**M X E**N X [0,1) 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,X,LAMBDA)/D LAMBDA , D RHO(A,X,LAMBDA)/DX] = N C C FOR ALL POINTS (X,LAMBDA) SUCH THAT RHO(A,X,LAMBDA)=0. IT IS C FURTHER ASSUMED THAT C C RANK [ D RHO(A,X0,0)/DX ] = N . C C WITH A FIXED, THE ZERO CURVE OF RHO(A,X,LAMBDA) EMANATING C FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY C SOLVING THE ORDINARY DIFFERENTIAL EQUATION C D RHO(A,X(S),LAMBDA(S))/DS = 0 FOR Y(S) = (X(S), LAMBDA(S)), C WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY C MAP RHO(A,X,LAMBDA) IS ASSUMED TO BE CONSTRUCTED SUCH THAT C C D LAMBDA(0)/DS > 0 . C C THIS CODE IS BASED ON THE ALGORITHM IN L. T. WATSON, A C GLOBALLY CONVERGENT ALGORITHM FOR COMPUTING FIXED POINTS OF C C2 MAPS, APPL. MATH. COMPUT., 5 (1979) 297-311. C C C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER C MUST SUPPLY A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X C AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE C FJACS(X,QR,LENQR,PIVOT) WHICH EVALUATES THE (SYMMETRIC) C JACOBIAN MATRIX OF F(X) AT X, AND RETURNS THE SYMMETRIC C JACOBIAN MATRIX IN PACKED SKYLINE STORAGE FORMAT IN QR. LENQR C AND PIVOT DESCRIBE THE DATA STRUCTURE IN QR. FOR THE CURVE C TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE C RHOA(V,LAMBDA,X,PAR,IPAR) WHICH GIVEN (X,LAMBDA) RETURNS A C PARAMETER VECTOR A IN V SUCH THAT RHO(A,X,LAMBDA)=0, AND A C SUBROUTINE RHOJS(A,LAMBDA,X,QR,LENQR,PIVOT,PP,PAR,IPAR) WHICH C RETURNS IN QR THE SYMMETRIC N X N JACOBIAN MATRIX [D RHO/DX] C EVALUATED AT (A,X,LAMBDA) AND STORED IN PACKED SKYLINE FORMAT, C AND RETURNS IN PP THE VECTOR -(D RHO/D LAMBDA) EVALUATED AT C (A,X,LAMBDA). LENQR AND PIVOT DESCRIBE THE DATA STRUCTURE IN QR. C *** NOTE THE MINUS SIGN IN THE DEFINITION OF PP. *** C C C FUNCTIONS AND SUBROUTINES DIRECTLY OR INDIRECTLY CALLED BY FIXPDS: C D1MACH , F (OR RHOA ), FJACS (OR RHOJS ), FODEDS , GMFADS , C MFACDS , MULTDS , PCGDS , QIMUDS , ROOT , SINTRP , SOLVDS , C STEPDS , AND THE BLAS FUNCTIONS DAXPY , DCOPY , DDOT , DNRM2 , C DSCAL , IDAMAX . ONLY D1MACH CONTAINS MACHINE DEPENDENT C CONSTANTS. NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. C C ***WARNING: THIS SUBROUTINE IS GENERALLY MORE ROBUST THAN FIXPNS C AND FIXPQS , BUT MAY BE SLOWER THAN THOSE SUBROUTINES BY A C FACTOR OF TWO. C C C ON INPUT: C C N IS THE DIMENSION OF X, F(X), AND RHO(A,X,LAMBDA). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(1),...,Y(N)) = A IS THE C STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND C ZERO FINDING PROBLEMS. (Y(1),...,Y(N)) = 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 FIXPDS FOR THE PROBLEM X=F(X), -1 FOR THE C PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,X,LAMBDA)=0. C IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPDS, C AND FIXPDS CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. C C ARCTOL IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN C FOLLOWING THE ZERO CURVE. IF ARCTOL .LE. 0.0 ON INPUT C IT IS RESET TO .5*DSQRT(EPS). NORMALLY ARCTOL SHOULD C BE CONSIDERABLY LARGER THAN EPS. C C EPS IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN VERY C NEAR THE FIXED POINT(ZERO). EPS IS APPROXIMATELY THE C MIXED ABSOLUTE AND RELATIVE ERROR IN THE COMPUTED FIXED C POINT(ZERO). 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:NDIMA) 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 HAS LENGTH NDIMA AND MUST BE INITIALIZED C BY THE USER. C C NDIMA IS THE DIMENSION OF A , AND IS USED ONLY FOR THE CURVE C TRACKING PROBLEM. C C YP(1:N+1) IS A WORK ARRAY CONTAINING THE CURRENT TANGENT C VECTOR TO THE ZERO CURVE. C C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS TANGENT C VECTOR TO THE ZERO CURVE. C C QR(1:LENQR) IS A WORK ARRAY CONTAINING THE (SYMMETRIC) JACOBIAN C MATRIX WITH RESPECT TO X, IN THE PACKED SKYLINE STORAGE FORMAT. C C LENQR IS THE DIMENSION OF QR . C C PIVOT(1:N+2), PP(1:N), AND WORK(1:6*(N+1)+LENQR) ARE ALL WORK C ARRAYS USED BY FODEDS TO CALCULATE THE TANGENT VECTOR YP. C C WT(1:N+1), PHI(1:N+1,1:16), AND P(1:N+1) ARE ALL WORK ARRAYS C USED BY THE ODE SUBROUTINE STEPDS . 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 RHOA, RHOJS. C C Y, ARCTOL, EPS, ARCLEN, NFE, AND IFLAG SHOULD ALL BE C VARIABLES IN THE CALLING PROGRAM. C C C ON OUTPUT: C C N AND TRACE ARE UNCHANGED. C C (Y(1),...,Y(N)) = X, Y(N+1) = LAMBDA, 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 FIXPDS TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,X,LAMBDA) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPDS TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPDS 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. EPS HAS BEEN C INCREASED TO A SUITABLE VALUE. TO CONTINUE, JUST CALL C FIXPDS AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPDS HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPDS AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK AND/OR THE CONJUGATE C GRADIENT ITERATION FOR THE KERNEL OF THE JACOBIAN MATRIX C FAILED TO CONVERGE. THE ALGORITHM HAS FAILED (THE ZERO C CURVE OF THE HOMOTOPY MAP CANNOT BE FOLLOWED ANY FURTHER). C C 5 EPS (OR ARCTOL ) IS TOO LARGE. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPDS WITH A SMALLER EPS (OR C ARCTOL ) AND IFLAG = 0 (-1, -2). C C 6 I - DF(X) IS NEARLY SINGULAR AT THE FIXED POINT (DF(X) IS C NEARLY SINGULAR AT THE ZERO, OR D RHO(A,X,LAMBDA)/DX IS C NEARLY SINGULAR AT LAMBDA = 1 ). ANSWER MAY NOT BE C ACCURATE. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ARCTOL = EPS AFTER A NORMAL RETURN (IFLAG = 1). C C EPS IS UNCHANGED AFTER A NORMAL RETURN (IFLAG = 1). IT IS C INCREASED TO AN APPROPRIATE VALUE ON THE RETURN IFLAG = 2. C C A WILL (NORMALLY) HAVE BEEN MODIFIED. 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 DOUBLE PRECISION AOLD,ARCLEN,ARCTOL,CURSW,CURTOL,EPS, 1 EPSSTP,EPST,H,HOLD,S,S99,SA,SB,SOUT,SQNP1,XOLD,Y1SOUT INTEGER IFLAG,IFLAGC,ITER,IVC,J,JW,K,KGI,KOLD, 1 KSTEPS,LCODE,LENQR,LIMIT,LIMITD,N,NDIMA,NFE,NFEC,NP1,TRACE LOGICAL START,CRASH,ST99 C C ***** ARRAY DECLARATIONS. ***** C C ARRAYS NEEDED BY THE ODE SUBROUTINE STEPDS . DOUBLE PRECISION Y(N+1),WT(N+1),PHI(N+1,16),P(N+1),YP(N+1), 1 ALPHAS(12),W(12),G(13),GI(11) INTEGER IV(10) C C ARRAYS NEEDED BY FIXPDS , FODEDS , AND PCGDS . DOUBLE PRECISION YPOLD(N+1),A(N),QR(LENQR),PP(N), 1 WORK(6*(N+1)+LENQR),PAR(1) INTEGER PIVOT(N+2),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C SAVE EXTERNAL FODEDS 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 C : : : : : : : : : : : : : : : : : : : : : IF (N .LE. 0 .OR. EPS .LE. 0.0 ) IFLAG=7 IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10 IF (IFLAG .EQ. 2) GO TO 35 IF (IFLAG .EQ. 3) GO TO 30 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 10 ARCLEN=0.0 S=0.0 IF (ARCTOL .LE. 0.0) ARCTOL=.5*SQRT(EPS) NFEC=0 IFLAGC=IFLAG NP1=N+1 SQNP1=SQRT(DBLE(NP1)) C C SWITCH FROM THE TOLERANCE ARCTOL TO THE (FINER) TOLERANCE EPS IF C THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. C CURSW=10.0 C ST99=.FALSE. START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 EPSSTP=ARCTOL KSTEPS=0 C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION. YPOLD(NP1)=1.0 YP(NP1)=1.0 Y(NP1)=0.0 WORK(2*NP1)=0.0 WORK(3*NP1)=0.0 DO 20 J=1,N YPOLD(J)=0.0 YP(J)=0.0 WORK(NP1+J)=0.0 WORK(2*NP1+J)=0.0 20 CONTINUE C LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. IF (IFLAGC .GE. -1) THEN CALL DCOPY(N,Y,1,A,1) ENDIF 30 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 35 DO 150 ITER=1,LIMIT IF (Y(NP1) .LT. 0.0) THEN 40 ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 50 IF (S .LE. 7.0*SQNP1) GO TO 80 C ARC LENGTH IS GETTING TOO LONG, THE PROBLEM WILL BE C RESTARTED WITH A DIFFERENT A VECTOR. ARCLEN=ARCLEN+S S=0.0 60 START=.TRUE. CRASH=.FALSE. C COMPUTE A NEW A VECTOR. IF (IFLAGC .EQ. -2) THEN DO 63 JW=1,NDIMA QR(JW)=A(JW) 63 CONTINUE CALL RHOA(A,Y(NP1),Y,PAR,IPAR) DO 65 JW=1,NDIMA AOLD=QR(JW) C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE. IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 65 CONTINUE ELSE CALL F(Y,YP) DO 70 JW=1,N AOLD=A(JW) IF (IFLAGC .EQ. -1) THEN A(JW)=Y(NP1)*YP(JW)/(1.0 - Y(NP1)) + Y(JW) ELSE A(JW)=(Y(JW) - Y(NP1)*YP(JW))/(1.0 - Y(NP1)) ENDIF C IF NEW AND OLD A DIFFER BY TOO MUCH, TRACKING SHOULD NOT CONTINUE. IF (ABS(A(JW)-AOLD) .GT. 1.0+ABS(AOLD)) THEN ARCLEN=ARCLEN+S IFLAG=5 RETURN ENDIF 70 CONTINUE ENDIF GO TO 100 80 IF (Y(NP1) .LE. .99 .OR. ST99) GO TO 100 C WHEN LAMBDA REACHES .99, THE PROBLEM WILL BE RESTARTED WITH C A NEW A VECTOR. 90 ST99=.TRUE. EPSSTP=EPS ARCTOL=EPS GO TO 60 C C SET DIFFERENT ERROR TOLERANCE FOR HIGH CURVATURE COMPONENTS OF THE C TRAJECTORY Y(S). 100 CURTOL=CURSW*HOLD EPST=EPS/EPSSTP DO 110 JW=1,NP1 IF (ABS(YP(JW)-YPOLD(JW)) .LE. CURTOL) THEN WT(JW)=(ABS(Y(JW))+1.0) ELSE WT(JW)=(ABS(Y(JW))+1.0)*EPST ENDIF 110 CONTINUE C C TAKE A STEP ALONG THE CURVE. CALL STEPDS(FODEDS,NP1,Y,S,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH, + PHI,P,YP,ALPHAS,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI, + YPOLD,A,QR,LENQR,PIVOT,PP,WORK,NFEC,IFLAGC,PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,117) ITER,NFEC,S,Y(NP1),(Y(JW),JW=1,N) 117 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 .EQ. 4) THEN ARCLEN=ARCLEN+S IFLAG=4 RETURN ENDIF 120 IF (CRASH) THEN C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. EPS=EPSSTP IF (ARCTOL .LT. EPS) ARCTOL=EPS C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMIT=LIMIT-ITER RETURN ENDIF C 130 IF (Y(NP1) .GE. 1.0) THEN IF (ST99) GO TO 160 C C IF LAMBDA .GE. 1.0 BUT THE PROBLEM HAS NOT BEEN RESTARTED C WITH A NEW A VECTOR, BACK UP AND RESTART. C S99=S-.5*HOLD C GET AN APPROXIMATE ZERO Y(S) WITH Y(NP1)=LAMBDA .LT. 1.0 . 135 CALL SINTRP(S,Y,S99,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) IF (WT(NP1) .LT. 1.0) GO TO 140 S99=.5*(S-HOLD+S99) GO TO 135 C 140 CALL DCOPY(NP1,WT,1,Y,1) CALL DCOPY(NP1,YP,1,YPOLD,1) S=S99 GO TO 90 ENDIF C 150 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 RETURN C C C USE INVERSE INTERPOLATION TO GET THE ANSWER AT LAMBDA = 1.0 . C 160 SA=S-HOLD SB=S LCODE=1 170 CALL ROOT(SOUT,Y1SOUT,SA,SB,EPS,EPS,LCODE) C ROOT FINDS S SUCH THAT Y(NP1)(S) = LAMBDA = 1 . IF (LCODE .GT. 0) GO TO 190 CALL SINTRP(S,Y,SOUT,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) Y1SOUT=WT(NP1)-1.0 GO TO 170 190 IFLAG=1 C SET IFLAG = 6 IF ROOT COULD NOT GET LAMBDA = 1.0 . IF (LCODE .GT. 2) IFLAG=6 ARCLEN=ARCLEN+SA C LAMBDA(SA) = 1.0 . CALL SINTRP(S,Y,SA,WT,YP,NP1,KOLD,PHI,IVC,IV,KGI,GI, $ ALPHAS,G,W,XOLD,P) C CALL DCOPY(NP1,WT,1,Y,1) C RETURN END