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[in,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. 
*
*> \date November 2011
*
*> \ingroup double_blas_level1
*
*  =====================================================================
      SUBROUTINE drotmg(DD1,DD2,DX1,DY1,DPARAM)
*
*  -- Reference BLAS level1 routine (version 3.4.0) --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. 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
           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