share/hompack/fortran/rootqf.f

   1         SUBROUTINE ROOTQF(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,
   2      $     YPOLD,A,QT,R,DZ,Z,W,T,F0,F1,PAR,IPAR)
   3 C
   4 C ROOTQF  FINDS THE POINT  YBAR = (1, XBAR)  ON THE ZERO CURVE OF THE
   5 C HOMOTOPY MAP.  IT STARTS WITH TWO POINTS  YOLD=(LAMBDAOLD,XOLD)  AND
   6 C Y=(LAMBDA,X)  SUCH THAT  LAMBDAOLD < 1 <= LAMBDA, AND ALTERNATES
   7 C BETWEEN USING A SECANT METHOD TO FIND A PREDICTED POINT ON THE
   8 C HYPERPLANE  LAMBDA=1, AND TAKING A QUASI-NEWTON STEP TO RETURN TO THE
   9 C ZERO CURVE OF THE HOMOTOPY MAP.
  10 C
  11 C
  12 C ON INPUT:
  13 C
  14 C N = DIMENSION OF X.
  15 C
  16 C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
  17 C
  18 C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
  19 C
  20 C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS
  21 C    CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(LAMBDA,X) IS FOUND
  22 C    SUCH THAT
  23 C
  24 C    |Y(1) - 1| <= RELERR + ABSERR              AND
  25 C
  26 C    ||DZ|| <= RELERR*||Y|| + ABSERR,           WHERE
  27 C
  28 C    DZ  IS THE QUASI-NEWTON STEP TO Y.
  29 C
  30 C Y(1:N+1) = POINT (LAMBDA(S), X(S)) ON ZERO CURVE OF HOMOTOPY MAP.
  31 C
  32 C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
  33 C    AT  Y.
  34 C
  35 C YOLD(1:N+1) = A POINT DIFFERENT FROM  Y  ON THE ZERO CURVE.
  36 C
  37 C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY
  38 C    MAP AT  YOLD.
  39 C
  40 C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
  41 C
  42 C QT(1:N+1,1:N+1) CONTAINS  Q  TRANSPOSE OF THE QR FACTORIZATION OF
  43 C    THE AUGMENTED JACOBIAN MATRIX EVALUATED AT THE POINT Y.
  44 C
  45 C R((N+1)*(N+2)/2) CONTAINS THE UPPER TRIANGLE OF THE R PART OF
  46 C    OF THE QR FACTORIZATION, STORED BY ROWS.
  47 C
  48 C DZ(1:N+1), Z(1:N+1), W(1:N+1), T(1:N+1), F0(1:N+1), F1(1:N+1)
  49 C    ARE WORK ARRAYS USED FOR THE QUASI-NEWTON STEP AND THE SECANT
  50 C    STEP.
  51 C
  52 C PAR(1:*)  AND  IPAR(1:*)  ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
  53 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
  54 C    RHO, RHOJAC.
  55 C
  56 C
  57 C ON OUTPUT:
  58 C
  59 C N, RELERR, ABSERR, AND A  ARE UNCHANGED.
  60 C
  61 C NFE  HAS BEEN UPDATED.
  62 C
  63 C IFLAG
  64 C    = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
  65 C
  66 C    = 4 IF A SINGULAR JACOBIAN MATRIX OCCURRED.  THE
  67 C        ITERATION WAS NOT COMPLETED.
  68 C
  69 C    = 6 IF THE ITERATION FAILED TO CONVERGE.  Y  AND  YOLD  CONTAIN
  70 C        THE LAST TWO POINTS OBTAINED BY QUASI-NEWTON STEPS, AND  YP
  71 C        CONTAINS A POINT OPPOSITE OF THE HYPERPLANE  LAMBDA=1  FROM
  72 C        Y.
  73 C
  74 C Y  IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT  LAMBDA = 1.
  75 C
  76 C YP  AND  YOLD  CONTAIN POINTS NEAR THE SOLUTION.
  77 C
  78 C CALLS  D1MACH, DAXPY, DCOPY, DDOT, DNRM2, F (OR RHO),
  79 C        QRSLQF, ROOT, UPQRQF.
  80 C
  81 C ***** DECLARATIONS *****
  82 C
  83 C     FUNCTION DECLARATIONS
  84 C
  85         DOUBLE PRECISION D1MACH, DDOT, DNRM2, QOFS
  86 C
  87 C     LOCAL VARIABLES
  88 C
  89         DOUBLE PRECISION AERR, DD001, DD0011, DD01, DD011, DELS, ETA,
  90      $     ONE, P0, P1, PP0, PP1, QSOUT, RERR, S, SA, SB, SOUT,
  91      $     U, ZERO
  92         INTEGER ISTEP, I, LCODE, LIMIT,NP1
  93         LOGICAL BRACK
  94 C
  95 C     SCALAR ARGUMENTS
  96 C
  97         DOUBLE PRECISION RELERR, ABSERR
  98         INTEGER N, NFE, IFLAG
  99 C
 100 C     ARRAY DECLARATIONS
 101 C
 102         DOUBLE PRECISION Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1), A(N),
 103      $     QT(N+1:N+1), R((N+1)*(N+2)/2), DZ(N+1), Z(N+1), W(N+1),
 104      $     T(N+1), F0(N+1), F1(N+1), PAR(1)
 105         INTEGER IPAR(1)
 106 C
 107 C ***** END OF DECLARATIONS *****
 108 C
 109 C
 110 C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
 111 C
 112         DD01(P0,P1,DELS)=(P1-P0)/DELS
 113         DD001(P0,PP0,P1,DELS)=(DD01(P0,P1,DELS)-PP0)/DELS
 114         DD011(P0,P1,PP1,DELS)=(PP1-DD01(P0,P1,DELS))/DELS
 115         DD0011(P0,PP0,P1,PP1,DELS)=(DD011(P0,P1,PP1,DELS) -
 116      $                              DD001(P0,PP0,P1,DELS))/DELS
 117         QOFS(P0,PP0,P1,PP1,DELS,S)=((DD0011(P0,PP0,P1,PP1,DELS)*
 118      $     (S-DELS) + DD001(P0,PP0,P1,DELS))*S + PP0)*S + P0
 119 C
 120 C ***** FIRST EXECUTABLE STATEMENT *****
 121 C
 122 C ***** INITIALIZATION *****
 123 C
 124 C ETA = PARAMETER FOR BROYDEN'S UPDATE.
 125 C LIMIT = MAXIMUM NUMBER OF ITERATIONS ALLOWED.
 126 C
 127         ONE=1.0
 128         ZERO=0.0
 129         U=D1MACH(4)
 130         RERR=MAX(RELERR,U)
 131         AERR=MAX(ABSERR,ZERO)
 132         NP1=N+1
 133         ETA = 100.0*U
 134         LIMIT = 2*(INT(-LOG10(AERR+RERR*DNRM2(NP1,Y,1)))+1)
 135 C
 136 C F0 = (RHO(Y), YP*Y) TRANSPOSE.
 137 C
 138         IF (IFLAG .EQ. -2) THEN
 139 C
 140 C         CURVE TRACKING PROBLEM.
 141 C
 142           CALL RHO(A,Y(1),Y(2),F0,PAR,IPAR)
 143         ELSE IF (IFLAG .EQ. -1) THEN
 144 C
 145 C         ZERO FINDING PROBLEM.
 146 C
 147           CALL F(Y(2),F0)
 148           DO 10 I=1,N
 149             F0(I) = Y(1)*F0(I) + (1.0-Y(1))*(Y(I+1)-A(I))
 150 10        CONTINUE
 151         ELSE
 152 C
 153 C         FIXED POINT PROBLEM.
 154 C
 155           CALL F(Y(2),F0)
 156           DO 20 I=1,N
 157             F0(I) = Y(1)*(A(I)-F0(I))+Y(I+1)-A(I)
 158 20        CONTINUE
 159         END IF
 160         F0(NP1) = DDOT(NP1,YP,1,Y,1)
 161 C
 162 C ***** END OF INITIALIZATION BLOCK *****
 163 C
 164 C ***** COMPUTE FIRST INTERPOLANT WITH A HERMITE CUBIC *****
 165 C
 166 C FIND DISTANCE BETWEEN Y AND YOLD.  DZ=||Y-YOLD||.
 167 C
 168         CALL DCOPY(NP1,Y,1,DZ,1)
 169         CALL DAXPY(NP1,-ONE,YOLD,1,DZ,1)
 170         DELS=DNRM2(NP1,DZ,1)
 171 C
 172 C USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT
 173 C THE HERMITE CUBIC INTERPOLANT Q(S).  THEN USE  ROOT  TO FIND THE S
 174 C CORRESPONDING TO  LAMBDA = 1.  THE TWO POINTS ON THE ZERO CURVE ARE
 175 C ALWAYS CHOSEN TO BRACKET  LAMBDA=1, WITH THE BRACKETING INTERVAL
 176 C ALWAYS BEING [0, DELS].
 177 C
 178         SA=0.0
 179         SB=DELS
 180         LCODE=1
 181 40      CALL ROOT(SOUT,QSOUT,SA,SB,RERR,AERR,LCODE)
 182           IF (LCODE .GT. 0) GO TO 50
 183           QSOUT=QOFS(YOLD(1),YPOLD(1),Y(1),YP(1),DELS,SOUT) - 1.0
 184           GO TO 40
 185 C
 186 C     IF  LAMBDA = 1  WERE BRACKETED,  ROOT  CANNOT FAIL.
 187 C
 188 50      IF (LCODE .GT. 2) THEN
 189           IFLAG=6
 190           RETURN
 191         ENDIF
 192 C
 193 C CALCULATE  Q(SA)  AS THE INITIAL POINT FOR A NEWTON ITERATION.
 194 C
 195         DO 60 I=1,NP1
 196           Z(I)=QOFS(YOLD(I),YPOLD(I),Y(I),YP(I),DELS,SA)
 197 60      CONTINUE
 198 C
 199 C CALCULATE DZ = Z-Y.
 200 C
 201         CALL DCOPY(NP1,Z,1,DZ,1)
 202         CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
 203 C
 204 C ***** END OF CALCULATION OF CUBIC INTERPOLANT *****
 205 C
 206 C TANGENT INFORMATION  YPOLD  IS NO LONGER NEEDED.  HEREAFTER,  YPOLD
 207 C REPRESENTS THE MOST RECENT POINT WHICH IS ON THE OPPOSITE SIDE OF
 208 C LAMBDA=1  FROM  Y.
 209 C
 210 C ***** PREPARE FOR MAIN LOOP *****
 211 C
 212         CALL DCOPY(NP1,YOLD,1,YPOLD,1)
 213 C
 214 C INITIALIZE  BRACK  TO INDICATE THAT THE POINTS  Y  AND  YOLD  BRACKET
 215 C LAMBDA=1,  THUS YOLD = YPOLD.
 216 C
 217         BRACK = .TRUE.
 218 C
 219 C ***** MAIN LOOP *****
 220 C
 221         DO 300 ISTEP=1,LIMIT
 222 C
 223 C UPDATE JACOBIAN MATRIX.
 224 C
 225 C       F1=(RHO(Z), YP*Z) TRANSPOSE.
 226 C
 227           IF (IFLAG .EQ. -2) THEN
 228             CALL RHO(A,Z(1),Z(2),F1,PAR,IPAR)
 229           ELSE IF (IFLAG .EQ. -1) THEN
 230             CALL F(Z(2),F1)
 231             DO 80 I=1,N
 232               F1(I) = Z(1)*F1(I) + (1-Z(1))*(Z(I+1)-A(I))
 233   80        CONTINUE
 234           ELSE
 235             CALL F(Z(2),F1)
 236             DO 90 I=1,N
 237               F1(I) = Z(1)*(A(I)-F1(I))+Z(I+1)-A(I)
 238   90        CONTINUE
 239           END IF
 240           F1(NP1) = DDOT(NP1,YP,1,Z,1)
 241 C
 242 C
 243 C PERFORM BROYDEN UPDATE.
 244 C
 245           CALL UPQRQF(NP1,ETA,DZ,F0,F1,QT,R,W,T)
 246 C
 247 C QUASI-NEWTON STEP.
 248 C
 249 C       COMPUTE NEWTON STEP.
 250 C
 251           CALL DCOPY(N,F1,1,DZ,1)
 252           CALL DSCAL(N,-ONE,DZ,1)
 253           DZ(NP1) = 0.0
 254           CALL QRSLQF(QT,R,DZ,W,NP1)
 255 C
 256 C       TAKE NEWTON STEP.
 257 C
 258           CALL DCOPY(NP1,Z,1,W,1)
 259           CALL DAXPY(NP1,ONE,DZ,1,Z,1)
 260 C
 261 C       CHECK FOR CONVERGENCE.
 262 C
 263           IF ((ABS(Z(1)-1.0) .LE. RERR+AERR) .AND.
 264      $        (DNRM2(NP1,DZ,1) .LE. RERR*DNRM2(N,Z(2),1)+AERR)) THEN
 265              CALL DCOPY(NP1,Z,1,Y,1)
 266              RETURN
 267           END IF
 268 C
 269 C PREPARE FOR NEXT ITERATION.
 270 C
 271 C       F0 = F1.
 272 C
 273           CALL DCOPY(NP1,F1,1,F0,1)
 274 C
 275 C       IF  Z(1) = 1.0  THEN PERFORM QUASI-NEWTON ITERATION AGAIN
 276 C       WITHOUT COMPUTING A NEW PREDICTOR.
 277 C
 278           IF (ABS(Z(1)-1.0) .LE. RERR+AERR) THEN
 279              CALL DCOPY(NP1,Z,1,DZ,1)
 280              CALL DAXPY(NP1,-ONE,W,1,DZ,1)
 281              GOTO 300
 282           END IF
 283 C
 284 C       UPDATE  Y  AND  YOLD.
 285 C
 286           CALL DCOPY(NP1,Y,1,YOLD,1)
 287           CALL DCOPY(NP1,Z,1,Y,1)
 288 C
 289 C       UPDATE  YPOLD  SUCH THAT  YPOLD  IS THE MOST RECENT POINT
 290 C       OPPOSITE OF  LAMBDA=1  FROM  Y.  SET  BRACK = .TRUE.  IFF
 291 C       Y & YOLD  BRACKET  LAMBDA=1  SO THAT  YPOLD=YOLD.
 292 C
 293           IF ((Y(1)-1.0)*(YOLD(1)-1.0) .GT. 0) THEN
 294             BRACK = .FALSE.
 295           ELSE
 296             BRACK = .TRUE.
 297             CALL DCOPY(NP1,YOLD,1,YPOLD,1)
 298           END IF
 299 C
 300 C       COMPUTE DELS = ||Y-YPOLD||.
 301 C
 302           CALL DCOPY(NP1,Y,1,DZ,1)
 303           CALL DAXPY(NP1,-ONE,YPOLD,1,DZ,1)
 304           DELS=DNRM2(NP1,DZ,1)
 305 C
 306 C       COMPUTE  DZ  FOR THE LINEAR PREDICTOR   Z = Y + DZ,
 307 C           WHERE  DZ = SA*(YOLD-Y).
 308 C
 309           SA = (1.0-Y(1))/(YOLD(1)-Y(1))
 310           CALL DCOPY(NP1,YOLD,1,DZ,1)
 311           CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
 312           CALL DSCAL(NP1,SA,DZ,1)
 313 C
 314 C       TO INSURE STABILITY, THE LINEAR PREDICTION MUST BE NO FARTHER
 315 C       FROM  Y  THAN  YPOLD  IS.  THIS IS GUARANTEED IF  BRACK = .TRUE.
 316 C       IF LINEAR PREDICTION IS TOO FAR AWAY, USE BRACKETING POINTS
 317 C       TO COMPUTE LINEAR PREDICTION.
 318 C
 319           IF (.NOT. BRACK) THEN
 320             IF (DNRM2(NP1,DZ,1) .GT. DELS) THEN
 321 C
 322 C             COMPUTE  DZ = SA*(YPOLD-Y).
 323 C
 324               SA = (1.0-Y(1))/(YPOLD(1)-Y(1))
 325               CALL DCOPY(NP1,YPOLD,1,DZ,1)
 326               CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
 327               CALL DSCAL(NP1,SA,DZ,1)
 328             END IF
 329           END IF
 330 C
 331 C       COMPUTE PREDICTOR Z = Y+DZ, AND DZ = NEW Z  - OLD Z (USED FOR
 332 C         QUASI-NEWTON UPDATE).
 333 C
 334           CALL DAXPY(NP1,ONE,DZ,1,Z,1)
 335           CALL DCOPY(NP1,Z,1,DZ,1)
 336           CALL DAXPY(NP1,-ONE,W,1,DZ,1)
 337   300   CONTINUE
 338 C
 339 C ***** END OF MAIN LOOP. *****
 340 C
 341 C THE ALTERNATING OSCULATORY LINEAR PREDICTION AND QUASI-NEWTON
 342 C CORRECTION HAS NOT CONVERGED IN  LIMIT  STEPS.  ERROR RETURN.
 343         IFLAG=6
 344         RETURN
 345 C
 346 C ***** END OF SUBROUTINE ROOTQF *****
 347         END