SUBROUTINE FIXPNS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A, $ NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,WORK,SSPAR, $ PAR,IPAR) C C SUBROUTINE FIXPNS 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)) 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,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)/DX , D RHO(A,X,LAMBDA)/D LAMBDA] = 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 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 RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO C AT (A,X,LAMBDA) AND RETURNS THE VECTOR RHO(A,X,LAMBDA) IN V, AND C A 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, AND C RETURNS IN PP THE VECTOR -(D RHO/D LAMBDA) EVALUATED AT C (A,X,LAMBDA). LENQR AND PIVOT DESCRIBE THE DATA STRUCTURE C 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 RHO ), FJACS (OR RHOJS ), GMFADS , MFACDS , C MULTDS , PCGDS , PCGNS , QIMUDS , ROOT , ROOTNS , SOLVDS , C STEPNS , TANGNS , AND THE BLAS FUNCTIONS DAXPY , DCOPY , DDOT , C DNRM2 , DSCAL , IDAMAX . ONLY D1MACH CONTAINS MACHINE DEPENDENT C CONSTANTS. 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,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 FIXPNS 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 FIXPNS, C AND FIXPNS 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 = (X, LAMBDA) C SATISFIES C C |Y(NP1) - 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:LENQR) IS A WORK ARRAY CONTAINING THE N X N SYMMETRIC C JACOBIAN MATRIX WITH RESPECT TO X STORED IN PACKED SKYLINE C STORAGE FORMAT. LENQR AND PIVOT DESCRIBE THE DATA C STRUCTURE IN QR . C C LENQR IS THE LENGTH OF THE ONE-DIMENSIONAL ARRAY QR . C C PIVOT(1:N+2) IS A WORK ARRAY CONTAINING THE INDICES OF THE C DIAGONAL ELEMENTS OF THE N X N SYMMETRIC JACOBIAN MATRIX C (WITH RESPECT TO X) WITHIN QR . C C WORK(1:13*(N+1)+2*N+LENQR) IS A WORK ARRAY SPLIT UP AND USED C FOR THE CALCULATION OF THE JACOBIAN MATRIX KERNEL, THE C NEWTON STEP, INTERPOLATION, AND THE ESTIMATION OF THE OPTIMAL C STEP SIZE H . 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 FIXPNS . OTHERWISE THE INPUT VALUE OF SSPAR(J) IS USED. C SEE THE COMMENTS BELOW AND IN STEPNS 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, RHOJS. C C C ON OUTPUT: C C N , TRACE , A ARE UNCHANGED. C C (Y(1),...,Y(N)) = X, Y(NP1) = 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 FIXPNS TO INITIALIZE EVERYTHING FOR THE PROBLEM C RHO(A,X,LAMBDA) = 0 (USE ON FIRST CALL). C C -1 CAUSES FIXPNS TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPNS 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 FIXPNS AGAIN C WITHOUT CHANGING ANY PARAMETERS. C C 3 STEPNS HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPNS AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 THE PRECONDITIONED CONJUGATE GRADIENT ITERATION FAILED TO C CONVERGE (PROBABLY BECAUSE THE JACOBIAN MATRIX DID NOT HAVE C FULL RANK). THE ALGORITHM HAS FAILED (THE ZERO CURVE OF C THE HOMOTOPY MAP CANNOT BE 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 FIXPNS WITH SMALLER ERROR TOLERANCES C AND IFLAG = 0 (-1, -2). C C 6 THE NORMAL FLOW NEWTON ITERATION IN STEPNS OR ROOTNS 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,IPP,IRHO,ITANGW,ITER,ITZ,IW,IWP, 1 IZ0,IZ1,JW,LENQR,LIMIT,LIMITD,N,NC,NFE,NFEC,NP1,TRACE LOGICAL START,CRASH C C ***** ARRAY DECLARATIONS. ***** C DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N), $ QR(LENQR),WORK(13*(N+1)+2*N+LENQR),SSPAR(8),PAR(1) INTEGER PIVOT(N+2),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 : : : : : : : : : : : : : : : : : : : : : : : : 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 INDICES FOR SPLITTING UP WORK ARRAY. IPP=1 IRHO=N+1 IW=IRHO+N IWP=IW+NP1 ITZ=IWP+NP1 IZ0=ITZ+NP1 IZ1=IZ0+NP1 ITANGW=IZ1+NP1 C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNS . START=.TRUE. CRASH=.FALSE. HOLD=1.0 H=.1 S=0.0 YPOLD(NP1)=1.0 YP(NP1)=1.0 Y(NP1)=0.0 DO 40 JW=1,N YPOLD(JW)=0.0 YP(JW)=0.0 40 CONTINUE DO 50 JW=ITANGW,ITANGW+NP1+N WORK(JW)=0.0 50 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 CALL DCOPY(N,Y,1,A,1) ENDIF 90 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 120 DO 400 ITER=1,LIMIT IF (Y(NP1) .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 STEPNS(NC,NFEC,IFLAGC,START,CRASH,HOLD,H,RELERR,ABSERR, + S,Y,YP,YOLD,YPOLD,A,QR,LENQR,PIVOT,WORK,SSPAR,PAR,IPAR) C PRINT LATEST POINT ON CURVE IF REQUESTED. IF (TRACE .GT. 0) THEN WRITE (TRACE,217) ITER,NFEC,S,Y(NP1),(Y(JW),JW=1,NC) 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(NP1) .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. CALL DCOPY(NP1,YOLD,1,WORK(IZ0),1) C CALL ROOTNS(NC,NFEC,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,YPOLD, $ A,QR,LENQR,PIVOT,WORK,PAR,IPAR) C NFE=NFEC IFLAG=1 C SET ERROR FLAG IF ROOTNS 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 WORK(JW)=Y(JW) - WORK(IZ0+JW-1) 290 CONTINUE ARCLEN=S - HOLD + DNRM2(NP1,WORK,1) RETURN 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 C END