Return to zlargv.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

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