Return to ztzrzf.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b ZTZRZF
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZTZRZF + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztzrzf.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztzrzf.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztzrzf.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
   22: * 
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, LDA, LWORK, M, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   28: *       ..
   29: *  
   30: *
   31: *> \par Purpose:
   32: *  =============
   33: *>
   34: *> \verbatim
   35: *>
   36: *> ZTZRZF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A
   37: *> to upper triangular form by means of unitary transformations.
   38: *>
   39: *> The upper trapezoidal matrix A is factored as
   40: *>
   41: *>    A = ( R  0 ) * Z,
   42: *>
   43: *> where Z is an N-by-N unitary matrix and R is an M-by-M upper
   44: *> triangular matrix.
   45: *> \endverbatim
   46: *
   47: *  Arguments:
   48: *  ==========
   49: *
   50: *> \param[in] M
   51: *> \verbatim
   52: *>          M is INTEGER
   53: *>          The number of rows of the matrix A.  M >= 0.
   54: *> \endverbatim
   55: *>
   56: *> \param[in] N
   57: *> \verbatim
   58: *>          N is INTEGER
   59: *>          The number of columns of the matrix A.  N >= M.
   60: *> \endverbatim
   61: *>
   62: *> \param[in,out] A
   63: *> \verbatim
   64: *>          A is COMPLEX*16 array, dimension (LDA,N)
   65: *>          On entry, the leading M-by-N upper trapezoidal part of the
   66: *>          array A must contain the matrix to be factorized.
   67: *>          On exit, the leading M-by-M upper triangular part of A
   68: *>          contains the upper triangular matrix R, and elements M+1 to
   69: *>          N of the first M rows of A, with the array TAU, represent the
   70: *>          unitary matrix Z as a product of M elementary reflectors.
   71: *> \endverbatim
   72: *>
   73: *> \param[in] LDA
   74: *> \verbatim
   75: *>          LDA is INTEGER
   76: *>          The leading dimension of the array A.  LDA >= max(1,M).
   77: *> \endverbatim
   78: *>
   79: *> \param[out] TAU
   80: *> \verbatim
   81: *>          TAU is COMPLEX*16 array, dimension (M)
   82: *>          The scalar factors of the elementary reflectors.
   83: *> \endverbatim
   84: *>
   85: *> \param[out] WORK
   86: *> \verbatim
   87: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
   88: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   89: *> \endverbatim
   90: *>
   91: *> \param[in] LWORK
   92: *> \verbatim
   93: *>          LWORK is INTEGER
   94: *>          The dimension of the array WORK.  LWORK >= max(1,M).
   95: *>          For optimum performance LWORK >= M*NB, where NB is
   96: *>          the optimal blocksize.
   97: *>
   98: *>          If LWORK = -1, then a workspace query is assumed; the routine
   99: *>          only calculates the optimal size of the WORK array, returns
  100: *>          this value as the first entry of the WORK array, and no error
  101: *>          message related to LWORK is issued by XERBLA.
  102: *> \endverbatim
  103: *>
  104: *> \param[out] INFO
  105: *> \verbatim
  106: *>          INFO is INTEGER
  107: *>          = 0:  successful exit
  108: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  109: *> \endverbatim
  110: *
  111: *  Authors:
  112: *  ========
  113: *
  114: *> \author Univ. of Tennessee 
  115: *> \author Univ. of California Berkeley 
  116: *> \author Univ. of Colorado Denver 
  117: *> \author NAG Ltd. 
  118: *
  119: *> \date April 2012
  120: *
  121: *> \ingroup complex16OTHERcomputational
  122: *
  123: *> \par Contributors:
  124: *  ==================
  125: *>
  126: *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  127: *
  128: *> \par Further Details:
  129: *  =====================
  130: *>
  131: *> \verbatim
  132: *>
  133: *>  The N-by-N matrix Z can be computed by
  134: *>
  135: *>     Z =  Z(1)*Z(2)* ... *Z(M)
  136: *>
  137: *>  where each N-by-N Z(k) is given by
  138: *>
  139: *>     Z(k) = I - tau(k)*v(k)*v(k)**H
  140: *>
  141: *>  with v(k) is the kth row vector of the M-by-N matrix
  142: *>
  143: *>     V = ( I   A(:,M+1:N) )
  144: *>
  145: *>  I is the M-by-M identity matrix, A(:,M+1:N) 
  146: *>  is the output stored in A on exit from DTZRZF,
  147: *>  and tau(k) is the kth element of the array TAU.
  148: *>
  149: *> \endverbatim
  150: *>
  151: *  =====================================================================
  152:       SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
  153: *
  154: *  -- LAPACK computational routine (version 3.4.1) --
  155: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  156: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  157: *     April 2012
  158: *
  159: *     .. Scalar Arguments ..
  160:       INTEGER            INFO, LDA, LWORK, M, N
  161: *     ..
  162: *     .. Array Arguments ..
  163:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  164: *     ..
  165: *
  166: *  =====================================================================
  167: *
  168: *     .. Parameters ..
  169:       COMPLEX*16         ZERO
  170:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  171: *     ..
  172: *     .. Local Scalars ..
  173:       LOGICAL            LQUERY
  174:       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
  175:      $                   M1, MU, NB, NBMIN, NX
  176: *     ..
  177: *     .. External Subroutines ..
  178:       EXTERNAL           XERBLA, ZLARZB, ZLARZT, ZLATRZ
  179: *     ..
  180: *     .. Intrinsic Functions ..
  181:       INTRINSIC          MAX, MIN
  182: *     ..
  183: *     .. External Functions ..
  184:       INTEGER            ILAENV
  185:       EXTERNAL           ILAENV
  186: *     ..
  187: *     .. Executable Statements ..
  188: *
  189: *     Test the input arguments
  190: *
  191:       INFO = 0
  192:       LQUERY = ( LWORK.EQ.-1 )
  193:       IF( M.LT.0 ) THEN
  194:          INFO = -1
  195:       ELSE IF( N.LT.M ) THEN
  196:          INFO = -2
  197:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  198:          INFO = -4
  199:       END IF
  200: *
  201:       IF( INFO.EQ.0 ) THEN
  202:          IF( M.EQ.0 .OR. M.EQ.N ) THEN
  203:             LWKOPT = 1
  204:             LWKMIN = 1
  205:          ELSE
  206: *
  207: *           Determine the block size.
  208: *
  209:             NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
  210:             LWKOPT = M*NB
  211:             LWKMIN = MAX( 1, M )
  212:          END IF
  213:          WORK( 1 ) = LWKOPT
  214: *
  215:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
  216:             INFO = -7
  217:          END IF
  218:       END IF
  219: *
  220:       IF( INFO.NE.0 ) THEN
  221:          CALL XERBLA( 'ZTZRZF', -INFO )
  222:          RETURN
  223:       ELSE IF( LQUERY ) THEN
  224:          RETURN
  225:       END IF
  226: *
  227: *     Quick return if possible
  228: *
  229:       IF( M.EQ.0 ) THEN
  230:          RETURN
  231:       ELSE IF( M.EQ.N ) THEN
  232:          DO 10 I = 1, N
  233:             TAU( I ) = ZERO
  234:    10    CONTINUE
  235:          RETURN
  236:       END IF
  237: *
  238:       NBMIN = 2
  239:       NX = 1
  240:       IWS = M
  241:       IF( NB.GT.1 .AND. NB.LT.M ) THEN
  242: *
  243: *        Determine when to cross over from blocked to unblocked code.
  244: *
  245:          NX = MAX( 0, ILAENV( 3, 'ZGERQF', ' ', M, N, -1, -1 ) )
  246:          IF( NX.LT.M ) THEN
  247: *
  248: *           Determine if workspace is large enough for blocked code.
  249: *
  250:             LDWORK = M
  251:             IWS = LDWORK*NB
  252:             IF( LWORK.LT.IWS ) THEN
  253: *
  254: *              Not enough workspace to use optimal NB:  reduce NB and
  255: *              determine the minimum value of NB.
  256: *
  257:                NB = LWORK / LDWORK
  258:                NBMIN = MAX( 2, ILAENV( 2, 'ZGERQF', ' ', M, N, -1,
  259:      $                 -1 ) )
  260:             END IF
  261:          END IF
  262:       END IF
  263: *
  264:       IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
  265: *
  266: *        Use blocked code initially.
  267: *        The last kk rows are handled by the block method.
  268: *
  269:          M1 = MIN( M+1, N )
  270:          KI = ( ( M-NX-1 ) / NB )*NB
  271:          KK = MIN( M, KI+NB )
  272: *
  273:          DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
  274:             IB = MIN( M-I+1, NB )
  275: *
  276: *           Compute the TZ factorization of the current block
  277: *           A(i:i+ib-1,i:n)
  278: *
  279:             CALL ZLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
  280:      $                   WORK )
  281:             IF( I.GT.1 ) THEN
  282: *
  283: *              Form the triangular factor of the block reflector
  284: *              H = H(i+ib-1) . . . H(i+1) H(i)
  285: *
  286:                CALL ZLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
  287:      $                      LDA, TAU( I ), WORK, LDWORK )
  288: *
  289: *              Apply H to A(1:i-1,i:n) from the right
  290: *
  291:                CALL ZLARZB( 'Right', 'No transpose', 'Backward',
  292:      $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
  293:      $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
  294:      $                      WORK( IB+1 ), LDWORK )
  295:             END IF
  296:    20    CONTINUE
  297:          MU = I + NB - 1
  298:       ELSE
  299:          MU = M
  300:       END IF
  301: *
  302: *     Use unblocked code to factor the last or only block
  303: *
  304:       IF( MU.GT.0 )
  305:      $   CALL ZLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
  306: *
  307:       WORK( 1 ) = LWKOPT
  308: *
  309:       RETURN
  310: *
  311: *     End of ZTZRZF
  312: *
  313:       END