construct the modified givens transformation matrix H
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=wp) | :: | dd1 | ||||
real(kind=wp) | :: | dd2 | ||||
real(kind=wp) | :: | dx1 | ||||
real(kind=wp) | :: | dy1 | ||||
real(kind=wp) | :: | dparam(5) |
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