File:  [local] / rpl / lapack / lapack / zlarfb.f
Revision 1.14: download - view: text, annotated - select for diffs - revision graph
Mon Jan 27 09:24:36 2014 UTC (10 years, 4 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Mise à jour de lapack vers la version 3.5.0.

    1: *> \brief \b ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZLARFB + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfb.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfb.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfb.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
   22: *                          T, LDT, C, LDC, WORK, LDWORK )
   23:    24: *       .. Scalar Arguments ..
   25: *       CHARACTER          DIRECT, SIDE, STOREV, TRANS
   26: *       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
   30: *      $                   WORK( LDWORK, * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> ZLARFB applies a complex block reflector H or its transpose H**H to a
   40: *> complex M-by-N matrix C, from either the left or the right.
   41: *> \endverbatim
   42: *
   43: *  Arguments:
   44: *  ==========
   45: *
   46: *> \param[in] SIDE
   47: *> \verbatim
   48: *>          SIDE is CHARACTER*1
   49: *>          = 'L': apply H or H**H from the Left
   50: *>          = 'R': apply H or H**H from the Right
   51: *> \endverbatim
   52: *>
   53: *> \param[in] TRANS
   54: *> \verbatim
   55: *>          TRANS is CHARACTER*1
   56: *>          = 'N': apply H (No transpose)
   57: *>          = 'C': apply H**H (Conjugate transpose)
   58: *> \endverbatim
   59: *>
   60: *> \param[in] DIRECT
   61: *> \verbatim
   62: *>          DIRECT is CHARACTER*1
   63: *>          Indicates how H is formed from a product of elementary
   64: *>          reflectors
   65: *>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
   66: *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
   67: *> \endverbatim
   68: *>
   69: *> \param[in] STOREV
   70: *> \verbatim
   71: *>          STOREV is CHARACTER*1
   72: *>          Indicates how the vectors which define the elementary
   73: *>          reflectors are stored:
   74: *>          = 'C': Columnwise
   75: *>          = 'R': Rowwise
   76: *> \endverbatim
   77: *>
   78: *> \param[in] M
   79: *> \verbatim
   80: *>          M is INTEGER
   81: *>          The number of rows of the matrix C.
   82: *> \endverbatim
   83: *>
   84: *> \param[in] N
   85: *> \verbatim
   86: *>          N is INTEGER
   87: *>          The number of columns of the matrix C.
   88: *> \endverbatim
   89: *>
   90: *> \param[in] K
   91: *> \verbatim
   92: *>          K is INTEGER
   93: *>          The order of the matrix T (= the number of elementary
   94: *>          reflectors whose product defines the block reflector).
   95: *> \endverbatim
   96: *>
   97: *> \param[in] V
   98: *> \verbatim
   99: *>          V is COMPLEX*16 array, dimension
  100: *>                                (LDV,K) if STOREV = 'C'
  101: *>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
  102: *>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
  103: *>          See Further Details.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] LDV
  107: *> \verbatim
  108: *>          LDV is INTEGER
  109: *>          The leading dimension of the array V.
  110: *>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
  111: *>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
  112: *>          if STOREV = 'R', LDV >= K.
  113: *> \endverbatim
  114: *>
  115: *> \param[in] T
  116: *> \verbatim
  117: *>          T is COMPLEX*16 array, dimension (LDT,K)
  118: *>          The triangular K-by-K matrix T in the representation of the
  119: *>          block reflector.
  120: *> \endverbatim
  121: *>
  122: *> \param[in] LDT
  123: *> \verbatim
  124: *>          LDT is INTEGER
  125: *>          The leading dimension of the array T. LDT >= K.
  126: *> \endverbatim
  127: *>
  128: *> \param[in,out] C
  129: *> \verbatim
  130: *>          C is COMPLEX*16 array, dimension (LDC,N)
  131: *>          On entry, the M-by-N matrix C.
  132: *>          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
  133: *> \endverbatim
  134: *>
  135: *> \param[in] LDC
  136: *> \verbatim
  137: *>          LDC is INTEGER
  138: *>          The leading dimension of the array C. LDC >= max(1,M).
  139: *> \endverbatim
  140: *>
  141: *> \param[out] WORK
  142: *> \verbatim
  143: *>          WORK is COMPLEX*16 array, dimension (LDWORK,K)
  144: *> \endverbatim
  145: *>
  146: *> \param[in] LDWORK
  147: *> \verbatim
  148: *>          LDWORK is INTEGER
  149: *>          The leading dimension of the array WORK.
  150: *>          If SIDE = 'L', LDWORK >= max(1,N);
  151: *>          if SIDE = 'R', LDWORK >= max(1,M).
  152: *> \endverbatim
  153: *
  154: *  Authors:
  155: *  ========
  156: *
  157: *> \author Univ. of Tennessee 
  158: *> \author Univ. of California Berkeley 
  159: *> \author Univ. of Colorado Denver 
  160: *> \author NAG Ltd. 
  161: *
  162: *> \date June 2013
  163: *
  164: *> \ingroup complex16OTHERauxiliary
  165: *
  166: *> \par Further Details:
  167: *  =====================
  168: *>
  169: *> \verbatim
  170: *>
  171: *>  The shape of the matrix V and the storage of the vectors which define
  172: *>  the H(i) is best illustrated by the following example with n = 5 and
  173: *>  k = 3. The elements equal to 1 are not stored; the corresponding
  174: *>  array elements are modified but restored on exit. The rest of the
  175: *>  array is not used.
  176: *>
  177: *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
  178: *>
  179: *>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
  180: *>                   ( v1  1    )                     (     1 v2 v2 v2 )
  181: *>                   ( v1 v2  1 )                     (        1 v3 v3 )
  182: *>                   ( v1 v2 v3 )
  183: *>                   ( v1 v2 v3 )
  184: *>
  185: *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
  186: *>
  187: *>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
  188: *>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
  189: *>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
  190: *>                   (     1 v3 )
  191: *>                   (        1 )
  192: *> \endverbatim
  193: *>
  194: *  =====================================================================
  195:       SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
  196:      $                   T, LDT, C, LDC, WORK, LDWORK )
  197: *
  198: *  -- LAPACK auxiliary routine (version 3.5.0) --
  199: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  200: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  201: *     June 2013
  202: *
  203: *     .. Scalar Arguments ..
  204:       CHARACTER          DIRECT, SIDE, STOREV, TRANS
  205:       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
  206: *     ..
  207: *     .. Array Arguments ..
  208:       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
  209:      $                   WORK( LDWORK, * )
  210: *     ..
  211: *
  212: *  =====================================================================
  213: *
  214: *     .. Parameters ..
  215:       COMPLEX*16         ONE
  216:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
  217: *     ..
  218: *     .. Local Scalars ..
  219:       CHARACTER          TRANST
  220:       INTEGER            I, J
  221: *     ..
  222: *     .. External Functions ..
  223:       LOGICAL            LSAME
  224:       EXTERNAL           LSAME
  225: *     ..
  226: *     .. External Subroutines ..
  227:       EXTERNAL           ZCOPY, ZGEMM, ZLACGV, ZTRMM
  228: *     ..
  229: *     .. Intrinsic Functions ..
  230:       INTRINSIC          DCONJG
  231: *     ..
  232: *     .. Executable Statements ..
  233: *
  234: *     Quick return if possible
  235: *
  236:       IF( M.LE.0 .OR. N.LE.0 )
  237:      $   RETURN
  238: *
  239:       IF( LSAME( TRANS, 'N' ) ) THEN
  240:          TRANST = 'C'
  241:       ELSE
  242:          TRANST = 'N'
  243:       END IF
  244: *
  245:       IF( LSAME( STOREV, 'C' ) ) THEN
  246: *
  247:          IF( LSAME( DIRECT, 'F' ) ) THEN
  248: *
  249: *           Let  V =  ( V1 )    (first K rows)
  250: *                     ( V2 )
  251: *           where  V1  is unit lower triangular.
  252: *
  253:             IF( LSAME( SIDE, 'L' ) ) THEN
  254: *
  255: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  256: *                                                    ( C2 )
  257: *
  258: *              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
  259: *
  260: *              W := C1**H
  261: *
  262:                DO 10 J = 1, K
  263:                   CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
  264:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  265:    10          CONTINUE
  266: *
  267: *              W := W * V1
  268: *
  269:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
  270:      $                     K, ONE, V, LDV, WORK, LDWORK )
  271:                IF( M.GT.K ) THEN
  272: *
  273: *                 W := W + C2**H * V2
  274: *
  275:                   CALL ZGEMM( 'Conjugate transpose', 'No transpose', N,
  276:      $                        K, M-K, ONE, C( K+1, 1 ), LDC,
  277:      $                        V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
  278:                END IF
  279: *
  280: *              W := W * T**H  or  W * T
  281: *
  282:                CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
  283:      $                     ONE, T, LDT, WORK, LDWORK )
  284: *
  285: *              C := C - V * W**H
  286: *
  287:                IF( M.GT.K ) THEN
  288: *
  289: *                 C2 := C2 - V2 * W**H
  290: *
  291:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose',
  292:      $                        M-K, N, K, -ONE, V( K+1, 1 ), LDV, WORK,
  293:      $                        LDWORK, ONE, C( K+1, 1 ), LDC )
  294:                END IF
  295: *
  296: *              W := W * V1**H
  297: *
  298:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  299:      $                     'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
  300: *
  301: *              C1 := C1 - W**H
  302: *
  303:                DO 30 J = 1, K
  304:                   DO 20 I = 1, N
  305:                      C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
  306:    20             CONTINUE
  307:    30          CONTINUE
  308: *
  309:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  310: *
  311: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  312: *
  313: *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
  314: *
  315: *              W := C1
  316: *
  317:                DO 40 J = 1, K
  318:                   CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
  319:    40          CONTINUE
  320: *
  321: *              W := W * V1
  322: *
  323:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
  324:      $                     K, ONE, V, LDV, WORK, LDWORK )
  325:                IF( N.GT.K ) THEN
  326: *
  327: *                 W := W + C2 * V2
  328: *
  329:                   CALL ZGEMM( 'No transpose', 'No transpose', M, K, N-K,
  330:      $                        ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
  331:      $                        ONE, WORK, LDWORK )
  332:                END IF
  333: *
  334: *              W := W * T  or  W * T**H
  335: *
  336:                CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
  337:      $                     ONE, T, LDT, WORK, LDWORK )
  338: *
  339: *              C := C - W * V**H
  340: *
  341:                IF( N.GT.K ) THEN
  342: *
  343: *                 C2 := C2 - W * V2**H
  344: *
  345:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  346:      $                        N-K, K, -ONE, WORK, LDWORK, V( K+1, 1 ),
  347:      $                        LDV, ONE, C( 1, K+1 ), LDC )
  348:                END IF
  349: *
  350: *              W := W * V1**H
  351: *
  352:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  353:      $                     'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
  354: *
  355: *              C1 := C1 - W
  356: *
  357:                DO 60 J = 1, K
  358:                   DO 50 I = 1, M
  359:                      C( I, J ) = C( I, J ) - WORK( I, J )
  360:    50             CONTINUE
  361:    60          CONTINUE
  362:             END IF
  363: *
  364:          ELSE
  365: *
  366: *           Let  V =  ( V1 )
  367: *                     ( V2 )    (last K rows)
  368: *           where  V2  is unit upper triangular.
  369: *
  370:             IF( LSAME( SIDE, 'L' ) ) THEN
  371: *
  372: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  373: *                                                    ( C2 )
  374: *
  375: *              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
  376: *
  377: *              W := C2**H
  378: *
  379:                DO 70 J = 1, K
  380:                   CALL ZCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
  381:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  382:    70          CONTINUE
  383: *
  384: *              W := W * V2
  385: *
  386:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
  387:      $                     K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
  388:                IF( M.GT.K ) THEN
  389: *
  390: *                 W := W + C1**H * V1
  391: *
  392:                   CALL ZGEMM( 'Conjugate transpose', 'No transpose', N,
  393:      $                        K, M-K, ONE, C, LDC, V, LDV, ONE, WORK,
  394:      $                        LDWORK )
  395:                END IF
  396: *
  397: *              W := W * T**H  or  W * T
  398: *
  399:                CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
  400:      $                     ONE, T, LDT, WORK, LDWORK )
  401: *
  402: *              C := C - V * W**H
  403: *
  404:                IF( M.GT.K ) THEN
  405: *
  406: *                 C1 := C1 - V1 * W**H
  407: *
  408:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose',
  409:      $                        M-K, N, K, -ONE, V, LDV, WORK, LDWORK,
  410:      $                        ONE, C, LDC )
  411:                END IF
  412: *
  413: *              W := W * V2**H
  414: *
  415:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  416:      $                     'Unit', N, K, ONE, V( M-K+1, 1 ), LDV, WORK,
  417:      $                     LDWORK )
  418: *
  419: *              C2 := C2 - W**H
  420: *
  421:                DO 90 J = 1, K
  422:                   DO 80 I = 1, N
  423:                      C( M-K+J, I ) = C( M-K+J, I ) -
  424:      $                               DCONJG( WORK( I, J ) )
  425:    80             CONTINUE
  426:    90          CONTINUE
  427: *
  428:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  429: *
  430: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  431: *
  432: *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
  433: *
  434: *              W := C2
  435: *
  436:                DO 100 J = 1, K
  437:                   CALL ZCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
  438:   100          CONTINUE
  439: *
  440: *              W := W * V2
  441: *
  442:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
  443:      $                     K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
  444:                IF( N.GT.K ) THEN
  445: *
  446: *                 W := W + C1 * V1
  447: *
  448:                   CALL ZGEMM( 'No transpose', 'No transpose', M, K, N-K,
  449:      $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
  450:                END IF
  451: *
  452: *              W := W * T  or  W * T**H
  453: *
  454:                CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
  455:      $                     ONE, T, LDT, WORK, LDWORK )
  456: *
  457: *              C := C - W * V**H
  458: *
  459:                IF( N.GT.K ) THEN
  460: *
  461: *                 C1 := C1 - W * V1**H
  462: *
  463:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  464:      $                        N-K, K, -ONE, WORK, LDWORK, V, LDV, ONE,
  465:      $                        C, LDC )
  466:                END IF
  467: *
  468: *              W := W * V2**H
  469: *
  470:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  471:      $                     'Unit', M, K, ONE, V( N-K+1, 1 ), LDV, WORK,
  472:      $                     LDWORK )
  473: *
  474: *              C2 := C2 - W
  475: *
  476:                DO 120 J = 1, K
  477:                   DO 110 I = 1, M
  478:                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
  479:   110             CONTINUE
  480:   120          CONTINUE
  481:             END IF
  482:          END IF
  483: *
  484:       ELSE IF( LSAME( STOREV, 'R' ) ) THEN
  485: *
  486:          IF( LSAME( DIRECT, 'F' ) ) THEN
  487: *
  488: *           Let  V =  ( V1  V2 )    (V1: first K columns)
  489: *           where  V1  is unit upper triangular.
  490: *
  491:             IF( LSAME( SIDE, 'L' ) ) THEN
  492: *
  493: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  494: *                                                    ( C2 )
  495: *
  496: *              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
  497: *
  498: *              W := C1**H
  499: *
  500:                DO 130 J = 1, K
  501:                   CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
  502:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  503:   130          CONTINUE
  504: *
  505: *              W := W * V1**H
  506: *
  507:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  508:      $                     'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
  509:                IF( M.GT.K ) THEN
  510: *
  511: *                 W := W + C2**H * V2**H
  512: *
  513:                   CALL ZGEMM( 'Conjugate transpose',
  514:      $                        'Conjugate transpose', N, K, M-K, ONE,
  515:      $                        C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
  516:      $                        WORK, LDWORK )
  517:                END IF
  518: *
  519: *              W := W * T**H  or  W * T
  520: *
  521:                CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
  522:      $                     ONE, T, LDT, WORK, LDWORK )
  523: *
  524: *              C := C - V**H * W**H
  525: *
  526:                IF( M.GT.K ) THEN
  527: *
  528: *                 C2 := C2 - V2**H * W**H
  529: *
  530:                   CALL ZGEMM( 'Conjugate transpose',
  531:      $                        'Conjugate transpose', M-K, N, K, -ONE,
  532:      $                        V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
  533:      $                        C( K+1, 1 ), LDC )
  534:                END IF
  535: *
  536: *              W := W * V1
  537: *
  538:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
  539:      $                     K, ONE, V, LDV, WORK, LDWORK )
  540: *
  541: *              C1 := C1 - W**H
  542: *
  543:                DO 150 J = 1, K
  544:                   DO 140 I = 1, N
  545:                      C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
  546:   140             CONTINUE
  547:   150          CONTINUE
  548: *
  549:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  550: *
  551: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  552: *
  553: *              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
  554: *
  555: *              W := C1
  556: *
  557:                DO 160 J = 1, K
  558:                   CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
  559:   160          CONTINUE
  560: *
  561: *              W := W * V1**H
  562: *
  563:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  564:      $                     'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
  565:                IF( N.GT.K ) THEN
  566: *
  567: *                 W := W + C2 * V2**H
  568: *
  569:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  570:      $                        K, N-K, ONE, C( 1, K+1 ), LDC,
  571:      $                        V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
  572:                END IF
  573: *
  574: *              W := W * T  or  W * T**H
  575: *
  576:                CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
  577:      $                     ONE, T, LDT, WORK, LDWORK )
  578: *
  579: *              C := C - W * V
  580: *
  581:                IF( N.GT.K ) THEN
  582: *
  583: *                 C2 := C2 - W * V2
  584: *
  585:                   CALL ZGEMM( 'No transpose', 'No transpose', M, N-K, K,
  586:      $                        -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
  587:      $                        C( 1, K+1 ), LDC )
  588:                END IF
  589: *
  590: *              W := W * V1
  591: *
  592:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
  593:      $                     K, ONE, V, LDV, WORK, LDWORK )
  594: *
  595: *              C1 := C1 - W
  596: *
  597:                DO 180 J = 1, K
  598:                   DO 170 I = 1, M
  599:                      C( I, J ) = C( I, J ) - WORK( I, J )
  600:   170             CONTINUE
  601:   180          CONTINUE
  602: *
  603:             END IF
  604: *
  605:          ELSE
  606: *
  607: *           Let  V =  ( V1  V2 )    (V2: last K columns)
  608: *           where  V2  is unit lower triangular.
  609: *
  610:             IF( LSAME( SIDE, 'L' ) ) THEN
  611: *
  612: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  613: *                                                    ( C2 )
  614: *
  615: *              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
  616: *
  617: *              W := C2**H
  618: *
  619:                DO 190 J = 1, K
  620:                   CALL ZCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
  621:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  622:   190          CONTINUE
  623: *
  624: *              W := W * V2**H
  625: *
  626:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  627:      $                     'Unit', N, K, ONE, V( 1, M-K+1 ), LDV, WORK,
  628:      $                     LDWORK )
  629:                IF( M.GT.K ) THEN
  630: *
  631: *                 W := W + C1**H * V1**H
  632: *
  633:                   CALL ZGEMM( 'Conjugate transpose',
  634:      $                        'Conjugate transpose', N, K, M-K, ONE, C,
  635:      $                        LDC, V, LDV, ONE, WORK, LDWORK )
  636:                END IF
  637: *
  638: *              W := W * T**H  or  W * T
  639: *
  640:                CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
  641:      $                     ONE, T, LDT, WORK, LDWORK )
  642: *
  643: *              C := C - V**H * W**H
  644: *
  645:                IF( M.GT.K ) THEN
  646: *
  647: *                 C1 := C1 - V1**H * W**H
  648: *
  649:                   CALL ZGEMM( 'Conjugate transpose',
  650:      $                        'Conjugate transpose', M-K, N, K, -ONE, V,
  651:      $                        LDV, WORK, LDWORK, ONE, C, LDC )
  652:                END IF
  653: *
  654: *              W := W * V2
  655: *
  656:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
  657:      $                     K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
  658: *
  659: *              C2 := C2 - W**H
  660: *
  661:                DO 210 J = 1, K
  662:                   DO 200 I = 1, N
  663:                      C( M-K+J, I ) = C( M-K+J, I ) -
  664:      $                               DCONJG( WORK( I, J ) )
  665:   200             CONTINUE
  666:   210          CONTINUE
  667: *
  668:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  669: *
  670: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  671: *
  672: *              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
  673: *
  674: *              W := C2
  675: *
  676:                DO 220 J = 1, K
  677:                   CALL ZCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
  678:   220          CONTINUE
  679: *
  680: *              W := W * V2**H
  681: *
  682:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  683:      $                     'Unit', M, K, ONE, V( 1, N-K+1 ), LDV, WORK,
  684:      $                     LDWORK )
  685:                IF( N.GT.K ) THEN
  686: *
  687: *                 W := W + C1 * V1**H
  688: *
  689:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  690:      $                        K, N-K, ONE, C, LDC, V, LDV, ONE, WORK,
  691:      $                        LDWORK )
  692:                END IF
  693: *
  694: *              W := W * T  or  W * T**H
  695: *
  696:                CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
  697:      $                     ONE, T, LDT, WORK, LDWORK )
  698: *
  699: *              C := C - W * V
  700: *
  701:                IF( N.GT.K ) THEN
  702: *
  703: *                 C1 := C1 - W * V1
  704: *
  705:                   CALL ZGEMM( 'No transpose', 'No transpose', M, N-K, K,
  706:      $                        -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
  707:                END IF
  708: *
  709: *              W := W * V2
  710: *
  711:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
  712:      $                     K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
  713: *
  714: *              C1 := C1 - W
  715: *
  716:                DO 240 J = 1, K
  717:                   DO 230 I = 1, M
  718:                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
  719:   230             CONTINUE
  720:   240          CONTINUE
  721: *
  722:             END IF
  723: *
  724:          END IF
  725:       END IF
  726: *
  727:       RETURN
  728: *
  729: *     End of ZLARFB
  730: *
  731:       END

CVSweb interface <joel.bertrand@systella.fr>