File:  [local] / rpl / lapack / lapack / zlarfb.f
Revision 1.21: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:31 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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: *>          If SIDE = 'L', M >= K >= 0;
   96: *>          if SIDE = 'R', N >= K >= 0.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] V
  100: *> \verbatim
  101: *>          V is COMPLEX*16 array, dimension
  102: *>                                (LDV,K) if STOREV = 'C'
  103: *>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
  104: *>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
  105: *>          See Further Details.
  106: *> \endverbatim
  107: *>
  108: *> \param[in] LDV
  109: *> \verbatim
  110: *>          LDV is INTEGER
  111: *>          The leading dimension of the array V.
  112: *>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
  113: *>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
  114: *>          if STOREV = 'R', LDV >= K.
  115: *> \endverbatim
  116: *>
  117: *> \param[in] T
  118: *> \verbatim
  119: *>          T is COMPLEX*16 array, dimension (LDT,K)
  120: *>          The triangular K-by-K matrix T in the representation of the
  121: *>          block reflector.
  122: *> \endverbatim
  123: *>
  124: *> \param[in] LDT
  125: *> \verbatim
  126: *>          LDT is INTEGER
  127: *>          The leading dimension of the array T. LDT >= K.
  128: *> \endverbatim
  129: *>
  130: *> \param[in,out] C
  131: *> \verbatim
  132: *>          C is COMPLEX*16 array, dimension (LDC,N)
  133: *>          On entry, the M-by-N matrix C.
  134: *>          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
  135: *> \endverbatim
  136: *>
  137: *> \param[in] LDC
  138: *> \verbatim
  139: *>          LDC is INTEGER
  140: *>          The leading dimension of the array C. LDC >= max(1,M).
  141: *> \endverbatim
  142: *>
  143: *> \param[out] WORK
  144: *> \verbatim
  145: *>          WORK is COMPLEX*16 array, dimension (LDWORK,K)
  146: *> \endverbatim
  147: *>
  148: *> \param[in] LDWORK
  149: *> \verbatim
  150: *>          LDWORK is INTEGER
  151: *>          The leading dimension of the array WORK.
  152: *>          If SIDE = 'L', LDWORK >= max(1,N);
  153: *>          if SIDE = 'R', LDWORK >= max(1,M).
  154: *> \endverbatim
  155: *
  156: *  Authors:
  157: *  ========
  158: *
  159: *> \author Univ. of Tennessee
  160: *> \author Univ. of California Berkeley
  161: *> \author Univ. of Colorado Denver
  162: *> \author NAG Ltd.
  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 --
  199: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  200: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  201: *
  202: *     .. Scalar Arguments ..
  203:       CHARACTER          DIRECT, SIDE, STOREV, TRANS
  204:       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
  205: *     ..
  206: *     .. Array Arguments ..
  207:       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
  208:      $                   WORK( LDWORK, * )
  209: *     ..
  210: *
  211: *  =====================================================================
  212: *
  213: *     .. Parameters ..
  214:       COMPLEX*16         ONE
  215:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
  216: *     ..
  217: *     .. Local Scalars ..
  218:       CHARACTER          TRANST
  219:       INTEGER            I, J
  220: *     ..
  221: *     .. External Functions ..
  222:       LOGICAL            LSAME
  223:       EXTERNAL           LSAME
  224: *     ..
  225: *     .. External Subroutines ..
  226:       EXTERNAL           ZCOPY, ZGEMM, ZLACGV, ZTRMM
  227: *     ..
  228: *     .. Intrinsic Functions ..
  229:       INTRINSIC          DCONJG
  230: *     ..
  231: *     .. Executable Statements ..
  232: *
  233: *     Quick return if possible
  234: *
  235:       IF( M.LE.0 .OR. N.LE.0 )
  236:      $   RETURN
  237: *
  238:       IF( LSAME( TRANS, 'N' ) ) THEN
  239:          TRANST = 'C'
  240:       ELSE
  241:          TRANST = 'N'
  242:       END IF
  243: *
  244:       IF( LSAME( STOREV, 'C' ) ) THEN
  245: *
  246:          IF( LSAME( DIRECT, 'F' ) ) THEN
  247: *
  248: *           Let  V =  ( V1 )    (first K rows)
  249: *                     ( V2 )
  250: *           where  V1  is unit lower triangular.
  251: *
  252:             IF( LSAME( SIDE, 'L' ) ) THEN
  253: *
  254: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  255: *                                                    ( C2 )
  256: *
  257: *              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
  258: *
  259: *              W := C1**H
  260: *
  261:                DO 10 J = 1, K
  262:                   CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
  263:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  264:    10          CONTINUE
  265: *
  266: *              W := W * V1
  267: *
  268:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
  269:      $                     K, ONE, V, LDV, WORK, LDWORK )
  270:                IF( M.GT.K ) THEN
  271: *
  272: *                 W := W + C2**H * V2
  273: *
  274:                   CALL ZGEMM( 'Conjugate transpose', 'No transpose', N,
  275:      $                        K, M-K, ONE, C( K+1, 1 ), LDC,
  276:      $                        V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
  277:                END IF
  278: *
  279: *              W := W * T**H  or  W * T
  280: *
  281:                CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
  282:      $                     ONE, T, LDT, WORK, LDWORK )
  283: *
  284: *              C := C - V * W**H
  285: *
  286:                IF( M.GT.K ) THEN
  287: *
  288: *                 C2 := C2 - V2 * W**H
  289: *
  290:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose',
  291:      $                        M-K, N, K, -ONE, V( K+1, 1 ), LDV, WORK,
  292:      $                        LDWORK, ONE, C( K+1, 1 ), LDC )
  293:                END IF
  294: *
  295: *              W := W * V1**H
  296: *
  297:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  298:      $                     'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
  299: *
  300: *              C1 := C1 - W**H
  301: *
  302:                DO 30 J = 1, K
  303:                   DO 20 I = 1, N
  304:                      C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
  305:    20             CONTINUE
  306:    30          CONTINUE
  307: *
  308:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  309: *
  310: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  311: *
  312: *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
  313: *
  314: *              W := C1
  315: *
  316:                DO 40 J = 1, K
  317:                   CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
  318:    40          CONTINUE
  319: *
  320: *              W := W * V1
  321: *
  322:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
  323:      $                     K, ONE, V, LDV, WORK, LDWORK )
  324:                IF( N.GT.K ) THEN
  325: *
  326: *                 W := W + C2 * V2
  327: *
  328:                   CALL ZGEMM( 'No transpose', 'No transpose', M, K, N-K,
  329:      $                        ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
  330:      $                        ONE, WORK, LDWORK )
  331:                END IF
  332: *
  333: *              W := W * T  or  W * T**H
  334: *
  335:                CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
  336:      $                     ONE, T, LDT, WORK, LDWORK )
  337: *
  338: *              C := C - W * V**H
  339: *
  340:                IF( N.GT.K ) THEN
  341: *
  342: *                 C2 := C2 - W * V2**H
  343: *
  344:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  345:      $                        N-K, K, -ONE, WORK, LDWORK, V( K+1, 1 ),
  346:      $                        LDV, ONE, C( 1, K+1 ), LDC )
  347:                END IF
  348: *
  349: *              W := W * V1**H
  350: *
  351:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  352:      $                     'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
  353: *
  354: *              C1 := C1 - W
  355: *
  356:                DO 60 J = 1, K
  357:                   DO 50 I = 1, M
  358:                      C( I, J ) = C( I, J ) - WORK( I, J )
  359:    50             CONTINUE
  360:    60          CONTINUE
  361:             END IF
  362: *
  363:          ELSE
  364: *
  365: *           Let  V =  ( V1 )
  366: *                     ( V2 )    (last K rows)
  367: *           where  V2  is unit upper triangular.
  368: *
  369:             IF( LSAME( SIDE, 'L' ) ) THEN
  370: *
  371: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  372: *                                                    ( C2 )
  373: *
  374: *              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
  375: *
  376: *              W := C2**H
  377: *
  378:                DO 70 J = 1, K
  379:                   CALL ZCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
  380:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  381:    70          CONTINUE
  382: *
  383: *              W := W * V2
  384: *
  385:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
  386:      $                     K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
  387:                IF( M.GT.K ) THEN
  388: *
  389: *                 W := W + C1**H * V1
  390: *
  391:                   CALL ZGEMM( 'Conjugate transpose', 'No transpose', N,
  392:      $                        K, M-K, ONE, C, LDC, V, LDV, ONE, WORK,
  393:      $                        LDWORK )
  394:                END IF
  395: *
  396: *              W := W * T**H  or  W * T
  397: *
  398:                CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
  399:      $                     ONE, T, LDT, WORK, LDWORK )
  400: *
  401: *              C := C - V * W**H
  402: *
  403:                IF( M.GT.K ) THEN
  404: *
  405: *                 C1 := C1 - V1 * W**H
  406: *
  407:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose',
  408:      $                        M-K, N, K, -ONE, V, LDV, WORK, LDWORK,
  409:      $                        ONE, C, LDC )
  410:                END IF
  411: *
  412: *              W := W * V2**H
  413: *
  414:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  415:      $                     'Unit', N, K, ONE, V( M-K+1, 1 ), LDV, WORK,
  416:      $                     LDWORK )
  417: *
  418: *              C2 := C2 - W**H
  419: *
  420:                DO 90 J = 1, K
  421:                   DO 80 I = 1, N
  422:                      C( M-K+J, I ) = C( M-K+J, I ) -
  423:      $                               DCONJG( WORK( I, J ) )
  424:    80             CONTINUE
  425:    90          CONTINUE
  426: *
  427:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  428: *
  429: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  430: *
  431: *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
  432: *
  433: *              W := C2
  434: *
  435:                DO 100 J = 1, K
  436:                   CALL ZCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
  437:   100          CONTINUE
  438: *
  439: *              W := W * V2
  440: *
  441:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
  442:      $                     K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
  443:                IF( N.GT.K ) THEN
  444: *
  445: *                 W := W + C1 * V1
  446: *
  447:                   CALL ZGEMM( 'No transpose', 'No transpose', M, K, N-K,
  448:      $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
  449:                END IF
  450: *
  451: *              W := W * T  or  W * T**H
  452: *
  453:                CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
  454:      $                     ONE, T, LDT, WORK, LDWORK )
  455: *
  456: *              C := C - W * V**H
  457: *
  458:                IF( N.GT.K ) THEN
  459: *
  460: *                 C1 := C1 - W * V1**H
  461: *
  462:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  463:      $                        N-K, K, -ONE, WORK, LDWORK, V, LDV, ONE,
  464:      $                        C, LDC )
  465:                END IF
  466: *
  467: *              W := W * V2**H
  468: *
  469:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  470:      $                     'Unit', M, K, ONE, V( N-K+1, 1 ), LDV, WORK,
  471:      $                     LDWORK )
  472: *
  473: *              C2 := C2 - W
  474: *
  475:                DO 120 J = 1, K
  476:                   DO 110 I = 1, M
  477:                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
  478:   110             CONTINUE
  479:   120          CONTINUE
  480:             END IF
  481:          END IF
  482: *
  483:       ELSE IF( LSAME( STOREV, 'R' ) ) THEN
  484: *
  485:          IF( LSAME( DIRECT, 'F' ) ) THEN
  486: *
  487: *           Let  V =  ( V1  V2 )    (V1: first K columns)
  488: *           where  V1  is unit upper triangular.
  489: *
  490:             IF( LSAME( SIDE, 'L' ) ) THEN
  491: *
  492: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  493: *                                                    ( C2 )
  494: *
  495: *              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
  496: *
  497: *              W := C1**H
  498: *
  499:                DO 130 J = 1, K
  500:                   CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
  501:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  502:   130          CONTINUE
  503: *
  504: *              W := W * V1**H
  505: *
  506:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  507:      $                     'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
  508:                IF( M.GT.K ) THEN
  509: *
  510: *                 W := W + C2**H * V2**H
  511: *
  512:                   CALL ZGEMM( 'Conjugate transpose',
  513:      $                        'Conjugate transpose', N, K, M-K, ONE,
  514:      $                        C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
  515:      $                        WORK, LDWORK )
  516:                END IF
  517: *
  518: *              W := W * T**H  or  W * T
  519: *
  520:                CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
  521:      $                     ONE, T, LDT, WORK, LDWORK )
  522: *
  523: *              C := C - V**H * W**H
  524: *
  525:                IF( M.GT.K ) THEN
  526: *
  527: *                 C2 := C2 - V2**H * W**H
  528: *
  529:                   CALL ZGEMM( 'Conjugate transpose',
  530:      $                        'Conjugate transpose', M-K, N, K, -ONE,
  531:      $                        V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
  532:      $                        C( K+1, 1 ), LDC )
  533:                END IF
  534: *
  535: *              W := W * V1
  536: *
  537:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
  538:      $                     K, ONE, V, LDV, WORK, LDWORK )
  539: *
  540: *              C1 := C1 - W**H
  541: *
  542:                DO 150 J = 1, K
  543:                   DO 140 I = 1, N
  544:                      C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
  545:   140             CONTINUE
  546:   150          CONTINUE
  547: *
  548:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  549: *
  550: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  551: *
  552: *              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
  553: *
  554: *              W := C1
  555: *
  556:                DO 160 J = 1, K
  557:                   CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
  558:   160          CONTINUE
  559: *
  560: *              W := W * V1**H
  561: *
  562:                CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
  563:      $                     'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
  564:                IF( N.GT.K ) THEN
  565: *
  566: *                 W := W + C2 * V2**H
  567: *
  568:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  569:      $                        K, N-K, ONE, C( 1, K+1 ), LDC,
  570:      $                        V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
  571:                END IF
  572: *
  573: *              W := W * T  or  W * T**H
  574: *
  575:                CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
  576:      $                     ONE, T, LDT, WORK, LDWORK )
  577: *
  578: *              C := C - W * V
  579: *
  580:                IF( N.GT.K ) THEN
  581: *
  582: *                 C2 := C2 - W * V2
  583: *
  584:                   CALL ZGEMM( 'No transpose', 'No transpose', M, N-K, K,
  585:      $                        -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
  586:      $                        C( 1, K+1 ), LDC )
  587:                END IF
  588: *
  589: *              W := W * V1
  590: *
  591:                CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
  592:      $                     K, ONE, V, LDV, WORK, LDWORK )
  593: *
  594: *              C1 := C1 - W
  595: *
  596:                DO 180 J = 1, K
  597:                   DO 170 I = 1, M
  598:                      C( I, J ) = C( I, J ) - WORK( I, J )
  599:   170             CONTINUE
  600:   180          CONTINUE
  601: *
  602:             END IF
  603: *
  604:          ELSE
  605: *
  606: *           Let  V =  ( V1  V2 )    (V2: last K columns)
  607: *           where  V2  is unit lower triangular.
  608: *
  609:             IF( LSAME( SIDE, 'L' ) ) THEN
  610: *
  611: *              Form  H * C  or  H**H * C  where  C = ( C1 )
  612: *                                                    ( C2 )
  613: *
  614: *              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
  615: *
  616: *              W := C2**H
  617: *
  618:                DO 190 J = 1, K
  619:                   CALL ZCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
  620:                   CALL ZLACGV( N, WORK( 1, J ), 1 )
  621:   190          CONTINUE
  622: *
  623: *              W := W * V2**H
  624: *
  625:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  626:      $                     'Unit', N, K, ONE, V( 1, M-K+1 ), LDV, WORK,
  627:      $                     LDWORK )
  628:                IF( M.GT.K ) THEN
  629: *
  630: *                 W := W + C1**H * V1**H
  631: *
  632:                   CALL ZGEMM( 'Conjugate transpose',
  633:      $                        'Conjugate transpose', N, K, M-K, ONE, C,
  634:      $                        LDC, V, LDV, ONE, WORK, LDWORK )
  635:                END IF
  636: *
  637: *              W := W * T**H  or  W * T
  638: *
  639:                CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
  640:      $                     ONE, T, LDT, WORK, LDWORK )
  641: *
  642: *              C := C - V**H * W**H
  643: *
  644:                IF( M.GT.K ) THEN
  645: *
  646: *                 C1 := C1 - V1**H * W**H
  647: *
  648:                   CALL ZGEMM( 'Conjugate transpose',
  649:      $                        'Conjugate transpose', M-K, N, K, -ONE, V,
  650:      $                        LDV, WORK, LDWORK, ONE, C, LDC )
  651:                END IF
  652: *
  653: *              W := W * V2
  654: *
  655:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
  656:      $                     K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
  657: *
  658: *              C2 := C2 - W**H
  659: *
  660:                DO 210 J = 1, K
  661:                   DO 200 I = 1, N
  662:                      C( M-K+J, I ) = C( M-K+J, I ) -
  663:      $                               DCONJG( WORK( I, J ) )
  664:   200             CONTINUE
  665:   210          CONTINUE
  666: *
  667:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  668: *
  669: *              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
  670: *
  671: *              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
  672: *
  673: *              W := C2
  674: *
  675:                DO 220 J = 1, K
  676:                   CALL ZCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
  677:   220          CONTINUE
  678: *
  679: *              W := W * V2**H
  680: *
  681:                CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
  682:      $                     'Unit', M, K, ONE, V( 1, N-K+1 ), LDV, WORK,
  683:      $                     LDWORK )
  684:                IF( N.GT.K ) THEN
  685: *
  686: *                 W := W + C1 * V1**H
  687: *
  688:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
  689:      $                        K, N-K, ONE, C, LDC, V, LDV, ONE, WORK,
  690:      $                        LDWORK )
  691:                END IF
  692: *
  693: *              W := W * T  or  W * T**H
  694: *
  695:                CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
  696:      $                     ONE, T, LDT, WORK, LDWORK )
  697: *
  698: *              C := C - W * V
  699: *
  700:                IF( N.GT.K ) THEN
  701: *
  702: *                 C1 := C1 - W * V1
  703: *
  704:                   CALL ZGEMM( 'No transpose', 'No transpose', M, N-K, K,
  705:      $                        -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
  706:                END IF
  707: *
  708: *              W := W * V2
  709: *
  710:                CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
  711:      $                     K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
  712: *
  713: *              C1 := C1 - W
  714: *
  715:                DO 240 J = 1, K
  716:                   DO 230 I = 1, M
  717:                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
  718:   230             CONTINUE
  719:   240          CONTINUE
  720: *
  721:             END IF
  722: *
  723:          END IF
  724:       END IF
  725: *
  726:       RETURN
  727: *
  728: *     End of ZLARFB
  729: *
  730:       END

CVSweb interface <joel.bertrand@systella.fr>