Annotation of rpl/lapack/lapack/ztzrzf.f, revision 1.17

1.9       bertrand    1: *> \brief \b ZTZRZF
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.16      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.16      bertrand    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">
1.9       bertrand   15: *> [TXT]</a>
1.16      bertrand   16: *> \endhtmlonly
1.9       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
1.16      bertrand   22: *
1.9       bertrand   23: *       .. Scalar Arguments ..
                     24: *       INTEGER            INFO, LDA, LWORK, M, N
                     25: *       ..
                     26: *       .. Array Arguments ..
                     27: *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
                     28: *       ..
1.16      bertrand   29: *
1.9       bertrand   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: *
1.16      bertrand  114: *> \author Univ. of Tennessee
                    115: *> \author Univ. of California Berkeley
                    116: *> \author Univ. of Colorado Denver
                    117: *> \author NAG Ltd.
1.9       bertrand  118: *
1.11      bertrand  119: *> \date April 2012
1.9       bertrand  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: *>
1.11      bertrand  133: *>  The N-by-N matrix Z can be computed by
1.9       bertrand  134: *>
1.11      bertrand  135: *>     Z =  Z(1)*Z(2)* ... *Z(M)
1.9       bertrand  136: *>
1.11      bertrand  137: *>  where each N-by-N Z(k) is given by
1.9       bertrand  138: *>
1.11      bertrand  139: *>     Z(k) = I - tau(k)*v(k)*v(k)**H
1.9       bertrand  140: *>
1.11      bertrand  141: *>  with v(k) is the kth row vector of the M-by-N matrix
1.9       bertrand  142: *>
1.11      bertrand  143: *>     V = ( I   A(:,M+1:N) )
1.9       bertrand  144: *>
1.16      bertrand  145: *>  I is the M-by-M identity matrix, A(:,M+1:N)
1.11      bertrand  146: *>  is the output stored in A on exit from DTZRZF,
                    147: *>  and tau(k) is the kth element of the array TAU.
1.9       bertrand  148: *>
                    149: *> \endverbatim
                    150: *>
                    151: *  =====================================================================
1.1       bertrand  152:       SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
                    153: *
1.16      bertrand  154: *  -- LAPACK computational routine (version 3.7.0) --
1.1       bertrand  155: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    156: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11      bertrand  157: *     April 2012
1.1       bertrand  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
1.8       bertrand  174:       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
                    175:      $                   M1, MU, NB, NBMIN, NX
1.1       bertrand  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
1.8       bertrand  204:             LWKMIN = 1
1.1       bertrand  205:          ELSE
                    206: *
                    207: *           Determine the block size.
                    208: *
                    209:             NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
                    210:             LWKOPT = M*NB
1.8       bertrand  211:             LWKMIN = MAX( 1, M )
1.1       bertrand  212:          END IF
                    213:          WORK( 1 ) = LWKOPT
                    214: *
1.8       bertrand  215:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
1.1       bertrand  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

CVSweb interface <joel.bertrand@systella.fr>