File:  [local] / rpl / lapack / lapack / zlartg.f
Revision 1.18: download - view: text, annotated - select for diffs - revision graph
Tue May 29 07:18:29 2018 UTC (5 years, 11 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, rpl-4_1_33, rpl-4_1_32, rpl-4_1_31, rpl-4_1_30, rpl-4_1_29, rpl-4_1_28, HEAD
Mise à jour de Lapack.

    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 December 2016
   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.7.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: *     December 2016
  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

CVSweb interface <joel.bertrand@systella.fr>