Return to zlartg.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

    1:       SUBROUTINE ZLARTG( F, G, CS, SN, R )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       DOUBLE PRECISION   CS
   10:       COMPLEX*16         F, G, R, SN
   11: *     ..
   12: *
   13: *  Purpose
   14: *  =======
   15: *
   16: *  ZLARTG generates a plane rotation so that
   17: *
   18: *     [  CS  SN  ]     [ F ]     [ R ]
   19: *     [  __      ]  .  [   ]  =  [   ]   where CS**2 + |SN|**2 = 1.
   20: *     [ -SN  CS  ]     [ G ]     [ 0 ]
   21: *
   22: *  This is a faster version of the BLAS1 routine ZROTG, except for
   23: *  the following differences:
   24: *     F and G are unchanged on return.
   25: *     If G=0, then CS=1 and SN=0.
   26: *     If F=0, then CS=0 and SN is chosen so that R is real.
   27: *
   28: *  Arguments
   29: *  =========
   30: *
   31: *  F       (input) COMPLEX*16
   32: *          The first component of vector to be rotated.
   33: *
   34: *  G       (input) COMPLEX*16
   35: *          The second component of vector to be rotated.
   36: *
   37: *  CS      (output) DOUBLE PRECISION
   38: *          The cosine of the rotation.
   39: *
   40: *  SN      (output) COMPLEX*16
   41: *          The sine of the rotation.
   42: *
   43: *  R       (output) COMPLEX*16
   44: *          The nonzero component of the rotated vector.
   45: *
   46: *  Further Details
   47: *  ======= =======
   48: *
   49: *  3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
   50: *
   51: *  This version has a few statements commented out for thread safety
   52: *  (machine parameters are computed on each entry). 10 feb 03, SJH.
   53: *
   54: *  =====================================================================
   55: *
   56: *     .. Parameters ..
   57:       DOUBLE PRECISION   TWO, ONE, ZERO
   58:       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
   59:       COMPLEX*16         CZERO
   60:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
   61: *     ..
   62: *     .. Local Scalars ..
   63: *     LOGICAL            FIRST
   64:       INTEGER            COUNT, I
   65:       DOUBLE PRECISION   D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
   66:      $                   SAFMN2, SAFMX2, SCALE
   67:       COMPLEX*16         FF, FS, GS
   68: *     ..
   69: *     .. External Functions ..
   70:       DOUBLE PRECISION   DLAMCH, DLAPY2
   71:       EXTERNAL           DLAMCH, DLAPY2
   72: *     ..
   73: *     .. Intrinsic Functions ..
   74:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
   75:      $                   MAX, SQRT
   76: *     ..
   77: *     .. Statement Functions ..
   78:       DOUBLE PRECISION   ABS1, ABSSQ
   79: *     ..
   80: *     .. Save statement ..
   81: *     SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2
   82: *     ..
   83: *     .. Data statements ..
   84: *     DATA               FIRST / .TRUE. /
   85: *     ..
   86: *     .. Statement Function definitions ..
   87:       ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
   88:       ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
   89: *     ..
   90: *     .. Executable Statements ..
   91: *
   92: *     IF( FIRST ) THEN
   93:          SAFMIN = DLAMCH( 'S' )
   94:          EPS = DLAMCH( 'E' )
   95:          SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
   96:      $            LOG( DLAMCH( 'B' ) ) / TWO )
   97:          SAFMX2 = ONE / SAFMN2
   98: *        FIRST = .FALSE.
   99: *     END IF
  100:       SCALE = MAX( ABS1( F ), ABS1( G ) )
  101:       FS = F
  102:       GS = G
  103:       COUNT = 0
  104:       IF( SCALE.GE.SAFMX2 ) THEN
  105:    10    CONTINUE
  106:          COUNT = COUNT + 1
  107:          FS = FS*SAFMN2
  108:          GS = GS*SAFMN2
  109:          SCALE = SCALE*SAFMN2
  110:          IF( SCALE.GE.SAFMX2 )
  111:      $      GO TO 10
  112:       ELSE IF( SCALE.LE.SAFMN2 ) THEN
  113:          IF( G.EQ.CZERO ) THEN
  114:             CS = ONE
  115:             SN = CZERO
  116:             R = F
  117:             RETURN
  118:          END IF
  119:    20    CONTINUE
  120:          COUNT = COUNT - 1
  121:          FS = FS*SAFMX2
  122:          GS = GS*SAFMX2
  123:          SCALE = SCALE*SAFMX2
  124:          IF( SCALE.LE.SAFMN2 )
  125:      $      GO TO 20
  126:       END IF
  127:       F2 = ABSSQ( FS )
  128:       G2 = ABSSQ( GS )
  129:       IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
  130: *
  131: *        This is a rare case: F is very small.
  132: *
  133:          IF( F.EQ.CZERO ) THEN
  134:             CS = ZERO
  135:             R = DLAPY2( DBLE( G ), DIMAG( G ) )
  136: *           Do complex/real division explicitly with two real divisions
  137:             D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
  138:             SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
  139:             RETURN
  140:          END IF
  141:          F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
  142: *        G2 and G2S are accurate
  143: *        G2 is at least SAFMIN, and G2S is at least SAFMN2
  144:          G2S = SQRT( G2 )
  145: *        Error in CS from underflow in F2S is at most
  146: *        UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
  147: *        If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
  148: *        and so CS .lt. sqrt(SAFMIN)
  149: *        If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
  150: *        and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
  151: *        Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
  152:          CS = F2S / G2S
  153: *        Make sure abs(FF) = 1
  154: *        Do complex/real division explicitly with 2 real divisions
  155:          IF( ABS1( F ).GT.ONE ) THEN
  156:             D = DLAPY2( DBLE( F ), DIMAG( F ) )
  157:             FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
  158:          ELSE
  159:             DR = SAFMX2*DBLE( F )
  160:             DI = SAFMX2*DIMAG( F )
  161:             D = DLAPY2( DR, DI )
  162:             FF = DCMPLX( DR / D, DI / D )
  163:          END IF
  164:          SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
  165:          R = CS*F + SN*G
  166:       ELSE
  167: *
  168: *        This is the most common case.
  169: *        Neither F2 nor F2/G2 are less than SAFMIN
  170: *        F2S cannot overflow, and it is accurate
  171: *
  172:          F2S = SQRT( ONE+G2 / F2 )
  173: *        Do the F2S(real)*FS(complex) multiply with two real multiplies
  174:          R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
  175:          CS = ONE / F2S
  176:          D = F2 + G2
  177: *        Do complex/real division explicitly with two real divisions
  178:          SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
  179:          SN = SN*DCONJG( GS )
  180:          IF( COUNT.NE.0 ) THEN
  181:             IF( COUNT.GT.0 ) THEN
  182:                DO 30 I = 1, COUNT
  183:                   R = R*SAFMX2
  184:    30          CONTINUE
  185:             ELSE
  186:                DO 40 I = 1, -COUNT
  187:                   R = R*SAFMN2
  188:    40          CONTINUE
  189:             END IF
  190:          END IF
  191:       END IF
  192:       RETURN
  193: *
  194: *     End of ZLARTG
  195: *
  196:       END