share/minpack/fortran/rwupdt.f

   1       subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
   2       integer n,ldr
   3       double precision alpha
   4       double precision r(ldr,n),w(n),b(n),cos(n),sin(n)
   5 c     **********
   6 c
   7 c     subroutine rwupdt
   8 c
   9 c     given an n by n upper triangular matrix r, this subroutine
  10 c     computes the qr decomposition of the matrix formed when a row
  11 c     is added to r. if the row is specified by the vector w, then
  12 c     rwupdt determines an orthogonal matrix q such that when the
  13 c     n+1 by n matrix composed of r augmented by w is premultiplied
  14 c     by (q transpose), the resulting matrix is upper trapezoidal.
  15 c     the matrix (q transpose) is the product of n transformations
  16 c
  17 c           g(n)*g(n-1)* ... *g(1)
  18 c
  19 c     where g(i) is a givens rotation in the (i,n+1) plane which
  20 c     eliminates elements in the (n+1)-st plane. rwupdt also
  21 c     computes the product (q transpose)*c where c is the
  22 c     (n+1)-vector (b,alpha). q itself is not accumulated, rather
  23 c     the information to recover the g rotations is supplied.
  24 c
  25 c     the subroutine statement is
  26 c
  27 c       subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
  28 c
  29 c     where
  30 c
  31 c       n is a positive integer input variable set to the order of r.
  32 c
  33 c       r is an n by n array. on input the upper triangular part of
  34 c         r must contain the matrix to be updated. on output r
  35 c         contains the updated triangular matrix.
  36 c
  37 c       ldr is a positive integer input variable not less than n
  38 c         which specifies the leading dimension of the array r.
  39 c
  40 c       w is an input array of length n which must contain the row
  41 c         vector to be added to r.
  42 c
  43 c       b is an array of length n. on input b must contain the
  44 c         first n elements of the vector c. on output b contains
  45 c         the first n elements of the vector (q transpose)*c.
  46 c
  47 c       alpha is a variable. on input alpha must contain the
  48 c         (n+1)-st element of the vector c. on output alpha contains
  49 c         the (n+1)-st element of the vector (q transpose)*c.
  50 c
  51 c       cos is an output array of length n which contains the
  52 c         cosines of the transforming givens rotations.
  53 c
  54 c       sin is an output array of length n which contains the
  55 c         sines of the transforming givens rotations.
  56 c
  57 c     subprograms called
  58 c
  59 c       fortran-supplied ... dabs,dsqrt
  60 c
  61 c     argonne national laboratory. minpack project. march 1980.
  62 c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
  63 c     jorge j. more
  64 c
  65 c     **********
  66       integer i,j,jm1
  67       double precision cotan,one,p5,p25,rowj,tan,temp,zero
  68       data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
  69 c
  70       do 60 j = 1, n
  71          rowj = w(j)
  72          jm1 = j - 1
  73 c
  74 c        apply the previous transformations to
  75 c        r(i,j), i=1,2,...,j-1, and to w(j).
  76 c
  77          if (jm1 .lt. 1) go to 20
  78          do 10 i = 1, jm1
  79             temp = cos(i)*r(i,j) + sin(i)*rowj
  80             rowj = -sin(i)*r(i,j) + cos(i)*rowj
  81             r(i,j) = temp
  82    10       continue
  83    20    continue
  84 c
  85 c        determine a givens rotation which eliminates w(j).
  86 c
  87          cos(j) = one
  88          sin(j) = zero
  89          if (rowj .eq. zero) go to 50
  90          if (dabs(r(j,j)) .ge. dabs(rowj)) go to 30
  91             cotan = r(j,j)/rowj
  92             sin(j) = p5/dsqrt(p25+p25*cotan**2)
  93             cos(j) = sin(j)*cotan
  94             go to 40
  95    30    continue
  96             tan = rowj/r(j,j)
  97             cos(j) = p5/dsqrt(p25+p25*tan**2)
  98             sin(j) = cos(j)*tan
  99    40    continue
 100 c
 101 c        apply the current transformation to r(j,j), b(j), and alpha.
 102 c
 103          r(j,j) = cos(j)*r(j,j) + sin(j)*rowj
 104          temp = cos(j)*b(j) + sin(j)*alpha
 105          alpha = -sin(j)*b(j) + cos(j)*alpha
 106          b(j) = temp
 107    50    continue
 108    60    continue
 109       return
 110 c
 111 c     last card of subroutine rwupdt.
 112 c
 113       end