R1MPYQ Subroutine

  SUBROUTINE R1MPYQ(m,n,a,lda,v,w)

    !  Given an M by N matrix A, this subroutine computes A*Q where
    !  Q is the product of 2*(N - 1) transformations
    !
    !  gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
    !
    !  and GV(I), GW(i) are Givens rotations in the (I,N) plane which
    !  eliminate elements in the I-th and N-th planes, respectively.
    !  Q itself is not given, rather the information to recover the
    !  GV, GW rotations is supplied.
    !
    !  The subroutine statement is
    !
    !  subroutine r1mpyq(m,n,a,lda,v,w)
    !
    !  where
    !
    !  M is a positive integer input variable set to the number
    !  of rows of A.
    !
    !  N is a positive integer input variable set to the number
    !  of columns of A.
    !
    !  A is an M by N array. On input A must contain the matrix
    !  to be postmultiplied by the orthogonal matrix Q
    !  described above. On output A*Q has replaced A.
    !
    !  LDA is a positive integer input variable not less than M
    !  which specifies the leading dimension of the array A.
    !
    !  V is an input array of length N. V(I) must contain the
    !  information necessary to recover the Givens rotation GV(I)
    !  described above.
    !
    !  W is an input array of length N. W(I) must contain the
    !  information necessary to recover the Givens rotation GW(I)
    !  described above.
    !
    !  Argonne National Laboratory. Minpack project. March 1980.
    !  Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More

    IMPLICIT NONE

    INTEGER m,n,lda,i,j,nmj,nm1

    real(dp) a(lda,n),v(n),w(n)
    real(dp) cos1,one,sin1,temp

    one = 1.0D0

    !  Apply the first set of givens rotations to a.

    nm1 = n - 1
    if (nm1 < 1) goto 50
    do nmj = 1, nm1
       j = n - nmj
       if (abs(v(j)) > one) cos1 = one/v(j)
       if (abs(v(j)) > one) sin1 = sqrt(one-cos1**2)
       if (abs(v(j)) <= one) sin1 = v(j)
       if (abs(v(j)) <= one) cos1 = sqrt(one-sin1**2)
       do i = 1, m
          temp  = cos1*a(i,j) - sin1*a(i,n)
          a(i,n) = sin1*a(i,j) + cos1*a(i,n)
          a(i,j) = temp
       end do
    end do

    !  Apply the second set of givens rotations to a.

    do j = 1, nm1
       if (abs(w(j)) > one) cos1 = one/w(j)
       if (abs(w(j)) > one) sin1 = sqrt(one-cos1**2)
       if (abs(w(j)) <= one) sin1 = w(j)
       if (abs(w(j)) <= one) cos1 = sqrt(one-sin1**2)
       do i = 1, m
          temp   =  cos1*a(i,j) + sin1*a(i,n)
          a(i,n) = -sin1*a(i,j) + cos1*a(i,n)
          a(i,j) = temp
       end do
    end do

50  continue

    return
  end SUBROUTINE R1MPYQ