drotmg.f

Go to the documentation of this file.

*> \brief \b DROTMG
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
*
*       .. Scalar Arguments ..
*       DOUBLE PRECISION DD1,DD2,DX1,DY1
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION DPARAM(5)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
*>    THE SECOND COMPONENT OF THE 2-VECTOR  (DSQRT(DD1)*DX1,DSQRT(DD2)*>    DY2)**T.
*>    WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
*>
*>    DFLAG=-1.D0     DFLAG=0.D0        DFLAG=1.D0     DFLAG=-2.D0
*>
*>      (DH11  DH12)    (1.D0  DH12)    (DH11  1.D0)    (1.D0  0.D0)
*>    H=(          )    (          )    (          )    (          )
*>      (DH21  DH22),   (DH21  1.D0),   (-1.D0 DH22),   (0.D0  1.D0).
*>    LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
*>    RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
*>    VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
*>
*>    THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
*>    INEXACT.  THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
*>    OF DD1 AND DD2.  ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in,out] DD1
*> \verbatim
*>          DD1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DD2
*> \verbatim
*>          DD2 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DX1
*> \verbatim
*>          DX1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] DY1
*> \verbatim
*>          DY1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] DPARAM
*> \verbatim
*>          DPARAM is DOUBLE PRECISION array, dimension (5)
*>     DPARAM(1)=DFLAG
*>     DPARAM(2)=DH11
*>     DPARAM(3)=DH21
*>     DPARAM(4)=DH12
*>     DPARAM(5)=DH22
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup rotmg
*
*  =====================================================================
      SUBROUTINE drotmg(DD1,DD2,DX1,DY1,DPARAM)
*
*  -- Reference BLAS level1 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION DD1,DD2,DX1,DY1
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION DPARAM(5)
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP,
     $                 DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC dabs
*     ..
*     .. Data statements ..
*
      DATA zero,one,two/0.d0,1.d0,2.d0/
      DATA gam,gamsq,rgamsq/4096.d0,16777216.d0,5.9604645d-8/
*     ..
 
      IF (dd1.LT.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.EQ.zero) THEN
            dflag = -two
            dparam(1) = dflag
            RETURN
         END IF
*        REGULAR-CASE..
         dp1 = dd1*dx1
         dq2 = dp2*dy1
         dq1 = dp1*dx1
*
         IF (dabs(dq1).GT.dabs(dq2)) THEN
            dh21 = -dy1/dx1
            dh12 = dp2/dp1
*
            du = one - dh12*dh21
*
           IF (du.GT.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.LT.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.NE.zero) THEN
            DO WHILE ((dd1.LE.rgamsq) .OR. (dd1.GE.gamsq))
               IF (dflag.EQ.zero) THEN
                  dh11 = one
                  dh22 = one
                  dflag = -one
               ELSE
                  dh21 = -one
                  dh12 = one
                  dflag = -one
               END IF
               IF (dd1.LE.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.NE.zero) THEN
            DO WHILE ( (dabs(dd2).LE.rgamsq) .OR. (dabs(dd2).GE.gamsq) )
               IF (dflag.EQ.zero) THEN
                  dh11 = one
                  dh22 = one
                  dflag = -one
               ELSE
                  dh21 = -one
                  dh12 = one
                  dflag = -one
               END IF
               IF (dabs(dd2).LE.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.LT.zero) THEN
         dparam(2) = dh11
         dparam(3) = dh21
         dparam(4) = dh12
         dparam(5) = dh22
      ELSE IF (dflag.EQ.zero) THEN
         dparam(3) = dh21
         dparam(4) = dh12
      ELSE
         dparam(2) = dh11
         dparam(5) = dh22
      END IF
 
      dparam(1) = dflag
      RETURN
*
*     End of DROTMG
*
      END