drotmg Subroutine

public subroutine drotmg(dd1, dd2, dx1, dy1, dparam)

construct the modified givens transformation matrix H

Arguments

Type IntentOptional Attributes Name
real(kind=wp) :: dd1
real(kind=wp) :: dd2
real(kind=wp) :: dx1
real(kind=wp) :: dy1
real(kind=wp) :: dparam(5)

Called by

proc~~drotmg~~CalledByGraph proc~drotmg bspline_blas_module::drotmg proc~dwnlit bspline_defc_module::dwnlit proc~dwnlit->proc~drotmg proc~dwnlsm bspline_defc_module::dwnlsm proc~dwnlsm->proc~drotmg proc~dwnlsm->proc~dwnlit proc~dwnnls bspline_defc_module::dwnnls proc~dwnnls->proc~dwnlsm proc~dlpdp bspline_defc_module::dlpdp proc~dlpdp->proc~dwnnls proc~dlsi bspline_defc_module::dlsi proc~dlsi->proc~dlpdp proc~dlsei bspline_defc_module::dlsei proc~dlsei->proc~dlsi proc~dfcmn bspline_defc_module::dfcmn proc~dfcmn->proc~dlsei proc~dfc bspline_defc_module::dfc proc~dfc->proc~dfcmn

Source Code

   subroutine drotmg(dd1,dd2,dx1,dy1,dparam)
      !! construct the modified givens transformation matrix H

      real(wp) :: dd1,dd2,dx1,dy1
      real(wp) :: dparam(5)

      real(wp) :: dflag,dh11,dh12,dh21,dh22,dp1,dp2,dq1,dq2,dtemp,&
                  du

      real(wp),parameter :: zero = 0.0_wp
      real(wp),parameter :: one = 1.0_wp
      real(wp),parameter :: two = 2.0_wp
      real(wp),parameter :: gam    = 4096.0_wp
      real(wp),parameter :: gamsq  = gam*gam    !! 16777216.0_wp
      real(wp),parameter :: rgamsq = one/gamsq  !! 5.9604645e-8_wp

         if (dd1<zero) then
            ! go zero-h-d-and-dx1..
            dflag = -one
            dh11 = zero
            dh12 = zero
            dh21 = zero
            dh22 = zero
            dd1 = zero
            dd2 = zero
            dx1 = zero
         else
            ! case-dd1-nonnegative
            dp2 = dd2*dy1
            if (dp2==zero) then
               dflag = -two
               dparam(1) = dflag
               return
            end if
            ! regular-case..
            dp1 = dd1*dx1
            dq2 = dp2*dy1
            dq1 = dp1*dx1
            if (abs(dq1)>abs(dq2)) then
               dh21 = -dy1/dx1
               dh12 = dp2/dp1
               du = one - dh12*dh21
               if (du>zero) then
                  dflag = zero
                  dd1 = dd1/du
                  dd2 = dd2/du
                  dx1 = dx1*du
               else
                  ! this code path if here for safety. we do not expect this
                  ! condition to ever hold except in edge cases with rounding
                  ! errors. see doi: 10.1145/355841.355847
                  dflag = -one
                  dh11 = zero
                  dh12 = zero
                  dh21 = zero
                  dh22 = zero
                  dd1 = zero
                  dd2 = zero
                  dx1 = zero
               end if
            else

               if (dq2<zero) then
                  ! go zero-h-d-and-dx1..
                  dflag = -one
                  dh11 = zero
                  dh12 = zero
                  dh21 = zero
                  dh22 = zero
                  dd1 = zero
                  dd2 = zero
                  dx1 = zero
               else
                  dflag = one
                  dh11 = dp1/dp2
                  dh22 = dx1/dy1
                  du = one + dh11*dh22
                  dtemp = dd2/du
                  dd2 = dd1/du
                  dd1 = dtemp
                  dx1 = dy1*du
               end if
            end if

            ! procedure..scale-check
            if (dd1/=zero) then
               do while ((dd1<=rgamsq) .or. (dd1>=gamsq))
                  if (dflag==zero) then
                     dh11 = one
                     dh22 = one
                     dflag = -one
                  else
                     dh21 = -one
                     dh12 = one
                     dflag = -one
                  end if
                  if (dd1<=rgamsq) then
                     dd1 = dd1*gam**2
                     dx1 = dx1/gam
                     dh11 = dh11/gam
                     dh12 = dh12/gam
                  else
                     dd1 = dd1/gam**2
                     dx1 = dx1*gam
                     dh11 = dh11*gam
                     dh12 = dh12*gam
                  end if
               enddo
            end if

            if (dd2/=zero) then
               do while ( (abs(dd2)<=rgamsq) .or. (abs(dd2)>=gamsq) )
                  if (dflag==zero) then
                     dh11 = one
                     dh22 = one
                     dflag = -one
                  else
                     dh21 = -one
                     dh12 = one
                     dflag = -one
                  end if
                  if (abs(dd2)<=rgamsq) then
                     dd2 = dd2*gam**2
                     dh21 = dh21/gam
                     dh22 = dh22/gam
                  else
                     dd2 = dd2/gam**2
                     dh21 = dh21*gam
                     dh22 = dh22*gam
                  end if
               end do
            end if

         end if

         if (dflag<zero) then
            dparam(2) = dh11
            dparam(3) = dh21
            dparam(4) = dh12
            dparam(5) = dh22
         else if (dflag==zero) then
            dparam(3) = dh21
            dparam(4) = dh12
         else
            dparam(2) = dh11
            dparam(5) = dh22
         end if

         dparam(1) = dflag

      end subroutine drotmg