SRC/zlartg.f

   1 *> \brief \b ZLARTG generates a plane rotation with real cosine and complex sine.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLARTG + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlartg.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlartg.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlartg.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLARTG( F, G, CS, SN, R )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       DOUBLE PRECISION   CS
  25 *       COMPLEX*16         F, G, R, SN
  26 *       ..
  27 *
  28 *
  29 *> \par Purpose:
  30 *  =============
  31 *>
  32 *> \verbatim
  33 *>
  34 *> ZLARTG generates a plane rotation so that
  35 *>
  36 *>    [  CS  SN  ]     [ F ]     [ R ]
  37 *>    [  __      ]  .  [   ]  =  [   ]   where CS**2 + |SN|**2 = 1.
  38 *>    [ -SN  CS  ]     [ G ]     [ 0 ]
  39 *>
  40 *> This is a faster version of the BLAS1 routine ZROTG, except for
  41 *> the following differences:
  42 *>    F and G are unchanged on return.
  43 *>    If G=0, then CS=1 and SN=0.
  44 *>    If F=0, then CS=0 and SN is chosen so that R is real.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] F
  51 *> \verbatim
  52 *>          F is COMPLEX*16
  53 *>          The first component of vector to be rotated.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] G
  57 *> \verbatim
  58 *>          G is COMPLEX*16
  59 *>          The second component of vector to be rotated.
  60 *> \endverbatim
  61 *>
  62 *> \param[out] CS
  63 *> \verbatim
  64 *>          CS is DOUBLE PRECISION
  65 *>          The cosine of the rotation.
  66 *> \endverbatim
  67 *>
  68 *> \param[out] SN
  69 *> \verbatim
  70 *>          SN is COMPLEX*16
  71 *>          The sine of the rotation.
  72 *> \endverbatim
  73 *>
  74 *> \param[out] R
  75 *> \verbatim
  76 *>          R is COMPLEX*16
  77 *>          The nonzero component of the rotated vector.
  78 *> \endverbatim
  79 *
  80 *  Authors:
  81 *  ========
  82 *
  83 *> \author Univ. of Tennessee
  84 *> \author Univ. of California Berkeley
  85 *> \author Univ. of Colorado Denver
  86 *> \author NAG Ltd.
  87 *
  88 *> \date November 2013
  89 *
  90 *> \ingroup complex16OTHERauxiliary
  91 *
  92 *> \par Further Details:
  93 *  =====================
  94 *>
  95 *> \verbatim
  96 *>
  97 *>  3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
  98 *>
  99 *>  This version has a few statements commented out for thread safety
 100 *>  (machine parameters are computed on each entry). 10 feb 03, SJH.
 101 *> \endverbatim
 102 *>
 103 *  =====================================================================
 104       SUBROUTINE ZLARTG( F, G, CS, SN, R )
 105 *
 106 *  -- LAPACK auxiliary routine (version 3.5.0) --
 107 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 108 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 109 *     November 2013
 110 *
 111 *     .. Scalar Arguments ..
 112       DOUBLE PRECISION   CS
 113       COMPLEX*16         F, G, R, SN
 114 *     ..
 115 *
 116 *  =====================================================================
 117 *
 118 *     .. Parameters ..
 119       DOUBLE PRECISION   TWO, ONE, ZERO
 120       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
 121       COMPLEX*16         CZERO
 122       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
 123 *     ..
 124 *     .. Local Scalars ..
 125 *     LOGICAL            FIRST
 126       INTEGER            COUNT, I
 127       DOUBLE PRECISION   D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
 128      $                   SAFMN2, SAFMX2, SCALE
 129       COMPLEX*16         FF, FS, GS
 130 *     ..
 131 *     .. External Functions ..
 132       DOUBLE PRECISION   DLAMCH, DLAPY2
 133       LOGICAL            DISNAN
 134       EXTERNAL           DLAMCH, DLAPY2, DISNAN
 135 *     ..
 136 *     .. Intrinsic Functions ..
 137       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
 138      $                   MAX, SQRT
 139 *     ..
 140 *     .. Statement Functions ..
 141       DOUBLE PRECISION   ABS1, ABSSQ
 142 *     ..
 143 *     .. Statement Function definitions ..
 144       ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
 145       ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
 146 *     ..
 147 *     .. Executable Statements ..
 148 *
 149       SAFMIN = DLAMCH( 'S' )
 150       EPS = DLAMCH( 'E' )
 151       SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
 152      $         LOG( DLAMCH( 'B' ) ) / TWO )
 153       SAFMX2 = ONE / SAFMN2
 154       SCALE = MAX( ABS1( F ), ABS1( G ) )
 155       FS = F
 156       GS = G
 157       COUNT = 0
 158       IF( SCALE.GE.SAFMX2 ) THEN
 159    10    CONTINUE
 160          COUNT = COUNT + 1
 161          FS = FS*SAFMN2
 162          GS = GS*SAFMN2
 163          SCALE = SCALE*SAFMN2
 164          IF( SCALE.GE.SAFMX2 )
 165      $      GO TO 10
 166       ELSE IF( SCALE.LE.SAFMN2 ) THEN
 167          IF( G.EQ.CZERO.OR.DISNAN( ABS( G ) ) ) THEN
 168             CS = ONE
 169             SN = CZERO
 170             R = F
 171             RETURN
 172          END IF
 173    20    CONTINUE
 174          COUNT = COUNT - 1
 175          FS = FS*SAFMX2
 176          GS = GS*SAFMX2
 177          SCALE = SCALE*SAFMX2
 178          IF( SCALE.LE.SAFMN2 )
 179      $      GO TO 20
 180       END IF
 181       F2 = ABSSQ( FS )
 182       G2 = ABSSQ( GS )
 183       IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
 184 *
 185 *        This is a rare case: F is very small.
 186 *
 187          IF( F.EQ.CZERO ) THEN
 188             CS = ZERO
 189             R = DLAPY2( DBLE( G ), DIMAG( G ) )
 190 *           Do complex/real division explicitly with two real divisions
 191             D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
 192             SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
 193             RETURN
 194          END IF
 195          F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
 196 *        G2 and G2S are accurate
 197 *        G2 is at least SAFMIN, and G2S is at least SAFMN2
 198          G2S = SQRT( G2 )
 199 *        Error in CS from underflow in F2S is at most
 200 *        UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
 201 *        If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
 202 *        and so CS .lt. sqrt(SAFMIN)
 203 *        If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
 204 *        and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
 205 *        Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
 206          CS = F2S / G2S
 207 *        Make sure abs(FF) = 1
 208 *        Do complex/real division explicitly with 2 real divisions
 209          IF( ABS1( F ).GT.ONE ) THEN
 210             D = DLAPY2( DBLE( F ), DIMAG( F ) )
 211             FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
 212          ELSE
 213             DR = SAFMX2*DBLE( F )
 214             DI = SAFMX2*DIMAG( F )
 215             D = DLAPY2( DR, DI )
 216             FF = DCMPLX( DR / D, DI / D )
 217          END IF
 218          SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
 219          R = CS*F + SN*G
 220       ELSE
 221 *
 222 *        This is the most common case.
 223 *        Neither F2 nor F2/G2 are less than SAFMIN
 224 *        F2S cannot overflow, and it is accurate
 225 *
 226          F2S = SQRT( ONE+G2 / F2 )
 227 *        Do the F2S(real)*FS(complex) multiply with two real multiplies
 228          R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
 229          CS = ONE / F2S
 230          D = F2 + G2
 231 *        Do complex/real division explicitly with two real divisions
 232          SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
 233          SN = SN*DCONJG( GS )
 234          IF( COUNT.NE.0 ) THEN
 235             IF( COUNT.GT.0 ) THEN
 236                DO 30 I = 1, COUNT
 237                   R = R*SAFMX2
 238    30          CONTINUE
 239             ELSE
 240                DO 40 I = 1, -COUNT
 241                   R = R*SAFMN2
 242    40          CONTINUE
 243             END IF
 244          END IF
 245       END IF
 246       RETURN
 247 *
 248 *     End of ZLARTG
 249 *
 250       END