SUBROUTINE ROOTNS(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,YPOLD, $ A,QR,LENQR,PIVOT,WORK,PAR,IPAR) C C ROOTNS FINDS THE POINT YBAR = (XBAR, 1) ON THE ZERO CURVE OF THE C HOMOTOPY MAP. IT STARTS WITH TWO POINTS YOLD=(XOLD,LAMBDAOLD) AND C Y=(X,LAMBDA) SUCH THAT LAMBDAOLD < 1 <= LAMBDA , AND ALTERNATES C BETWEEN HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION UNTIL C CONVERGENCE. C C ON INPUT: C C N = DIMENSION OF X AND THE HOMOTOPY MAP. C C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS. C C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. C C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES. THE ITERATION IS C CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(X,LAMBDA) IS FOUND C SUCH THAT C C |Y(NP1) - 1| <= RELERR + ABSERR AND C C ||Z|| <= RELERR*||X|| + ABSERR , WHERE C C (Z,?) IS THE NEWTON STEP TO Y=(X,LAMBDA). C C Y(1:N+1) = POINT (X(S), LAMBDA(S)) ON ZERO CURVE OF HOMOTOPY MAP. C C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP C AT Y . C C YOLD(1:N+1) = A POINT DIFFERENT FROM Y ON THE ZERO CURVE. C C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY C MAP AT YOLD . C C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. C C QR(1:LENQR) = THE N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X C STORED IN PACKED SKYLINE STORAGE FORMAT. LENQR AND PIVOT C DESCRIBE THE DATA STRUCTURE IN QR . C C LENQR = LENGTH OF THE ONE-DIMENSIONAL ARRAY QR USED TO CONTAIN THE C N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X IN PACKED C SKYLINE STORAGE FORMAT. C C PIVOT(1:N+2) = INDICES OF THE DIAGONAL ELEMENTS OF THE N X N SYMMETRIC C JACOBIAN MATRIX (WITH RESPECT TO X) WITHIN QR . C C WORK(1:13*(N+1)+2*N+LENQR) = WORK ARRAY SPLIT UP AND USED FOR THE C CALCULATION OF THE JACOBIAN MATRIX KERNEL, THE NEWTON STEP, C AND INTERPOLATION. 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 ON OUTPUT: C C N , RELERR , ABSERR , A ARE UNCHANGED. C C NFE HAS BEEN UPDATED. C C IFLAG C = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN. C C = 4 IF THE PRECONDITIONED CONJUGATE GRADIENT ITERATION FAILED TO C CONVERGE (MOST LIKELY DUE TO A JACOBIAN MATRIX WITH RANK < N). C THE ITERATION WAS NOT COMPLETED. C C = 6 IF THE INTERPOLATION/NEWTON ITERATION FAILED TO CONVERGE. C Y AND YOLD CONTAIN THE LAST TWO POINTS FOUND ON THE C ZERO CURVE. C C Y IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT LAMBDA = 1 . C C C CALLS D1MACH , DAXPY , DCOPY , DNRM2 , ROOT , TANGNS . C DOUBLE PRECISION ABSERR,AERR,D1MACH,DD001,DD0011,DD01,DD011, $ DELS,DNRM2,F0,F1,FP0,FP1,QOFS,QSOUT,RELERR,RERR,S,SA,SB, $ SOUT,U INTEGER IFLAG,IPP,IRHO,ITANGW,ITZ,IW,IWP,IZ0,IZ1,JUDY,JW, $ LCODE,LENQR,LIMIT,N,NFE,NP1 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),PAR(1) INTEGER PIVOT(N+2),IPAR(1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C C THE LIMIT ON THE NUMBER OF ITERATIONS ALLOWED MAY BE CHANGED BY C CHANGING THE FOLLOWING PARAMETER STATEMENT: PARAMETER (LIMIT=20) C C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES. C DD01(F0,F1,DELS)=(F1-F0)/DELS DD001(F0,FP0,F1,DELS)=(DD01(F0,F1,DELS)-FP0)/DELS DD011(F0,F1,FP1,DELS)=(FP1-DD01(F0,F1,DELS))/DELS DD0011(F0,FP0,F1,FP1,DELS)=(DD011(F0,F1,FP1,DELS) - $ DD001(F0,FP0,F1,DELS))/DELS QOFS(F0,FP0,F1,FP1,DELS,S)=((DD0011(F0,FP0,F1,FP1,DELS)*(S-DELS) + $ DD001(F0,FP0,F1,DELS))*S + FP0)*S + F0 C C U=D1MACH(4) RERR=MAX(RELERR,U) AERR=MAX(ABSERR,0.0D0) NP1=N+1 IPP=1 IRHO=N+1 IW=IRHO+N IWP=IW+NP1 ITZ=IWP+NP1 IZ0=ITZ+NP1 IZ1=IZ0+NP1 ITANGW=IZ1+NP1 C C ***** MAIN LOOP. ***** C 100 DO 300 JUDY=1,LIMIT DO 110 JW=1,NP1 WORK(ITZ+JW-1)=Y(JW)-YOLD(JW) 110 CONTINUE DELS=DNRM2(NP1,WORK(ITZ),1) C C USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT C THE HERMITE CUBIC INTERPOLANT Q(S). THEN USE ROOT TO FIND THE S C CORRESPONDING TO LAMBDA = 1 . THE TWO POINTS ON THE ZERO CURVE ARE C ALWAYS CHOSEN TO BRACKET LAMBDA=1, WITH THE BRACKETING INTERVAL C ALWAYS BEING [0, DELS]. C SA=0.0 SB=DELS LCODE=1 130 CALL ROOT(SOUT,QSOUT,SA,SB,RERR,AERR,LCODE) IF (LCODE .GT. 0) GO TO 140 QSOUT=QOFS(YOLD(NP1),YPOLD(NP1),Y(NP1),YP(NP1),DELS,SOUT) - 1.0 GO TO 130 C IF LAMBDA = 1 WERE BRACKETED, ROOT CANNOT FAIL. 140 IF (LCODE .GT. 2) THEN IFLAG=6 RETURN ENDIF C C CALCULATE Q(SA) AS THE INITIAL POINT FOR A NEWTON ITERATION. DO 150 JW=1,NP1 WORK(IW+JW-1)=QOFS(YOLD(JW),YPOLD(JW),Y(JW),YP(JW),DELS,SA) 150 CONTINUE C CALCULATE NEWTON STEP AT Q(SA). CALL TANGNS(SA,WORK(IW),WORK(IWP),WORK(ITZ),YPOLD,A,QR,LENQR, $ PIVOT,WORK(IPP),WORK(IRHO),WORK(ITANGW),NFE,N,IFLAG, $ PAR,IPAR) IF (IFLAG .GT. 0) RETURN C NEXT POINT = CURRENT POINT + NEWTON STEP. CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1) C GET THE TANGENT WP AT W AND THE NEXT NEWTON STEP IN TZ . CALL TANGNS(SA,WORK(IW),WORK(IWP),WORK(ITZ),YPOLD,A,QR,LENQR, $ PIVOT,WORK(IPP),WORK(IRHO),WORK(ITANGW),NFE,N,IFLAG, $ PAR,IPAR) IF (IFLAG .GT. 0) RETURN C TAKE NEWTON STEP AND CHECK CONVERGENCE. CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1) IF ((ABS(WORK(IW+N)-1.0) .LE. RERR+AERR) .AND. $ (DNRM2(N,WORK(ITZ),1) .LE. RERR*DNRM2(N,WORK(IW),1)+AERR)) $ THEN CALL DCOPY(NP1,WORK(IW),1,Y,1) RETURN ENDIF C IF THE ITERATION HAS NOT CONVERGED, DISCARD ONE OF THE OLD POINTS C SUCH THAT LAMBDA = 1 IS STILL BRACKETED. IF ((YOLD(NP1)-1.0)*(WORK(IW+N)-1.0) .GT. 0.0) THEN CALL DCOPY(NP1,WORK(IW),1,YOLD,1) CALL DCOPY(NP1,WORK(IWP),1,YPOLD,1) ELSE CALL DCOPY(NP1,WORK(IW),1,Y,1) CALL DCOPY(NP1,WORK(IWP),1,YP,1) ENDIF 300 CONTINUE C C ***** END OF MAIN LOOP. ***** C C THE ALTERNATING OSCULATORY CUBIC INTERPOLATION AND NEWTON ITERATION C HAS NOT CONVERGED IN LIMIT STEPS. ERROR RETURN. IFLAG=6 RETURN END