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

1.9     ! bertrand    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 November 2011
        !           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 factorization is obtained by Householder's method.  The kth
        !           134: *>  transformation matrix, Z( k ), which is used to introduce zeros into
        !           135: *>  the ( m - k + 1 )th row of A, is given in the form
        !           136: *>
        !           137: *>     Z( k ) = ( I     0   ),
        !           138: *>              ( 0  T( k ) )
        !           139: *>
        !           140: *>  where
        !           141: *>
        !           142: *>     T( k ) = I - tau*u( k )*u( k )**H,   u( k ) = (   1    ),
        !           143: *>                                                 (   0    )
        !           144: *>                                                 ( z( k ) )
        !           145: *>
        !           146: *>  tau is a scalar and z( k ) is an ( n - m ) element vector.
        !           147: *>  tau and z( k ) are chosen to annihilate the elements of the kth row
        !           148: *>  of X.
        !           149: *>
        !           150: *>  The scalar tau is returned in the kth element of TAU and the vector
        !           151: *>  u( k ) in the kth row of A, such that the elements of z( k ) are
        !           152: *>  in  a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
        !           153: *>  the upper triangular part of A.
        !           154: *>
        !           155: *>  Z is given by
        !           156: *>
        !           157: *>     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
        !           158: *> \endverbatim
        !           159: *>
        !           160: *  =====================================================================
1.1       bertrand  161:       SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
                    162: *
1.9     ! bertrand  163: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  164: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    165: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9     ! bertrand  166: *     November 2011
1.1       bertrand  167: *
                    168: *     .. Scalar Arguments ..
                    169:       INTEGER            INFO, LDA, LWORK, M, N
                    170: *     ..
                    171: *     .. Array Arguments ..
                    172:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
                    173: *     ..
                    174: *
                    175: *  =====================================================================
                    176: *
                    177: *     .. Parameters ..
                    178:       COMPLEX*16         ZERO
                    179:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
                    180: *     ..
                    181: *     .. Local Scalars ..
                    182:       LOGICAL            LQUERY
1.8       bertrand  183:       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
                    184:      $                   M1, MU, NB, NBMIN, NX
1.1       bertrand  185: *     ..
                    186: *     .. External Subroutines ..
                    187:       EXTERNAL           XERBLA, ZLARZB, ZLARZT, ZLATRZ
                    188: *     ..
                    189: *     .. Intrinsic Functions ..
                    190:       INTRINSIC          MAX, MIN
                    191: *     ..
                    192: *     .. External Functions ..
                    193:       INTEGER            ILAENV
                    194:       EXTERNAL           ILAENV
                    195: *     ..
                    196: *     .. Executable Statements ..
                    197: *
                    198: *     Test the input arguments
                    199: *
                    200:       INFO = 0
                    201:       LQUERY = ( LWORK.EQ.-1 )
                    202:       IF( M.LT.0 ) THEN
                    203:          INFO = -1
                    204:       ELSE IF( N.LT.M ) THEN
                    205:          INFO = -2
                    206:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
                    207:          INFO = -4
                    208:       END IF
                    209: *
                    210:       IF( INFO.EQ.0 ) THEN
                    211:          IF( M.EQ.0 .OR. M.EQ.N ) THEN
                    212:             LWKOPT = 1
1.8       bertrand  213:             LWKMIN = 1
1.1       bertrand  214:          ELSE
                    215: *
                    216: *           Determine the block size.
                    217: *
                    218:             NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
                    219:             LWKOPT = M*NB
1.8       bertrand  220:             LWKMIN = MAX( 1, M )
1.1       bertrand  221:          END IF
                    222:          WORK( 1 ) = LWKOPT
                    223: *
1.8       bertrand  224:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
1.1       bertrand  225:             INFO = -7
                    226:          END IF
                    227:       END IF
                    228: *
                    229:       IF( INFO.NE.0 ) THEN
                    230:          CALL XERBLA( 'ZTZRZF', -INFO )
                    231:          RETURN
                    232:       ELSE IF( LQUERY ) THEN
                    233:          RETURN
                    234:       END IF
                    235: *
                    236: *     Quick return if possible
                    237: *
                    238:       IF( M.EQ.0 ) THEN
                    239:          RETURN
                    240:       ELSE IF( M.EQ.N ) THEN
                    241:          DO 10 I = 1, N
                    242:             TAU( I ) = ZERO
                    243:    10    CONTINUE
                    244:          RETURN
                    245:       END IF
                    246: *
                    247:       NBMIN = 2
                    248:       NX = 1
                    249:       IWS = M
                    250:       IF( NB.GT.1 .AND. NB.LT.M ) THEN
                    251: *
                    252: *        Determine when to cross over from blocked to unblocked code.
                    253: *
                    254:          NX = MAX( 0, ILAENV( 3, 'ZGERQF', ' ', M, N, -1, -1 ) )
                    255:          IF( NX.LT.M ) THEN
                    256: *
                    257: *           Determine if workspace is large enough for blocked code.
                    258: *
                    259:             LDWORK = M
                    260:             IWS = LDWORK*NB
                    261:             IF( LWORK.LT.IWS ) THEN
                    262: *
                    263: *              Not enough workspace to use optimal NB:  reduce NB and
                    264: *              determine the minimum value of NB.
                    265: *
                    266:                NB = LWORK / LDWORK
                    267:                NBMIN = MAX( 2, ILAENV( 2, 'ZGERQF', ' ', M, N, -1,
                    268:      $                 -1 ) )
                    269:             END IF
                    270:          END IF
                    271:       END IF
                    272: *
                    273:       IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
                    274: *
                    275: *        Use blocked code initially.
                    276: *        The last kk rows are handled by the block method.
                    277: *
                    278:          M1 = MIN( M+1, N )
                    279:          KI = ( ( M-NX-1 ) / NB )*NB
                    280:          KK = MIN( M, KI+NB )
                    281: *
                    282:          DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
                    283:             IB = MIN( M-I+1, NB )
                    284: *
                    285: *           Compute the TZ factorization of the current block
                    286: *           A(i:i+ib-1,i:n)
                    287: *
                    288:             CALL ZLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
                    289:      $                   WORK )
                    290:             IF( I.GT.1 ) THEN
                    291: *
                    292: *              Form the triangular factor of the block reflector
                    293: *              H = H(i+ib-1) . . . H(i+1) H(i)
                    294: *
                    295:                CALL ZLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
                    296:      $                      LDA, TAU( I ), WORK, LDWORK )
                    297: *
                    298: *              Apply H to A(1:i-1,i:n) from the right
                    299: *
                    300:                CALL ZLARZB( 'Right', 'No transpose', 'Backward',
                    301:      $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
                    302:      $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
                    303:      $                      WORK( IB+1 ), LDWORK )
                    304:             END IF
                    305:    20    CONTINUE
                    306:          MU = I + NB - 1
                    307:       ELSE
                    308:          MU = M
                    309:       END IF
                    310: *
                    311: *     Use unblocked code to factor the last or only block
                    312: *
                    313:       IF( MU.GT.0 )
                    314:      $   CALL ZLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
                    315: *
                    316:       WORK( 1 ) = LWKOPT
                    317: *
                    318:       RETURN
                    319: *
                    320: *     End of ZTZRZF
                    321: *
                    322:       END

CVSweb interface <joel.bertrand@systella.fr>