share/hompack/fortran/r1upqf.f

   1         SUBROUTINE R1UPQF(N,S,T,QT,R,W)
   2 C ***** DECLARATIONS *****
   3 C
   4 C     FUNCTION DECLARATIONS
   5 C
   6         DOUBLE PRECISION DDOT, DNRM2
   7 C
   8 C     LOCAL VARIABLES
   9 C
  10         DOUBLE PRECISION C, DEN, ONE, SS, WW, YY, temp
  11         INTEGER I, INDEXR, INDXR2, J, K
  12         LOGICAL SKIPUP
  13 C
  14 C     SCALAR ARGUMENTS
  15 C
  16         DOUBLE PRECISION ETA
  17         INTEGER N
  18 C
  19 C     ARRAY DECLARATIONS
  20 C
  21         DOUBLE PRECISION  S(N), QT(N,N), R(N*(N+1)/2),
  22      $    W(N), T(N), TT(2)
  23 C
  24 C ***** END OF DECLARATIONS *****
  25 C
  26 C ***** COMPUTE THE QR FACTORIZATION Q- R- OF (R + T*S).  THEN,  *****
  27 C       Q+ = Q*Q-,  AND  R+ = R-.
  28 C
  29 C FIND THE LARGEST  K  SUCH THAT  T(K) .NE. 0.
  30 C
  31         K = N
  32   50    IF (T(K) .NE. 0.0 .OR. K .LE. 1) GOTO 60
  33           K=K-1
  34           GOTO 50
  35   60    CONTINUE
  36 C
  37 C COMPUTE THE INDEX OF R(K-1,K-1).
  38 C
  39         INDEXR = (N + N - K + 3)*(K - 2) / 2 + 1
  40 C
  41 C ***** TRANSFORM R+T*ST INTO AN UPPER HESSENBERG MATRIX *****
  42 C
  43 C DETERMINE JACOBI ROTATIONS WHICH WILL ZERO OUT ROWS
  44 C N, N-1,...,2  OF THE MATRIX  T*ST,  AND APPLY THESE
  45 C ROTATIONS TO  R.  (THIS IS EQUIVALENT TO APPLYING THE
  46 C SAME ROTATIONS TO  R+T*ST, EXCEPT FOR THE FIRST ROW.
  47 C THUS, AFTER AN ADJUSTMENT FOR THE FIRST ROW, THE
  48 C RESULT IS AN UPPER HESSENBERG MATRIX.  THE
  49 C SUBDIAGONAL ELEMENTS OF WHICH WILL BE STORED IN  W.
  50 C
  51 C NOTE:  ROWS N,N-1,...,K+1 ARE ALREADY ALL ZERO.
  52 C
  53         DO 90 I=K-1,1,-1
  54 C
  55 C         DETERMINE THE JACOBI ROTATION WHICH WILL ZERO OUT
  56 C         ROW  I+1  OF THE  T*ST  MATRIX.
  57 C
  58           IF (T(I) .EQ. 0.0) THEN
  59             C = 0.0
  60 C         SS = SIGN(-T(I+1))= -T(I+1)/|T(I+1)|
  61             SS = -SIGN(ONE,T(I+1))
  62           ELSE
  63             DEN = DNRM2(2,T(I),1)
  64             C = T(I) / DEN
  65             SS = -T(I+1)/DEN
  66           END IF
  67 C
  68 C         PREMULTIPLY  R  BY THE JACOBI ROTATION.
  69 C
  70           YY = R(INDEXR)
  71           WW = 0.0
  72           R(INDEXR) = C*YY - SS*WW
  73           W(I+1) = SS*YY + C*WW
  74           INDEXR = INDEXR + 1
  75           INDXR2 = INDEXR + N - I
  76           DO 70 J= I+1,N
  77 C           YY = R(I,J)
  78 C           WW = R(I+1,J)
  79               YY = R(INDEXR)
  80               WW = R(INDXR2)
  81 C           R(I,J) = C*YY - SS*WW
  82 C           R(I+1,J) = SS*YY + C*WW
  83               R(INDEXR) = C*YY - SS*WW
  84               R(INDXR2) = SS*YY + C*WW
  85             INDEXR = INDEXR + 1
  86             INDXR2 = INDXR2 + 1
  87   70      CONTINUE
  88 C
  89 C         PREMULTIPLY  QT  BY THE JACOBI ROTATION.
  90 C
  91           DO 80 J=1,N
  92             YY = QT(I,J)
  93             WW = QT(I+1,J)
  94             QT(I,J) = C*YY - SS*WW
  95             QT(I+1,J) = SS*YY + C*WW
  96   80      CONTINUE
  97 C
  98 C         UPDATE  T(I)  SO THAT  T(I)*ST(J)  IS THE  (I,J)TH  COMPONENT
  99 C         OF  T*ST, PREMULTIPLIED BY ALL OF THE JACOBI ROTATIONS SO
 100 C         FAR.
 101 C
 102           IF (T(I) .EQ. 0.0) THEN
 103             T(I) = ABS(T(I+1))
 104           ELSE
 105             T(I) = DNRM2(2,T(I),1)
 106           END IF
 107 C
 108 C         LET INDEXR = THE INDEX OF R(I-1,I-1).
 109 C
 110           INDEXR = INDEXR - 2*(N - I) - 3
 111 C
 112   90    CONTINUE
 113 C
 114 C     UPDATE THE FIRST ROW OF  R  SO THAT  R  HOLDS  (R+T*ST)
 115 C     PREMULTIPLIED BY ALL OF THE ABOVE JACOBI ROTATIONS.
 116 C
 117         temp = T(1)
 118         CALL DAXPY(N,temp,S,1,R,1)
 119 C
 120 C ***** END OF TRANSFORMATION TO UPPER HESSENBERG *****
 121 C
 122 C
 123 C ***** TRANSFORM UPPER HESSENBERG MATRIX INTO UPPER *****
 124 C       TRIANGULAR MATRIX.
 125 C
 126 C       INDEXR = INDEX OF R(1,1).
 127 C
 128           INDEXR = 1
 129           DO 120 I=1,K-1
 130 C
 131 C           DETERMINE APPROPRIATE JACOBI ROTATION TO ZERO OUT
 132 C           R(I+1,I).
 133 C
 134             IF (R(INDEXR) .EQ. 0.0) THEN
 135               C = 0.0
 136               SS = -SIGN(ONE,W(I+1))
 137             ELSE
 138               TT(1) = R(INDEXR)
 139               TT(2) = W(I+1)
 140               DEN = DNRM2(2,TT,1)
 141               C = R(INDEXR) / DEN
 142               SS = -W(I+1)/DEN
 143             END IF
 144 C
 145 C           PREMULTIPLY  R  BY JACOBI ROTATION.
 146 C
 147             YY = R(INDEXR)
 148             WW = W(I+1)
 149             R(INDEXR) = C*YY - SS*WW
 150             W(I+1) = 0.0
 151             INDEXR = INDEXR + 1
 152             INDXR2 = INDEXR + N - I
 153             DO 100 J= I+1,N
 154 C             YY = R(I,J)
 155 C             WW = R(I+1,J)
 156                 YY = R(INDEXR)
 157                 WW = R(INDXR2)
 158 C             R(I,J) = C*YY -SS*WW
 159 C             R(I+1,J) = SS*YY + C*WW
 160                 R(INDEXR) = C*YY - SS*WW
 161                 R(INDXR2) = SS*YY + C*WW
 162               INDEXR = INDEXR + 1
 163               INDXR2 = INDXR2 + 1
 164   100       CONTINUE
 165 C
 166 C           PREMULTIPLY  QT  BY JACOBI ROTATION.
 167 C
 168             DO 110 J=1,N
 169               YY = QT(I,J)
 170               WW = QT(I+1,J)
 171               QT(I,J) = C*YY - SS*WW
 172               QT(I+1,J) = SS*YY + C*WW
 173   110       CONTINUE
 174   120     CONTINUE
 175 C
 176 C ***** END OF TRANSFORMATION TO UPPER TRIANGULAR *****
 177 C
 178 C
 179 C ***** END OF UPDATE *****
 180 C
 181 C
 182         RETURN
 183 C
 184 C ***** END OF SUBROUTINE UPQRQF *****
 185         END