share/hompack/fortran/rootns.f

   1       SUBROUTINE ROOTNS(N,NFE,IFLAG,RELERR,ABSERR,Y,YP,YOLD,YPOLD,
   2      $   A,QR,LENQR,PIVOT,WORK,PAR,IPAR)
   3 C
   4 C ROOTNS  FINDS THE POINT  YBAR = (XBAR, 1)  ON THE ZERO CURVE OF THE
   5 C HOMOTOPY MAP.  IT STARTS WITH TWO POINTS YOLD=(XOLD,LAMBDAOLD) AND
   6 C Y=(X,LAMBDA) SUCH THAT  LAMBDAOLD < 1 <= LAMBDA , AND ALTERNATES
   7 C BETWEEN HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION UNTIL
   8 C CONVERGENCE.
   9 C
  10 C ON INPUT:
  11 C
  12 C N = DIMENSION OF X AND THE HOMOTOPY MAP.
  13 C
  14 C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
  15 C
  16 C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
  17 C
  18 C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS
  19 C    CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(X,LAMBDA) IS FOUND
  20 C    SUCH THAT
  21 C
  22 C    |Y(NP1) - 1| <= RELERR + ABSERR              AND
  23 C
  24 C    ||Z|| <= RELERR*||X|| + ABSERR  ,          WHERE
  25 C
  26 C    (Z,?) IS THE NEWTON STEP TO Y=(X,LAMBDA).
  27 C
  28 C Y(1:N+1) = POINT (X(S), LAMBDA(S)) ON ZERO CURVE OF HOMOTOPY MAP.
  29 C
  30 C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
  31 C    AT  Y .
  32 C
  33 C YOLD(1:N+1) = A POINT DIFFERENT FROM  Y  ON THE ZERO CURVE.
  34 C
  35 C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY
  36 C    MAP AT  YOLD .
  37 C
  38 C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
  39 C
  40 C QR(1:LENQR) = THE N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X
  41 C    STORED IN PACKED SKYLINE STORAGE FORMAT.  LENQR  AND  PIVOT
  42 C    DESCRIBE THE DATA STRUCTURE IN  QR .
  43 C
  44 C LENQR = LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR  USED TO CONTAIN THE
  45 C    N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X IN PACKED
  46 C    SKYLINE STORAGE FORMAT.
  47 C
  48 C PIVOT(1:N+2) = INDICES OF THE DIAGONAL ELEMENTS OF THE N X N SYMMETRIC
  49 C    JACOBIAN MATRIX (WITH RESPECT TO X) WITHIN  QR .
  50 C
  51 C WORK(1:13*(N+1)+2*N+LENQR) = WORK ARRAY SPLIT UP AND USED FOR THE
  52 C    CALCULATION OF THE JACOBIAN MATRIX KERNEL, THE NEWTON STEP,
  53 C    AND INTERPOLATION.
  54 C
  55 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
  56 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
  57 C    RHO, RHOJS.
  58 C
  59 C ON OUTPUT:
  60 C
  61 C N , RELERR , ABSERR , A  ARE UNCHANGED.
  62 C
  63 C NFE  HAS BEEN UPDATED.
  64 C
  65 C IFLAG
  66 C    = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
  67 C
  68 C    = 4 IF THE PRECONDITIONED CONJUGATE GRADIENT ITERATION FAILED TO
  69 C        CONVERGE (MOST LIKELY DUE TO A JACOBIAN MATRIX WITH RANK < N).
  70 C        THE ITERATION WAS NOT COMPLETED.
  71 C
  72 C    = 6 IF THE INTERPOLATION/NEWTON ITERATION FAILED TO CONVERGE.
  73 C        Y  AND  YOLD  CONTAIN THE LAST TWO POINTS FOUND ON THE
  74 C        ZERO CURVE.
  75 C
  76 C Y  IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT  LAMBDA = 1 .
  77 C
  78 C
  79 C CALLS  D1MACH , DAXPY , DCOPY , DNRM2 , ROOT , TANGNS .
  80 C
  81       DOUBLE PRECISION ABSERR,AERR,D1MACH,DD001,DD0011,DD01,DD011,
  82      $   DELS,DNRM2,F0,F1,FP0,FP1,QOFS,QSOUT,RELERR,RERR,S,SA,SB,
  83      $   SOUT,U
  84       INTEGER IFLAG,IPP,IRHO,ITANGW,ITZ,IW,IWP,IZ0,IZ1,JUDY,JW,
  85      $   LCODE,LENQR,LIMIT,N,NFE,NP1
  86 C
  87 C ***** ARRAY DECLARATIONS. *****
  88 C
  89       DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N),
  90      $   QR(LENQR),WORK(13*(N+1)+2*N+LENQR),PAR(1)
  91       INTEGER PIVOT(N+2),IPAR(1)
  92 C
  93 C ***** END OF DIMENSIONAL INFORMATION. *****
  94 C
  95 C THE LIMIT ON THE NUMBER OF ITERATIONS ALLOWED MAY BE CHANGED BY
  96 C CHANGING THE FOLLOWING PARAMETER STATEMENT:
  97       PARAMETER (LIMIT=20)
  98 C
  99 C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
 100 C
 101       DD01(F0,F1,DELS)=(F1-F0)/DELS
 102       DD001(F0,FP0,F1,DELS)=(DD01(F0,F1,DELS)-FP0)/DELS
 103       DD011(F0,F1,FP1,DELS)=(FP1-DD01(F0,F1,DELS))/DELS
 104       DD0011(F0,FP0,F1,FP1,DELS)=(DD011(F0,F1,FP1,DELS) -
 105      $                            DD001(F0,FP0,F1,DELS))/DELS
 106       QOFS(F0,FP0,F1,FP1,DELS,S)=((DD0011(F0,FP0,F1,FP1,DELS)*(S-DELS) +
 107      $   DD001(F0,FP0,F1,DELS))*S + FP0)*S + F0
 108 C
 109 C
 110       U=D1MACH(4)
 111       RERR=MAX(RELERR,U)
 112       AERR=MAX(ABSERR,0.0D0)
 113       NP1=N+1
 114       IPP=1
 115       IRHO=N+1
 116       IW=IRHO+N
 117       IWP=IW+NP1
 118       ITZ=IWP+NP1
 119       IZ0=ITZ+NP1
 120       IZ1=IZ0+NP1
 121       ITANGW=IZ1+NP1
 122 C
 123 C *****  MAIN LOOP.  *****
 124 C
 125 100   DO 300 JUDY=1,LIMIT
 126       DO 110 JW=1,NP1
 127         WORK(ITZ+JW-1)=Y(JW)-YOLD(JW)
 128 110   CONTINUE
 129       DELS=DNRM2(NP1,WORK(ITZ),1)
 130 C
 131 C USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT
 132 C THE HERMITE CUBIC INTERPOLANT Q(S).  THEN USE  ROOT  TO FIND THE S
 133 C CORRESPONDING TO  LAMBDA = 1 .  THE TWO POINTS ON THE ZERO CURVE ARE
 134 C ALWAYS CHOSEN TO BRACKET LAMBDA=1, WITH THE BRACKETING INTERVAL
 135 C ALWAYS BEING [0, DELS].
 136 C
 137       SA=0.0
 138       SB=DELS
 139       LCODE=1
 140 130   CALL ROOT(SOUT,QSOUT,SA,SB,RERR,AERR,LCODE)
 141       IF (LCODE .GT. 0) GO TO 140
 142       QSOUT=QOFS(YOLD(NP1),YPOLD(NP1),Y(NP1),YP(NP1),DELS,SOUT) - 1.0
 143       GO TO 130
 144 C IF LAMBDA = 1 WERE BRACKETED,  ROOT  CANNOT FAIL.
 145 140   IF (LCODE .GT. 2) THEN
 146         IFLAG=6
 147         RETURN
 148       ENDIF
 149 C
 150 C CALCULATE Q(SA) AS THE INITIAL POINT FOR A NEWTON ITERATION.
 151       DO 150 JW=1,NP1
 152         WORK(IW+JW-1)=QOFS(YOLD(JW),YPOLD(JW),Y(JW),YP(JW),DELS,SA)
 153 150   CONTINUE
 154 C CALCULATE NEWTON STEP AT Q(SA).
 155       CALL TANGNS(SA,WORK(IW),WORK(IWP),WORK(ITZ),YPOLD,A,QR,LENQR,
 156      $   PIVOT,WORK(IPP),WORK(IRHO),WORK(ITANGW),NFE,N,IFLAG,
 157      $   PAR,IPAR)
 158       IF (IFLAG .GT. 0) RETURN
 159 C NEXT POINT = CURRENT POINT + NEWTON STEP.
 160       CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1)
 161 C GET THE TANGENT  WP  AT  W  AND THE NEXT NEWTON STEP IN  TZ .
 162       CALL TANGNS(SA,WORK(IW),WORK(IWP),WORK(ITZ),YPOLD,A,QR,LENQR,
 163      $   PIVOT,WORK(IPP),WORK(IRHO),WORK(ITANGW),NFE,N,IFLAG,
 164      $   PAR,IPAR)
 165       IF (IFLAG .GT. 0) RETURN
 166 C TAKE NEWTON STEP AND CHECK CONVERGENCE.
 167       CALL DAXPY(NP1,1.0D0,WORK(ITZ),1,WORK(IW),1)
 168       IF ((ABS(WORK(IW+N)-1.0) .LE. RERR+AERR) .AND.
 169      $   (DNRM2(N,WORK(ITZ),1) .LE. RERR*DNRM2(N,WORK(IW),1)+AERR))
 170      $   THEN
 171         CALL DCOPY(NP1,WORK(IW),1,Y,1)
 172         RETURN
 173       ENDIF
 174 C IF THE ITERATION HAS NOT CONVERGED, DISCARD ONE OF THE OLD POINTS
 175 C SUCH THAT  LAMBDA = 1  IS STILL BRACKETED.
 176       IF ((YOLD(NP1)-1.0)*(WORK(IW+N)-1.0) .GT. 0.0) THEN
 177         CALL DCOPY(NP1,WORK(IW),1,YOLD,1)
 178         CALL DCOPY(NP1,WORK(IWP),1,YPOLD,1)
 179       ELSE
 180         CALL DCOPY(NP1,WORK(IW),1,Y,1)
 181         CALL DCOPY(NP1,WORK(IWP),1,YP,1)
 182       ENDIF
 183 300   CONTINUE
 184 C
 185 C ***** END OF MAIN LOOP. *****
 186 C
 187 C THE ALTERNATING OSCULATORY CUBIC INTERPOLATION AND NEWTON ITERATION
 188 C HAS NOT CONVERGED IN  LIMIT  STEPS.  ERROR RETURN.
 189       IFLAG=6
 190       RETURN
 191       END