website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/rwupdt.f

   1 *DECK RWUPDT
   2       SUBROUTINE RWUPDT (N, R, LDR, W, B, ALPHA, COS, SIN)
   3 C***BEGIN PROLOGUE  RWUPDT
   4 C***SUBSIDIARY
   5 C***PURPOSE  Subsidiary to SNLS1 and SNLS1E
   6 C***LIBRARY   SLATEC
   7 C***TYPE      SINGLE PRECISION (RWUPDT-S, DWUPDT-D)
   8 C***AUTHOR  (UNKNOWN)
   9 C***DESCRIPTION
  10 C
  11 C     Given an N by N upper triangular matrix R, this subroutine
  12 C     computes the QR decomposition of the matrix formed when a row
  13 C     is added to R. If the row is specified by the vector W, then
  14 C     RWUPDT determines an orthogonal matrix Q such that when the
  15 C     N+1 by N matrix composed of R augmented by W is premultiplied
  16 C     by (Q TRANSPOSE), the resulting matrix is upper trapezoidal.
  17 C     The orthogonal matrix Q is the product of N transformations
  18 C
  19 C           G(1)*G(2)* ... *G(N)
  20 C
  21 C     where G(I) is a Givens rotation in the (I,N+1) plane which
  22 C     eliminates elements in the I-th plane. RWUPDT also
  23 C     computes the product (Q TRANSPOSE)*C where C is the
  24 C     (N+1)-vector (b,alpha). Q itself is not accumulated, rather
  25 C     the information to recover the G rotations is supplied.
  26 C
  27 C     The subroutine statement is
  28 C
  29 C       SUBROUTINE RWUPDT(N,R,LDR,W,B,ALPHA,COS,SIN)
  30 C
  31 C     where
  32 C
  33 C       N is a positive integer input variable set to the order of R.
  34 C
  35 C       R is an N by N array. On input the upper triangular part of
  36 C         R must contain the matrix to be updated. On output R
  37 C         contains the updated triangular matrix.
  38 C
  39 C       LDR is a positive integer input variable not less than N
  40 C         which specifies the leading dimension of the array R.
  41 C
  42 C       W is an input array of length N which must contain the row
  43 C         vector to be added to R.
  44 C
  45 C       B is an array of length N. On input B must contain the
  46 C         first N elements of the vector C. On output B contains
  47 C         the first N elements of the vector (Q TRANSPOSE)*C.
  48 C
  49 C       ALPHA is a variable. On input ALPHA must contain the
  50 C         (N+1)-st element of the vector C. On output ALPHA contains
  51 C         the (N+1)-st element of the vector (Q TRANSPOSE)*C.
  52 C
  53 C       COS is an output array of length N which contains the
  54 C         cosines of the transforming Givens rotations.
  55 C
  56 C       SIN is an output array of length N which contains the
  57 C         sines of the transforming Givens rotations.
  58 C
  59 C***SEE ALSO  SNLS1, SNLS1E
  60 C***ROUTINES CALLED  (NONE)
  61 C***REVISION HISTORY  (YYMMDD)
  62 C   800301  DATE WRITTEN
  63 C   890831  Modified array declarations.  (WRB)
  64 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  65 C   900326  Removed duplicate information from DESCRIPTION section.
  66 C           (WRB)
  67 C   900328  Added TYPE section.  (WRB)
  68 C***END PROLOGUE  RWUPDT
  69       INTEGER N,LDR
  70       REAL ALPHA
  71       REAL R(LDR,*),W(*),B(*),COS(*),SIN(*)
  72       INTEGER I,J,JM1
  73       REAL COTAN,ONE,P5,P25,ROWJ,TAN,TEMP,ZERO
  74       SAVE ONE, P5, P25, ZERO
  75       DATA ONE,P5,P25,ZERO /1.0E0,5.0E-1,2.5E-1,0.0E0/
  76 C***FIRST EXECUTABLE STATEMENT  RWUPDT
  77       DO 60 J = 1, N
  78          ROWJ = W(J)
  79          JM1 = J - 1
  80 C
  81 C        APPLY THE PREVIOUS TRANSFORMATIONS TO
  82 C        R(I,J), I=1,2,...,J-1, AND TO W(J).
  83 C
  84          IF (JM1 .LT. 1) GO TO 20
  85          DO 10 I = 1, JM1
  86             TEMP = COS(I)*R(I,J) + SIN(I)*ROWJ
  87             ROWJ = -SIN(I)*R(I,J) + COS(I)*ROWJ
  88             R(I,J) = TEMP
  89    10       CONTINUE
  90    20    CONTINUE
  91 C
  92 C        DETERMINE A GIVENS ROTATION WHICH ELIMINATES W(J).
  93 C
  94          COS(J) = ONE
  95          SIN(J) = ZERO
  96          IF (ROWJ .EQ. ZERO) GO TO 50
  97          IF (ABS(R(J,J)) .GE. ABS(ROWJ)) GO TO 30
  98             COTAN = R(J,J)/ROWJ
  99             SIN(J) = P5/SQRT(P25+P25*COTAN**2)
 100             COS(J) = SIN(J)*COTAN
 101             GO TO 40
 102    30    CONTINUE
 103             TAN = ROWJ/R(J,J)
 104             COS(J) = P5/SQRT(P25+P25*TAN**2)
 105             SIN(J) = COS(J)*TAN
 106    40    CONTINUE
 107 C
 108 C        APPLY THE CURRENT TRANSFORMATION TO R(J,J), B(J), AND ALPHA.
 109 C
 110          R(J,J) = COS(J)*R(J,J) + SIN(J)*ROWJ
 111          TEMP = COS(J)*B(J) + SIN(J)*ALPHA
 112          ALPHA = -SIN(J)*B(J) + COS(J)*ALPHA
 113          B(J) = TEMP
 114    50    CONTINUE
 115    60    CONTINUE
 116       RETURN
 117 C
 118 C     LAST CARD OF SUBROUTINE RWUPDT.
 119 C
 120       END