File:  [local] / rpl / lapack / lapack / zlarfb.f
Revision 1.20: download - view: text, annotated - select for diffs - revision graph
Thu May 21 21:46:09 2020 UTC (4 years ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_33, rpl-4_1_32, HEAD
Mise à jour de Lapack.

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

CVSweb interface <joel.bertrand@systella.fr>