qrupdt Subroutine

  subroutine qrupdt(r,qt,u,v)
    real(8), dimension(:,:), intent(inout) :: r,qt
    real(8), dimension(:), intent(inout) :: u
    real(8), dimension(:), intent(in) :: v
    integer :: i,k,n
    n=assert_eqn((/size(r,1),size(r,2),size(qt,1),size(qt,2),size(u),&
         size(v)/),'qrupdt')
    k=n+1-ifirstloc(u(n:1:-1) /= 0.0)
    if (k < 1) k=1
    do i=k-1,1,-1
       call rotate(r,qt,i,u(i),-u(i+1))
       u(i)=pythag(u(i),u(i+1))
    end do
    r(1,:)=r(1,:)+u(1)*v
    do i=1,k-1
       call rotate(r,qt,i,r(i,i),-r(i+1,i))
    end do
  contains
    subroutine rotate(r,qt,i,a,b)
      real(8), dimension(:,:), target, intent(inout) :: r,qt
      integer, intent(in) :: i
      real(8), intent(in) :: a,b
      real(8), dimension(size(r,1)) :: temp
      integer :: n
      real(8) :: c,fact,s
      n=assert_eq4(size(r,1),size(r,2),size(qt,1),size(qt,2),'rotate')
      if (a == 0.0) then
         c=0.0
         s=sign(1.0d0,b)
      else if (abs(a) > abs(b)) then
         fact=b/a
         c=sign(1.0d0/sqrt(1.0d0+fact**2),a)
         s=fact*c
      else
         fact=a/b
         s=sign(1.0d0/sqrt(1.0d0+fact**2),b)
         c=fact*s
      end if
      temp(i:n)=r(i,i:n)
      r(i,i:n)=c*temp(i:n)-s*r(i+1,i:n)
      r(i+1,i:n)=s*temp(i:n)+c*r(i+1,i:n)
      temp=qt(i,:)
      qt(i,:)=c*temp-s*qt(i+1,:)
      qt(i+1,:)=s*temp+c*qt(i+1,:)
    end subroutine rotate
    !    
    function pythag(a,b)
      real(8), intent(in) :: a,b
      real(8) :: pythag
      real(8) :: absa,absb
      absa=abs(a)
      absb=abs(b)
      if (absa > absb) then
         pythag=absa*sqrt(1.0d0+(absb/absa)**2)
      else
         if (absb == 0.0) then
            pythag=0.0
         else
            pythag=absb*sqrt(1.0d0+(absa/absb)**2)
         end if
      end if
    end function pythag
  end subroutine qrupdt