share/hompack/fortran/upqrqf.f

   1         SUBROUTINE UPQRQF(N,ETA,S,F0,F1,QT,R,W,T)
   2 C
   3 C SUBROUTINE  UPQRQF  PERFORMS A BROYDEN UPDATE ON THE  Q R
   4 C FACTORIZATION OF A MATRIX  A, (AN APPROXIMATION TO J(X0)),
   5 C RESULTING IN THE FACTORIZATION  Q+ R+ OF
   6 C
   7 C       A+  =  A  +  (Y - A*S) (ST)/(ST * S),
   8 C
   9 C (AN APPROXIMATION TO J(X1))
  10 C WHERE S = X1 - X0, ST = S TRANSPOSE,  Y = F(X1) - F(X0).
  11 C
  12 C THE ENTRY POINT  R1UPQF  PERFORMS THE RANK ONE UPDATE ON THE QR
  13 C FACTORIZATION OF
  14 C
  15 C       A+ =  A + Q*(T*ST).
  16 C
  17 C
  18 C ON INPUT:
  19 C
  20 C N  IS THE DIMENSION OF X AND F(X).
  21 C
  22 C ETA  IS A NOISE PARAMETER.  IF (Y-A*S)(I) .LE. ETA*(|F1(I)|+|F0(I)|)
  23 C    FOR 1 .LE. I .LE. N, THEN NO UPDATE IS PERFORMED.
  24 C
  25 C S(1:N) = X1 - X0   (OR S FOR THE ENTRY POINT R1UPQF).
  26 C
  27 C F0(1:N) = F(X0).
  28 C
  29 C F1(1:N) = F(X1).
  30 C
  31 C QT(1:N,1:N)  CONTAINS THE OLD Q TRANSPOSE, WHERE  A = Q*R .
  32 C
  33 C R(1:N*(N+1)/2)  CONTAINS THE OLD R, STORED BY ROWS.
  34 C
  35 C W(1:N), T(1:N)  ARE WORK ARRAYS ( T  CONTAINS THE VECTOR T FOR THE
  36 C    ENTRY POINT  R1UPQF ).
  37 C
  38 C
  39 C ON OUTPUT:
  40 C
  41 C N  AND  ETA  ARE UNCHANGED.
  42 C
  43 C QT  CONTAINS Q+ TRANSPOSE.
  44 C
  45 C R   CONTAINS R+, STORED BY ROWS.
  46 C
  47 C S, F0, F1, W, AND T  HAVE ALL BEEN CHANGED.
  48 C
  49 C
  50 C CALLS  DAXPY, DDOT, AND DNRM2.
  51 C
  52 C ***** DECLARATIONS *****
  53 C
  54 C     FUNCTION DECLARATIONS
  55 C
  56         DOUBLE PRECISION DDOT, DNRM2
  57 C
  58 C     LOCAL VARIABLES
  59 C
  60         DOUBLE PRECISION C, DEN, ONE, SS, WW, YY
  61         INTEGER I, INDEXR, INDXR2, J, K
  62         LOGICAL SKIPUP
  63 C
  64 C     SCALAR ARGUMENTS
  65 C
  66         DOUBLE PRECISION ETA
  67         INTEGER N
  68 C
  69 C     ARRAY DECLARATIONS
  70 C
  71         DOUBLE PRECISION  S(N), F0(N), F1(N), QT(N,N), R(N*(N+1)/2),
  72      $    W(N), T(N), TT(2)
  73 C
  74 C ***** END OF DECLARATIONS *****
  75 C
  76 C ***** FIRST EXECUTABLE STATEMENT *****
  77 C
  78         ONE = 1.0
  79         SKIPUP = .TRUE.
  80 C
  81 C ***** DEFINE T AND S SUCH THAT *****
  82 C
  83 C           A+ = Q*(R + T*ST).
  84 C
  85 C T = R*S.
  86 C
  87         INDEXR = 1
  88         DO 10 I=1,N
  89           T(I) = DDOT(N-I+1,R(INDEXR),1,S(I),1)
  90           INDEXR = INDEXR + N - I + 1
  91   10    CONTINUE
  92 C
  93 C W = Y - Q*T  = Y - A*S.
  94 C
  95         DO 20 I=1,N
  96           W(I) = F1(I) - F0(I) - DDOT(N,QT(1,I),1,T,1)
  97 C
  98 C         IF W(I) IS NOT SMALL, THEN UPDATE MUST BE PERFORMED,
  99 C         OTHERWISE SET W(I) TO 0.
 100 C
 101           IF (ABS(W(I)) .GT. ETA*(ABS(F1(I)) + ABS(F0(I)))) THEN
 102             SKIPUP = .FALSE.
 103           ELSE
 104             W(I) = 0.0
 105           END IF
 106   20    CONTINUE
 107 C
 108 C  IF NO UPDATE IS NECESSARY, THEN RETURN.
 109 C
 110         IF (SKIPUP) RETURN
 111 C
 112 C T = QT*W = QT*Y - R*S.
 113 C
 114         DO 30 I=1,N
 115           T(I) = DDOT(N,QT(I,1),N,W,1)
 116   30    CONTINUE
 117 C
 118 C S = S/(ST*S).
 119 C
 120         DEN = 1.0/DDOT(N,S,1,S,1)
 121         CALL DSCAL(N,DEN,S,1)
 122 C
 123 C ***** END OF COMPUTATION OF  T & S      *****
 124 C       AT THIS POINT,  A+ = Q*(R + T*ST).
 125 C
 126         CALL r1upqf(n, s, t, qt, r, w)
 127         RETURN
 128         END