dh12 Subroutine

private subroutine dh12(Mode, Lpivot, l1, m, u, Iue, Up, c, Ice, Icv, Ncv)

Construction and/or application of a single Householder transformation. Q = I + U*(U**T)/B

Reference

  • C.L.Lawson and R.J.Hanson, Jet Propulsion Laboratory, 1973 Jun 12 to appear in 'Solving Least Squares Problems', Prentice-Hall, 1974

Revision history

  • 790101 DATE WRITTEN
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 890831 Modified array declarations. (WRB)
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900328 Added TYPE section. (WRB)
  • 900911 Added DDOT to real(wp) statement. (WRB)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Mode

1 or 2 to select algorithm H1 or H2 .
integer, intent(in) :: Lpivot

the index of the pivot element.
integer, intent(in) :: l1

If L1 <= M the transformation will be constructed to zero elements indexed from L1 through M. If L1 > M the subroutine does an identity transformation.
integer, intent(in) :: m

see l1
real(kind=wp), intent(inout) :: u(Iue,*)

On entry to H1 U() contains the pivot vector. On exit from H1 U() and UP contain quantities defining the vector U of the Householder transformation. On entry to H2 U() and UP should contain quantities previously computed by H1. These will not be modified by H2.
integer, intent(in) :: Iue

the storage increment between elements of U.
real(kind=wp), intent(inout) :: Up

see u
real(kind=wp), intent(inout) :: c(*)

On entry to H1 or H2 C() contains a matrix which will be regarded as a set of vectors to which the Householder transformation is to be applied. On exit C() contains the set of transformed vectors.
integer, intent(in) :: Ice

Storage increment between elements of vectors in C().
integer, intent(in) :: Icv

Storage increment between vectors in C().
integer, intent(in) :: Ncv

Number of vectors in C() to be transformed. If NCV <= 0 no operations will be done on C().

Calls

proc~~dh12~~CallsGraph proc~dh12 bspline_defc_module::dh12 proc~daxpy bspline_blas_module::daxpy proc~dh12->proc~daxpy proc~ddot bspline_blas_module::ddot proc~dh12->proc~ddot proc~dswap bspline_blas_module::dswap proc~dh12->proc~dswap

Called by

proc~~dh12~~CalledByGraph proc~dh12 bspline_defc_module::dh12 proc~dbndac bspline_defc_module::dbndac proc~dbndac->proc~dh12 proc~dhfti bspline_defc_module::dhfti proc~dhfti->proc~dh12 proc~dlsei bspline_defc_module::dlsei proc~dlsei->proc~dh12 proc~dlsi bspline_defc_module::dlsi proc~dlsei->proc~dlsi proc~dlsi->proc~dh12 proc~dlsi->proc~dhfti proc~dlpdp bspline_defc_module::dlpdp proc~dlsi->proc~dlpdp proc~dwnlit bspline_defc_module::dwnlit proc~dwnlit->proc~dh12 proc~dwnlsm bspline_defc_module::dwnlsm proc~dwnlsm->proc~dh12 proc~dwnlsm->proc~dwnlit proc~defcmn bspline_defc_module::defcmn proc~defcmn->proc~dbndac proc~dfcmn bspline_defc_module::dfcmn proc~dfcmn->proc~dbndac proc~dfcmn->proc~dlsei proc~dwnnls bspline_defc_module::dwnnls proc~dwnnls->proc~dwnlsm proc~defc bspline_defc_module::defc proc~defc->proc~defcmn proc~dfc bspline_defc_module::dfc proc~dfc->proc~dfcmn proc~dlpdp->proc~dwnnls

Source Code

    subroutine dh12(Mode, Lpivot, l1, m, u, Iue, Up, c, Ice, Icv, Ncv)

    integer,intent(in) :: Mode !! 1 or 2   to select algorithm  H1  or  H2 .
    integer,intent(in) :: Lpivot !! the index of the pivot element.
    integer,intent(in) :: l1 !! If `L1 <= M` the transformation will be constructed to
                             !! zero elements indexed from `L1` through `M`. If `L1 > M`
                             !! the subroutine does an identity transformation.
    integer,intent(in) :: m !! see `l1`
    integer,intent(in) :: Iue !! the storage increment between elements of `U`.
    real(wp),intent(inout) :: u(Iue, *) !! On entry to H1 `U()` contains the pivot vector.
                                        !! On exit from H1 `U()` and `UP`
                                        !! contain quantities defining the vector `U` of the
                                        !! Householder transformation.   On entry to H2 `U()`
                                        !! and `UP` should contain quantities previously computed
                                        !! by H1.  These will not be modified by H2.
    real(wp),intent(inout)  :: Up !! see `u`
    real(wp),intent(inout) :: c(*) !! On entry to H1 or H2 `C()` contains a matrix which will be
                                   !! regarded as a set of vectors to which the Householder
                                   !! transformation is to be applied.  On exit `C()` contains the
                                   !! set of transformed vectors.
    integer,intent(in) :: Ice !! Storage increment between elements of vectors in `C()`.
    integer,intent(in) :: Icv !! Storage increment between vectors in `C()`.
    integer,intent(in) :: Ncv !! Number of vectors in `C()` to be transformed. If `NCV <= 0`
                              !! no operations will be done on `C()`.

      integer :: i, i2, i3, i4, incr, j, kl1, &
                 kl2, klp, l1m1, mml1p2
      real(wp) :: b, cl, clinv, ul1m1, sm

      real(wp), parameter :: one = 1.0_wp

      if (0 < Lpivot .and. Lpivot < l1 .and. l1 <= m) then
         cl = abs(u(1, Lpivot))
         if (Mode /= 2) then
            ! ****** CONSTRUCT THE TRANSFORMATION. ******
            do j = l1, m
               cl = max(abs(u(1, j)), cl)
            end do
            if (cl <= 0.0_wp) return
            clinv = one/cl
            sm = (u(1, Lpivot)*clinv)**2
            do j = l1, m
               sm = sm + (u(1, j)*clinv)**2
            end do
            cl = cl*sqrt(sm)
            if (u(1, Lpivot) > 0.0_wp) cl = -cl
            Up = u(1, Lpivot) - cl
            u(1, Lpivot) = cl
            ! ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ******
         elseif (cl <= 0.0_wp) then
            return
         end if

         if (Ncv > 0) then
            b = Up*u(1, Lpivot)
            ! B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN.
            if (b < 0.0_wp) then
               b = one/b
               mml1p2 = m - l1 + 2
               if (mml1p2 <= 20) then
                  i2 = 1 - Icv + Ice*(Lpivot - 1)
                  incr = Ice*(l1 - Lpivot)
                  do j = 1, Ncv
                     i2 = i2 + Icv
                     i3 = i2 + incr
                     i4 = i3
                     sm = c(i2)*Up
                     do i = l1, m
                        sm = sm + c(i3)*u(1, i)
                        i3 = i3 + Ice
                     end do
                     if (sm /= 0.0_wp) then
                        sm = sm*b
                        c(i2) = c(i2) + sm*Up
                        do i = l1, m
                           c(i4) = c(i4) + sm*u(1, i)
                           i4 = i4 + Ice
                        end do
                     end if
                  end do
               else
                  l1m1 = l1 - 1
                  kl1 = 1 + (l1m1 - 1)*Ice
                  kl2 = kl1
                  klp = 1 + (Lpivot - 1)*Ice
                  ul1m1 = u(1, l1m1)
                  u(1, l1m1) = Up
                  if (Lpivot /= l1m1) call dswap(Ncv, c(kl1), Icv, c(klp), Icv)
                  do j = 1, Ncv
                     sm = ddot(mml1p2, u(1, l1m1), Iue, c(kl1), Ice)
                     sm = sm*b
                     call daxpy(mml1p2, sm, u(1, l1m1), Iue, c(kl1), Ice)
                     kl1 = kl1 + Icv
                  end do
                  u(1, l1m1) = ul1m1
                  if (Lpivot /= l1m1) then
                     kl1 = kl2
                     call dswap(Ncv, c(kl1), Icv, c(klp), Icv)
                  end if
               end if
            end if
         end if
      end if

   end subroutine dh12