File:  [local] / rpl / lapack / lapack / zunmrz.f
Revision 1.9: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 20:43:24 2011 UTC (12 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Mise à jour de Lapack.

    1: *> \brief \b ZUNMRZ
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZUNMRZ + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunmrz.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunmrz.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunmrz.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZUNMRZ( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
   22: *                          WORK, LWORK, INFO )
   23:    24: *       .. Scalar Arguments ..
   25: *       CHARACTER          SIDE, TRANS
   26: *       INTEGER            INFO, K, L, LDA, LDC, LWORK, M, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       COMPLEX*16         A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> ZUNMRZ overwrites the general complex M-by-N matrix C with
   39: *>
   40: *>                 SIDE = 'L'     SIDE = 'R'
   41: *> TRANS = 'N':      Q * C          C * Q
   42: *> TRANS = 'C':      Q**H * C       C * Q**H
   43: *>
   44: *> where Q is a complex unitary matrix defined as the product of k
   45: *> elementary reflectors
   46: *>
   47: *>       Q = H(1) H(2) . . . H(k)
   48: *>
   49: *> as returned by ZTZRZF. Q is of order M if SIDE = 'L' and of order N
   50: *> if SIDE = 'R'.
   51: *> \endverbatim
   52: *
   53: *  Arguments:
   54: *  ==========
   55: *
   56: *> \param[in] SIDE
   57: *> \verbatim
   58: *>          SIDE is CHARACTER*1
   59: *>          = 'L': apply Q or Q**H from the Left;
   60: *>          = 'R': apply Q or Q**H from the Right.
   61: *> \endverbatim
   62: *>
   63: *> \param[in] TRANS
   64: *> \verbatim
   65: *>          TRANS is CHARACTER*1
   66: *>          = 'N':  No transpose, apply Q;
   67: *>          = 'C':  Conjugate transpose, apply Q**H.
   68: *> \endverbatim
   69: *>
   70: *> \param[in] M
   71: *> \verbatim
   72: *>          M is INTEGER
   73: *>          The number of rows of the matrix C. M >= 0.
   74: *> \endverbatim
   75: *>
   76: *> \param[in] N
   77: *> \verbatim
   78: *>          N is INTEGER
   79: *>          The number of columns of the matrix C. N >= 0.
   80: *> \endverbatim
   81: *>
   82: *> \param[in] K
   83: *> \verbatim
   84: *>          K is INTEGER
   85: *>          The number of elementary reflectors whose product defines
   86: *>          the matrix Q.
   87: *>          If SIDE = 'L', M >= K >= 0;
   88: *>          if SIDE = 'R', N >= K >= 0.
   89: *> \endverbatim
   90: *>
   91: *> \param[in] L
   92: *> \verbatim
   93: *>          L is INTEGER
   94: *>          The number of columns of the matrix A containing
   95: *>          the meaningful part of the Householder reflectors.
   96: *>          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] A
  100: *> \verbatim
  101: *>          A is COMPLEX*16 array, dimension
  102: *>                               (LDA,M) if SIDE = 'L',
  103: *>                               (LDA,N) if SIDE = 'R'
  104: *>          The i-th row must contain the vector which defines the
  105: *>          elementary reflector H(i), for i = 1,2,...,k, as returned by
  106: *>          ZTZRZF in the last k rows of its array argument A.
  107: *>          A is modified by the routine but restored on exit.
  108: *> \endverbatim
  109: *>
  110: *> \param[in] LDA
  111: *> \verbatim
  112: *>          LDA is INTEGER
  113: *>          The leading dimension of the array A. LDA >= max(1,K).
  114: *> \endverbatim
  115: *>
  116: *> \param[in] TAU
  117: *> \verbatim
  118: *>          TAU is COMPLEX*16 array, dimension (K)
  119: *>          TAU(i) must contain the scalar factor of the elementary
  120: *>          reflector H(i), as returned by ZTZRZF.
  121: *> \endverbatim
  122: *>
  123: *> \param[in,out] C
  124: *> \verbatim
  125: *>          C is COMPLEX*16 array, dimension (LDC,N)
  126: *>          On entry, the M-by-N matrix C.
  127: *>          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
  128: *> \endverbatim
  129: *>
  130: *> \param[in] LDC
  131: *> \verbatim
  132: *>          LDC is INTEGER
  133: *>          The leading dimension of the array C. LDC >= max(1,M).
  134: *> \endverbatim
  135: *>
  136: *> \param[out] WORK
  137: *> \verbatim
  138: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
  139: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  140: *> \endverbatim
  141: *>
  142: *> \param[in] LWORK
  143: *> \verbatim
  144: *>          LWORK is INTEGER
  145: *>          The dimension of the array WORK.
  146: *>          If SIDE = 'L', LWORK >= max(1,N);
  147: *>          if SIDE = 'R', LWORK >= max(1,M).
  148: *>          For optimum performance LWORK >= N*NB if SIDE = 'L', and
  149: *>          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
  150: *>          blocksize.
  151: *>
  152: *>          If LWORK = -1, then a workspace query is assumed; the routine
  153: *>          only calculates the optimal size of the WORK array, returns
  154: *>          this value as the first entry of the WORK array, and no error
  155: *>          message related to LWORK is issued by XERBLA.
  156: *> \endverbatim
  157: *>
  158: *> \param[out] INFO
  159: *> \verbatim
  160: *>          INFO is INTEGER
  161: *>          = 0:  successful exit
  162: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  163: *> \endverbatim
  164: *
  165: *  Authors:
  166: *  ========
  167: *
  168: *> \author Univ. of Tennessee 
  169: *> \author Univ. of California Berkeley 
  170: *> \author Univ. of Colorado Denver 
  171: *> \author NAG Ltd. 
  172: *
  173: *> \date November 2011
  174: *
  175: *> \ingroup complex16OTHERcomputational
  176: *
  177: *> \par Contributors:
  178: *  ==================
  179: *>
  180: *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  181: *
  182: *> \par Further Details:
  183: *  =====================
  184: *>
  185: *> \verbatim
  186: *> \endverbatim
  187: *>
  188: *  =====================================================================
  189:       SUBROUTINE ZUNMRZ( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
  190:      $                   WORK, LWORK, INFO )
  191: *
  192: *  -- LAPACK computational routine (version 3.4.0) --
  193: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  194: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  195: *     November 2011
  196: *
  197: *     .. Scalar Arguments ..
  198:       CHARACTER          SIDE, TRANS
  199:       INTEGER            INFO, K, L, LDA, LDC, LWORK, M, N
  200: *     ..
  201: *     .. Array Arguments ..
  202:       COMPLEX*16         A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
  203: *     ..
  204: *
  205: *  =====================================================================
  206: *
  207: *     .. Parameters ..
  208:       INTEGER            NBMAX, LDT
  209:       PARAMETER          ( NBMAX = 64, LDT = NBMAX+1 )
  210: *     ..
  211: *     .. Local Scalars ..
  212:       LOGICAL            LEFT, LQUERY, NOTRAN
  213:       CHARACTER          TRANST
  214:       INTEGER            I, I1, I2, I3, IB, IC, IINFO, IWS, JA, JC,
  215:      $                   LDWORK, LWKOPT, MI, NB, NBMIN, NI, NQ, NW
  216: *     ..
  217: *     .. Local Arrays ..
  218:       COMPLEX*16         T( LDT, NBMAX )
  219: *     ..
  220: *     .. External Functions ..
  221:       LOGICAL            LSAME
  222:       INTEGER            ILAENV
  223:       EXTERNAL           LSAME, ILAENV
  224: *     ..
  225: *     .. External Subroutines ..
  226:       EXTERNAL           XERBLA, ZLARZB, ZLARZT, ZUNMR3
  227: *     ..
  228: *     .. Intrinsic Functions ..
  229:       INTRINSIC          MAX, MIN
  230: *     ..
  231: *     .. Executable Statements ..
  232: *
  233: *     Test the input arguments
  234: *
  235:       INFO = 0
  236:       LEFT = LSAME( SIDE, 'L' )
  237:       NOTRAN = LSAME( TRANS, 'N' )
  238:       LQUERY = ( LWORK.EQ.-1 )
  239: *
  240: *     NQ is the order of Q and NW is the minimum dimension of WORK
  241: *
  242:       IF( LEFT ) THEN
  243:          NQ = M
  244:          NW = MAX( 1, N )
  245:       ELSE
  246:          NQ = N
  247:          NW = MAX( 1, M )
  248:       END IF
  249:       IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
  250:          INFO = -1
  251:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
  252:          INFO = -2
  253:       ELSE IF( M.LT.0 ) THEN
  254:          INFO = -3
  255:       ELSE IF( N.LT.0 ) THEN
  256:          INFO = -4
  257:       ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
  258:          INFO = -5
  259:       ELSE IF( L.LT.0 .OR. ( LEFT .AND. ( L.GT.M ) ) .OR.
  260:      $         ( .NOT.LEFT .AND. ( L.GT.N ) ) ) THEN
  261:          INFO = -6
  262:       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
  263:          INFO = -8
  264:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  265:          INFO = -11
  266:       END IF
  267: *
  268:       IF( INFO.EQ.0 ) THEN
  269:          IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  270:             LWKOPT = 1
  271:          ELSE
  272: *
  273: *           Determine the block size.  NB may be at most NBMAX, where
  274: *           NBMAX is used to define the local array T.
  275: *
  276:             NB = MIN( NBMAX, ILAENV( 1, 'ZUNMRQ', SIDE // TRANS, M, N,
  277:      $                               K, -1 ) )
  278:             LWKOPT = NW*NB
  279:          END IF
  280:          WORK( 1 ) = LWKOPT
  281: *
  282:          IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
  283:             INFO = -13
  284:          END IF
  285:       END IF
  286: *
  287:       IF( INFO.NE.0 ) THEN
  288:          CALL XERBLA( 'ZUNMRZ', -INFO )
  289:          RETURN
  290:       ELSE IF( LQUERY ) THEN
  291:          RETURN
  292:       END IF
  293: *
  294: *     Quick return if possible
  295: *
  296:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  297:          RETURN
  298:       END IF
  299: *
  300: *     Determine the block size.  NB may be at most NBMAX, where NBMAX
  301: *     is used to define the local array T.
  302: *
  303:       NB = MIN( NBMAX, ILAENV( 1, 'ZUNMRQ', SIDE // TRANS, M, N, K,
  304:      $     -1 ) )
  305:       NBMIN = 2
  306:       LDWORK = NW
  307:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
  308:          IWS = NW*NB
  309:          IF( LWORK.LT.IWS ) THEN
  310:             NB = LWORK / LDWORK
  311:             NBMIN = MAX( 2, ILAENV( 2, 'ZUNMRQ', SIDE // TRANS, M, N, K,
  312:      $              -1 ) )
  313:          END IF
  314:       ELSE
  315:          IWS = NW
  316:       END IF
  317: *
  318:       IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
  319: *
  320: *        Use unblocked code
  321: *
  322:          CALL ZUNMR3( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC,
  323:      $                WORK, IINFO )
  324:       ELSE
  325: *
  326: *        Use blocked code
  327: *
  328:          IF( ( LEFT .AND. .NOT.NOTRAN ) .OR.
  329:      $       ( .NOT.LEFT .AND. NOTRAN ) ) THEN
  330:             I1 = 1
  331:             I2 = K
  332:             I3 = NB
  333:          ELSE
  334:             I1 = ( ( K-1 ) / NB )*NB + 1
  335:             I2 = 1
  336:             I3 = -NB
  337:          END IF
  338: *
  339:          IF( LEFT ) THEN
  340:             NI = N
  341:             JC = 1
  342:             JA = M - L + 1
  343:          ELSE
  344:             MI = M
  345:             IC = 1
  346:             JA = N - L + 1
  347:          END IF
  348: *
  349:          IF( NOTRAN ) THEN
  350:             TRANST = 'C'
  351:          ELSE
  352:             TRANST = 'N'
  353:          END IF
  354: *
  355:          DO 10 I = I1, I2, I3
  356:             IB = MIN( NB, K-I+1 )
  357: *
  358: *           Form the triangular factor of the block reflector
  359: *           H = H(i+ib-1) . . . H(i+1) H(i)
  360: *
  361:             CALL ZLARZT( 'Backward', 'Rowwise', L, IB, A( I, JA ), LDA,
  362:      $                   TAU( I ), T, LDT )
  363: *
  364:             IF( LEFT ) THEN
  365: *
  366: *              H or H**H is applied to C(i:m,1:n)
  367: *
  368:                MI = M - I + 1
  369:                IC = I
  370:             ELSE
  371: *
  372: *              H or H**H is applied to C(1:m,i:n)
  373: *
  374:                NI = N - I + 1
  375:                JC = I
  376:             END IF
  377: *
  378: *           Apply H or H**H
  379: *
  380:             CALL ZLARZB( SIDE, TRANST, 'Backward', 'Rowwise', MI, NI,
  381:      $                   IB, L, A( I, JA ), LDA, T, LDT, C( IC, JC ),
  382:      $                   LDC, WORK, LDWORK )
  383:    10    CONTINUE
  384: *
  385:       END IF
  386: *
  387:       WORK( 1 ) = LWKOPT
  388: *
  389:       RETURN
  390: *
  391: *     End of ZUNMRZ
  392: *
  393:       END

CVSweb interface <joel.bertrand@systella.fr>