srotmg.f

Go to the documentation of this file.

*> \brief \b SROTMG
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SROTMG(SD1,SD2,SX1,SY1,SPARAM)
*
*       .. Scalar Arguments ..
*       REAL SD1,SD2,SX1,SY1
*       ..
*       .. Array Arguments ..
*       REAL SPARAM(5)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
*>    THE SECOND COMPONENT OF THE 2-VECTOR  (SQRT(SD1)*SX1,SQRT(SD2)*>    SY2)**T.
*>    WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
*>
*>    SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
*>
*>      (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
*>    H=(          )    (          )    (          )    (          )
*>      (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
*>    LOCATIONS 2-4 OF SPARAM CONTAIN SH11,SH21,SH12, AND SH22
*>    RESPECTIVELY. (VALUES OF 1.E0, -1.E0, OR 0.E0 IMPLIED BY THE
*>    VALUE OF SPARAM(1) ARE NOT STORED IN SPARAM.)
*>
*>    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 SD1 AND SD2.  ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in,out] SD1
*> \verbatim
*>          SD1 is REAL
*> \endverbatim
*>
*> \param[in,out] SD2
*> \verbatim
*>          SD2 is REAL
*> \endverbatim
*>
*> \param[in,out] SX1
*> \verbatim
*>          SX1 is REAL
*> \endverbatim
*>
*> \param[in] SY1
*> \verbatim
*>          SY1 is REAL
*> \endverbatim
*>
*> \param[out] SPARAM
*> \verbatim
*>          SPARAM is REAL array, dimension (5)
*>     SPARAM(1)=SFLAG
*>     SPARAM(2)=SH11
*>     SPARAM(3)=SH21
*>     SPARAM(4)=SH12
*>     SPARAM(5)=SH22
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup rotmg
*
*  =====================================================================
      SUBROUTINE srotmg(SD1,SD2,SX1,SY1,SPARAM)
*
*  -- 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 ..
      REAL SD1,SD2,SX1,SY1
*     ..
*     .. Array Arguments ..
      REAL SPARAM(5)
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      REAL GAM,GAMSQ,ONE,RGAMSQ,SFLAG,SH11,SH12,SH21,SH22,SP1,SP2,SQ1,
     $     SQ2,STEMP,SU,TWO,ZERO
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC abs
*     ..
*     .. Data statements ..
*
      DATA zero,one,two/0.e0,1.e0,2.e0/
      DATA gam,gamsq,rgamsq/4096.e0,1.67772e7,5.96046e-8/
*     ..
 
      IF (sd1.LT.zero) THEN
*        GO ZERO-H-D-AND-SX1..
         sflag = -one
         sh11 = zero
         sh12 = zero
         sh21 = zero
         sh22 = zero
*
         sd1 = zero
         sd2 = zero
         sx1 = zero
      ELSE
*        CASE-SD1-NONNEGATIVE
         sp2 = sd2*sy1
         IF (sp2.EQ.zero) THEN
            sflag = -two
            sparam(1) = sflag
            RETURN
         END IF
*        REGULAR-CASE..
         sp1 = sd1*sx1
         sq2 = sp2*sy1
         sq1 = sp1*sx1
*
         IF (abs(sq1).GT.abs(sq2)) THEN
            sh21 = -sy1/sx1
            sh12 = sp2/sp1
*
            su = one - sh12*sh21
*
           IF (su.GT.zero) THEN
             sflag = zero
             sd1 = sd1/su
             sd2 = sd2/su
             sx1 = sx1*su
           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
             sflag = -one
             sh11 = zero
             sh12 = zero
             sh21 = zero
             sh22 = zero
*
             sd1 = zero
             sd2 = zero
             sx1 = zero
           END IF
         ELSE
 
            IF (sq2.LT.zero) THEN
*              GO ZERO-H-D-AND-SX1..
               sflag = -one
               sh11 = zero
               sh12 = zero
               sh21 = zero
               sh22 = zero
*
               sd1 = zero
               sd2 = zero
               sx1 = zero
            ELSE
               sflag = one
               sh11 = sp1/sp2
               sh22 = sx1/sy1
               su = one + sh11*sh22
               stemp = sd2/su
               sd2 = sd1/su
               sd1 = stemp
               sx1 = sy1*su
            END IF
         END IF
 
*     PROCEDURE..SCALE-CHECK
         IF (sd1.NE.zero) THEN
            DO WHILE ((sd1.LE.rgamsq) .OR. (sd1.GE.gamsq))
               IF (sflag.EQ.zero) THEN
                  sh11 = one
                  sh22 = one
                  sflag = -one
               ELSE
                  sh21 = -one
                  sh12 = one
                  sflag = -one
               END IF
               IF (sd1.LE.rgamsq) THEN
                  sd1 = sd1*gam**2
                  sx1 = sx1/gam
                  sh11 = sh11/gam
                  sh12 = sh12/gam
               ELSE
                  sd1 = sd1/gam**2
                  sx1 = sx1*gam
                  sh11 = sh11*gam
                  sh12 = sh12*gam
               END IF
            ENDDO
         END IF
 
         IF (sd2.NE.zero) THEN
            DO WHILE ( (abs(sd2).LE.rgamsq) .OR. (abs(sd2).GE.gamsq) )
               IF (sflag.EQ.zero) THEN
                  sh11 = one
                  sh22 = one
                  sflag = -one
               ELSE
                  sh21 = -one
                  sh12 = one
                  sflag = -one
               END IF
               IF (abs(sd2).LE.rgamsq) THEN
                  sd2 = sd2*gam**2
                  sh21 = sh21/gam
                  sh22 = sh22/gam
               ELSE
                  sd2 = sd2/gam**2
                  sh21 = sh21*gam
                  sh22 = sh22*gam
               END IF
            END DO
         END IF
 
      END IF
 
      IF (sflag.LT.zero) THEN
         sparam(2) = sh11
         sparam(3) = sh21
         sparam(4) = sh12
         sparam(5) = sh22
      ELSE IF (sflag.EQ.zero) THEN
         sparam(3) = sh21
         sparam(4) = sh12
      ELSE
         sparam(2) = sh11
         sparam(5) = sh22
      END IF
 
      sparam(1) = sflag
      RETURN
*
*     End of SROTMG
*
      END