share/minpack/fortran/r1mpyq.f

   1       subroutine r1mpyq(m,n,a,lda,v,w)
   2       integer m,n,lda
   3       double precision a(lda,n),v(n),w(n)
   4 c     **********
   5 c
   6 c     subroutine r1mpyq
   7 c
   8 c     given an m by n matrix a, this subroutine computes a*q where
   9 c     q is the product of 2*(n - 1) transformations
  10 c
  11 c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
  12 c
  13 c     and gv(i), gw(i) are givens rotations in the (i,n) plane which
  14 c     eliminate elements in the i-th and n-th planes, respectively.
  15 c     q itself is not given, rather the information to recover the
  16 c     gv, gw rotations is supplied.
  17 c
  18 c     the subroutine statement is
  19 c
  20 c       subroutine r1mpyq(m,n,a,lda,v,w)
  21 c
  22 c     where
  23 c
  24 c       m is a positive integer input variable set to the number
  25 c         of rows of a.
  26 c
  27 c       n is a positive integer input variable set to the number
  28 c         of columns of a.
  29 c
  30 c       a is an m by n array. on input a must contain the matrix
  31 c         to be postmultiplied by the orthogonal matrix q
  32 c         described above. on output a*q has replaced a.
  33 c
  34 c       lda is a positive integer input variable not less than m
  35 c         which specifies the leading dimension of the array a.
  36 c
  37 c       v is an input array of length n. v(i) must contain the
  38 c         information necessary to recover the givens rotation gv(i)
  39 c         described above.
  40 c
  41 c       w is an input array of length n. w(i) must contain the
  42 c         information necessary to recover the givens rotation gw(i)
  43 c         described above.
  44 c
  45 c     subroutines called
  46 c
  47 c       fortran-supplied ... dabs,dsqrt
  48 c
  49 c     argonne national laboratory. minpack project. march 1980.
  50 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  51 c
  52 c     **********
  53       integer i,j,nmj,nm1
  54       double precision cos,one,sin,temp
  55       data one /1.0d0/
  56 c
  57 c     apply the first set of givens rotations to a.
  58 c
  59       nm1 = n - 1
  60       if (nm1 .lt. 1) go to 50
  61       do 20 nmj = 1, nm1
  62          j = n - nmj
  63          if (dabs(v(j)) .gt. one) cos = one/v(j)
  64          if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
  65          if (dabs(v(j)) .le. one) sin = v(j)
  66          if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
  67          do 10 i = 1, m
  68             temp = cos*a(i,j) - sin*a(i,n)
  69             a(i,n) = sin*a(i,j) + cos*a(i,n)
  70             a(i,j) = temp
  71    10       continue
  72    20    continue
  73 c
  74 c     apply the second set of givens rotations to a.
  75 c
  76       do 40 j = 1, nm1
  77          if (dabs(w(j)) .gt. one) cos = one/w(j)
  78          if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
  79          if (dabs(w(j)) .le. one) sin = w(j)
  80          if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
  81          do 30 i = 1, m
  82             temp = cos*a(i,j) + sin*a(i,n)
  83             a(i,n) = -sin*a(i,j) + cos*a(i,n)
  84             a(i,j) = temp
  85    30       continue
  86    40    continue
  87    50 continue
  88       return
  89 c
  90 c     last card of subroutine r1mpyq.
  91 c
  92       end