File:  [local] / rpl / lapack / lapack / dlarfb.f
Revision 1.8: download - view: text, annotated - select for diffs - revision graph
Fri Jul 22 07:38:07 2011 UTC (12 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_3, rpl-4_1_2, rpl-4_1_1, HEAD
En route vers la 4.4.1.

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

CVSweb interface <joel.bertrand@systella.fr>