share/minpack/fortran/r1updt.f

   1       subroutine r1updt(m,n,s,ls,u,v,w,sing)
   2       integer m,n,ls
   3       logical sing
   4       double precision s(ls),u(m),v(n),w(m)
   5 c     **********
   6 c
   7 c     subroutine r1updt
   8 c
   9 c     given an m by n lower trapezoidal matrix s, an m-vector u,
  10 c     and an n-vector v, the problem is to determine an
  11 c     orthogonal matrix q such that
  12 c
  13 c                   t
  14 c           (s + u*v )*q
  15 c
  16 c     is again lower trapezoidal.
  17 c
  18 c     this subroutine determines q as the product of 2*(n - 1)
  19 c     transformations
  20 c
  21 c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
  22 c
  23 c     where gv(i), gw(i) are givens rotations in the (i,n) plane
  24 c     which eliminate elements in the i-th and n-th planes,
  25 c     respectively. q itself is not accumulated, rather the
  26 c     information to recover the gv, gw rotations is returned.
  27 c
  28 c     the subroutine statement is
  29 c
  30 c       subroutine r1updt(m,n,s,ls,u,v,w,sing)
  31 c
  32 c     where
  33 c
  34 c       m is a positive integer input variable set to the number
  35 c         of rows of s.
  36 c
  37 c       n is a positive integer input variable set to the number
  38 c         of columns of s. n must not exceed m.
  39 c
  40 c       s is an array of length ls. on input s must contain the lower
  41 c         trapezoidal matrix s stored by columns. on output s contains
  42 c         the lower trapezoidal matrix produced as described above.
  43 c
  44 c       ls is a positive integer input variable not less than
  45 c         (n*(2*m-n+1))/2.
  46 c
  47 c       u is an input array of length m which must contain the
  48 c         vector u.
  49 c
  50 c       v is an array of length n. on input v must contain the vector
  51 c         v. on output v(i) contains the information necessary to
  52 c         recover the givens rotation gv(i) described above.
  53 c
  54 c       w is an output array of length m. w(i) contains information
  55 c         necessary to recover the givens rotation gw(i) described
  56 c         above.
  57 c
  58 c       sing is a logical output variable. sing is set true if any
  59 c         of the diagonal elements of the output s are zero. otherwise
  60 c         sing is set false.
  61 c
  62 c     subprograms called
  63 c
  64 c       minpack-supplied ... dpmpar
  65 c
  66 c       fortran-supplied ... dabs,dsqrt
  67 c
  68 c     argonne national laboratory. minpack project. march 1980.
  69 c     burton s. garbow, kenneth e. hillstrom, jorge j. more,
  70 c     john l. nazareth
  71 c
  72 c     **********
  73       integer i,j,jj,l,nmj,nm1
  74       double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp,
  75      *                 zero
  76       double precision dpmpar
  77       data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
  78 c
  79 c     giant is the largest magnitude.
  80 c
  81       giant = dpmpar(3)
  82 c
  83 c     initialize the diagonal element pointer.
  84 c
  85       jj = (n*(2*m - n + 1))/2 - (m - n)
  86 c
  87 c     move the nontrivial part of the last column of s into w.
  88 c
  89       l = jj
  90       do 10 i = n, m
  91          w(i) = s(l)
  92          l = l + 1
  93    10    continue
  94 c
  95 c     rotate the vector v into a multiple of the n-th unit vector
  96 c     in such a way that a spike is introduced into w.
  97 c
  98       nm1 = n - 1
  99       if (nm1 .lt. 1) go to 70
 100       do 60 nmj = 1, nm1
 101          j = n - nmj
 102          jj = jj - (m - j + 1)
 103          w(j) = zero
 104          if (v(j) .eq. zero) go to 50
 105 c
 106 c        determine a givens rotation which eliminates the
 107 c        j-th element of v.
 108 c
 109          if (dabs(v(n)) .ge. dabs(v(j))) go to 20
 110             cotan = v(n)/v(j)
 111             sin = p5/dsqrt(p25+p25*cotan**2)
 112             cos = sin*cotan
 113             tau = one
 114             if (dabs(cos)*giant .gt. one) tau = one/cos
 115             go to 30
 116    20    continue
 117             tan = v(j)/v(n)
 118             cos = p5/dsqrt(p25+p25*tan**2)
 119             sin = cos*tan
 120             tau = sin
 121    30    continue
 122 c
 123 c        apply the transformation to v and store the information
 124 c        necessary to recover the givens rotation.
 125 c
 126          v(n) = sin*v(j) + cos*v(n)
 127          v(j) = tau
 128 c
 129 c        apply the transformation to s and extend the spike in w.
 130 c
 131          l = jj
 132          do 40 i = j, m
 133             temp = cos*s(l) - sin*w(i)
 134             w(i) = sin*s(l) + cos*w(i)
 135             s(l) = temp
 136             l = l + 1
 137    40       continue
 138    50    continue
 139    60    continue
 140    70 continue
 141 c
 142 c     add the spike from the rank 1 update to w.
 143 c
 144       do 80 i = 1, m
 145          w(i) = w(i) + v(n)*u(i)
 146    80    continue
 147 c
 148 c     eliminate the spike.
 149 c
 150       sing = .false.
 151       if (nm1 .lt. 1) go to 140
 152       do 130 j = 1, nm1
 153          if (w(j) .eq. zero) go to 120
 154 c
 155 c        determine a givens rotation which eliminates the
 156 c        j-th element of the spike.
 157 c
 158          if (dabs(s(jj)) .ge. dabs(w(j))) go to 90
 159             cotan = s(jj)/w(j)
 160             sin = p5/dsqrt(p25+p25*cotan**2)
 161             cos = sin*cotan
 162             tau = one
 163             if (dabs(cos)*giant .gt. one) tau = one/cos
 164             go to 100
 165    90    continue
 166             tan = w(j)/s(jj)
 167             cos = p5/dsqrt(p25+p25*tan**2)
 168             sin = cos*tan
 169             tau = sin
 170   100    continue
 171 c
 172 c        apply the transformation to s and reduce the spike in w.
 173 c
 174          l = jj
 175          do 110 i = j, m
 176             temp = cos*s(l) + sin*w(i)
 177             w(i) = -sin*s(l) + cos*w(i)
 178             s(l) = temp
 179             l = l + 1
 180   110       continue
 181 c
 182 c        store the information necessary to recover the
 183 c        givens rotation.
 184 c
 185          w(j) = tau
 186   120    continue
 187 c
 188 c        test for zero diagonal elements in the output s.
 189 c
 190          if (s(jj) .eq. zero) sing = .true.
 191          jj = jj + (m - j + 1)
 192   130    continue
 193   140 continue
 194 c
 195 c     move w back into the last column of the output s.
 196 c
 197       l = jj
 198       do 150 i = n, m
 199          s(l) = w(i)
 200          l = l + 1
 201   150    continue
 202       if (s(jj) .eq. zero) sing = .true.
 203       return
 204 c
 205 c     last card of subroutine r1updt.
 206 c
 207       end