[BACK]Return to dlarfb.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dlarfb.f between versions 1.2 and 1.12

version 1.2, 2010/04/21 13:45:18 version 1.12, 2012/12/14 12:30:24
Line 1 Line 1
   *> \brief \b DLARFB applies a block reflector or its transpose to a general rectangular matrix.
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at 
   *            http://www.netlib.org/lapack/explore-html/ 
   *
   *> \htmlonly
   *> Download DLARFB + dependencies 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb.f"> 
   *> [TGZ]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb.f"> 
   *> [ZIP]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb.f"> 
   *> [TXT]</a>
   *> \endhtmlonly 
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
   *                          T, LDT, C, LDC, WORK, LDWORK )
   * 
   *       .. Scalar Arguments ..
   *       CHARACTER          DIRECT, SIDE, STOREV, TRANS
   *       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
   *      $                   WORK( LDWORK, * )
   *       ..
   *  
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DLARFB applies a real block reflector H or its transpose H**T to a
   *> real m by n matrix C, from either the left or the right.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] SIDE
   *> \verbatim
   *>          SIDE is CHARACTER*1
   *>          = 'L': apply H or H**T from the Left
   *>          = 'R': apply H or H**T from the Right
   *> \endverbatim
   *>
   *> \param[in] TRANS
   *> \verbatim
   *>          TRANS is CHARACTER*1
   *>          = 'N': apply H (No transpose)
   *>          = 'T': apply H**T (Transpose)
   *> \endverbatim
   *>
   *> \param[in] DIRECT
   *> \verbatim
   *>          DIRECT is CHARACTER*1
   *>          Indicates how H is formed from a product of elementary
   *>          reflectors
   *>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
   *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
   *> \endverbatim
   *>
   *> \param[in] STOREV
   *> \verbatim
   *>          STOREV is CHARACTER*1
   *>          Indicates how the vectors which define the elementary
   *>          reflectors are stored:
   *>          = 'C': Columnwise
   *>          = 'R': Rowwise
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The number of rows of the matrix C.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The number of columns of the matrix C.
   *> \endverbatim
   *>
   *> \param[in] K
   *> \verbatim
   *>          K is INTEGER
   *>          The order of the matrix T (= the number of elementary
   *>          reflectors whose product defines the block reflector).
   *> \endverbatim
   *>
   *> \param[in] V
   *> \verbatim
   *>          V is DOUBLE PRECISION array, dimension
   *>                                (LDV,K) if STOREV = 'C'
   *>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
   *>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
   *>          The matrix V. See Further Details.
   *> \endverbatim
   *>
   *> \param[in] LDV
   *> \verbatim
   *>          LDV is INTEGER
   *>          The leading dimension of the array V.
   *>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
   *>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
   *>          if STOREV = 'R', LDV >= K.
   *> \endverbatim
   *>
   *> \param[in] T
   *> \verbatim
   *>          T is DOUBLE PRECISION array, dimension (LDT,K)
   *>          The triangular k by k matrix T in the representation of the
   *>          block reflector.
   *> \endverbatim
   *>
   *> \param[in] LDT
   *> \verbatim
   *>          LDT is INTEGER
   *>          The leading dimension of the array T. LDT >= K.
   *> \endverbatim
   *>
   *> \param[in,out] C
   *> \verbatim
   *>          C is DOUBLE PRECISION array, dimension (LDC,N)
   *>          On entry, the m by n matrix C.
   *>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
   *> \endverbatim
   *>
   *> \param[in] LDC
   *> \verbatim
   *>          LDC is INTEGER
   *>          The leading dimension of the array C. LDC >= max(1,M).
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)
   *> \endverbatim
   *>
   *> \param[in] LDWORK
   *> \verbatim
   *>          LDWORK is INTEGER
   *>          The leading dimension of the array WORK.
   *>          If SIDE = 'L', LDWORK >= max(1,N);
   *>          if SIDE = 'R', LDWORK >= max(1,M).
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee 
   *> \author Univ. of California Berkeley 
   *> \author Univ. of Colorado Denver 
   *> \author NAG Ltd. 
   *
   *> \date September 2012
   *
   *> \ingroup doubleOTHERauxiliary
   *
   *> \par Further Details:
   *  =====================
   *>
   *> \verbatim
   *>
   *>  The shape of the matrix V and the storage of the vectors which define
   *>  the H(i) is best illustrated by the following example with n = 5 and
   *>  k = 3. The elements equal to 1 are not stored; the corresponding
   *>  array elements are modified but restored on exit. The rest of the
   *>  array is not used.
   *>
   *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
   *>
   *>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
   *>                   ( v1  1    )                     (     1 v2 v2 v2 )
   *>                   ( v1 v2  1 )                     (        1 v3 v3 )
   *>                   ( v1 v2 v3 )
   *>                   ( v1 v2 v3 )
   *>
   *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
   *>
   *>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
   *>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
   *>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
   *>                   (     1 v3 )
   *>                   (        1 )
   *> \endverbatim
   *>
   *  =====================================================================
       SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,        SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
      $                   T, LDT, C, LDC, WORK, LDWORK )       $                   T, LDT, C, LDC, WORK, LDWORK )
       IMPLICIT NONE  
 *  *
 *  -- LAPACK auxiliary routine (version 3.2) --  *  -- LAPACK auxiliary routine (version 3.4.2) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2006  *     September 2012
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          DIRECT, SIDE, STOREV, TRANS        CHARACTER          DIRECT, SIDE, STOREV, TRANS
Line 16 Line 209
      $                   WORK( LDWORK, * )       $                   WORK( LDWORK, * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DLARFB applies a real block reflector H or its transpose H' to a  
 *  real m by n matrix C, from either the left or the right.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  SIDE    (input) CHARACTER*1  
 *          = 'L': apply H or H' from the Left  
 *          = 'R': apply H or H' from the Right  
 *  
 *  TRANS   (input) CHARACTER*1  
 *          = 'N': apply H (No transpose)  
 *          = 'T': apply H' (Transpose)  
 *  
 *  DIRECT  (input) CHARACTER*1  
 *          Indicates how H is formed from a product of elementary  
 *          reflectors  
 *          = 'F': H = H(1) H(2) . . . H(k) (Forward)  
 *          = 'B': H = H(k) . . . H(2) H(1) (Backward)  
 *  
 *  STOREV  (input) CHARACTER*1  
 *          Indicates how the vectors which define the elementary  
 *          reflectors are stored:  
 *          = 'C': Columnwise  
 *          = 'R': Rowwise  
 *  
 *  M       (input) INTEGER  
 *          The number of rows of the matrix C.  
 *  
 *  N       (input) INTEGER  
 *          The number of columns of the matrix C.  
 *  
 *  K       (input) INTEGER  
 *          The order of the matrix T (= the number of elementary  
 *          reflectors whose product defines the block reflector).  
 *  
 *  V       (input) DOUBLE PRECISION array, dimension  
 *                                (LDV,K) if STOREV = 'C'  
 *                                (LDV,M) if STOREV = 'R' and SIDE = 'L'  
 *                                (LDV,N) if STOREV = 'R' and SIDE = 'R'  
 *          The matrix V. See further details.  
 *  
 *  LDV     (input) INTEGER  
 *          The leading dimension of the array V.  
 *          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);  
 *          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);  
 *          if STOREV = 'R', LDV >= K.  
 *  
 *  T       (input) DOUBLE PRECISION array, dimension (LDT,K)  
 *          The triangular k by k matrix T in the representation of the  
 *          block reflector.  
 *  
 *  LDT     (input) INTEGER  
 *          The leading dimension of the array T. LDT >= K.  
 *  
 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)  
 *          On entry, the m by n matrix C.  
 *          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.  
 *  
 *  LDC     (input) INTEGER  
 *          The leading dimension of the array C. LDA >= max(1,M).  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)  
 *  
 *  LDWORK  (input) INTEGER  
 *          The leading dimension of the array WORK.  
 *          If SIDE = 'L', LDWORK >= max(1,N);  
 *          if SIDE = 'R', LDWORK >= max(1,M).  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 96 Line 217
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       CHARACTER          TRANST        CHARACTER          TRANST
       INTEGER            I, J, LASTV, LASTC        INTEGER            I, J, LASTV, LASTC, lastv2
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
Line 129 Line 250
 *  *
             IF( LSAME( SIDE, 'L' ) ) THEN              IF( LSAME( SIDE, 'L' ) ) THEN
 *  *
 *              Form  H * C  or  H' * C  where  C = ( C1 )  *              Form  H * C  or  H**T * C  where  C = ( C1 )
 *                                                  ( C2 )  *                                                    ( C2 )
 *  *
                LASTV = MAX( K, ILADLR( M, K, V, LDV ) )                 LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
                LASTC = ILADLC( LASTV, N, C, LDC )                 LASTC = ILADLC( LASTV, N, C, LDC )
 *  *
 *              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)  *              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
 *  *
 *              W := C1'  *              W := C1**T
 *  *
                DO 10 J = 1, K                 DO 10 J = 1, K
                   CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )                    CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
Line 149 Line 270
      $              LASTC, K, ONE, V, LDV, WORK, LDWORK )       $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( LASTV.GT.K ) THEN
 *  *
 *                 W := W + C2'*V2  *                 W := W + C2**T *V2
 *  *
                   CALL DGEMM( 'Transpose', 'No transpose',                    CALL DGEMM( 'Transpose', 'No transpose',
      $                 LASTC, K, LASTV-K,       $                 LASTC, K, LASTV-K,
Line 157 Line 278
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T'  or  W * T  *              W := W * T**T  or  W * T
 *  *
                CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',                 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - V * W'  *              C := C - V * W**T
 *  *
                IF( LASTV.GT.K ) THEN                 IF( LASTV.GT.K ) THEN
 *  *
 *                 C2 := C2 - V2 * W'  *                 C2 := C2 - V2 * W**T
 *  *
                   CALL DGEMM( 'No transpose', 'Transpose',                    CALL DGEMM( 'No transpose', 'Transpose',
      $                 LASTV-K, LASTC, K,       $                 LASTV-K, LASTC, K,
Line 174 Line 295
      $                 C( K+1, 1 ), LDC )       $                 C( K+1, 1 ), LDC )
                END IF                 END IF
 *  *
 *              W := W * V1'  *              W := W * V1**T
 *  *
                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V, LDV, WORK, LDWORK )       $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 *  *
 *              C1 := C1 - W'  *              C1 := C1 - W**T
 *  *
                DO 30 J = 1, K                 DO 30 J = 1, K
                   DO 20 I = 1, LASTC                    DO 20 I = 1, LASTC
Line 189 Line 310
 *  *
             ELSE IF( LSAME( SIDE, 'R' ) ) THEN              ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 *  *
 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )  *              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 *  *
                LASTV = MAX( K, ILADLR( N, K, V, LDV ) )                 LASTV = MAX( K, ILADLR( N, K, V, LDV ) )
                LASTC = ILADLR( M, LASTV, C, LDC )                 LASTC = ILADLR( M, LASTV, C, LDC )
Line 216 Line 337
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T  or  W * T'  *              W := W * T  or  W * T**T
 *  *
                CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',                 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - W * V'  *              C := C - W * V**T
 *  *
                IF( LASTV.GT.K ) THEN                 IF( LASTV.GT.K ) THEN
 *  *
 *                 C2 := C2 - W * V2'  *                 C2 := C2 - W * V2**T
 *  *
                   CALL DGEMM( 'No transpose', 'Transpose',                    CALL DGEMM( 'No transpose', 'Transpose',
      $                 LASTC, LASTV-K, K,       $                 LASTC, LASTV-K, K,
Line 233 Line 354
      $                 C( 1, K+1 ), LDC )       $                 C( 1, K+1 ), LDC )
                END IF                 END IF
 *  *
 *              W := W * V1'  *              W := W * V1**T
 *  *
                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V, LDV, WORK, LDWORK )       $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
Line 255 Line 376
 *  *
             IF( LSAME( SIDE, 'L' ) ) THEN              IF( LSAME( SIDE, 'L' ) ) THEN
 *  *
 *              Form  H * C  or  H' * C  where  C = ( C1 )  *              Form  H * C  or  H**T * C  where  C = ( C1 )
 *                                                  ( C2 )  *                                                    ( C2 )
 *  *
                LASTV = MAX( K, ILADLR( M, K, V, LDV ) )                 LASTC = ILADLC( M, N, C, LDC )
                LASTC = ILADLC( LASTV, N, C, LDC )  
 *  *
 *              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)  *              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
 *  *
 *              W := C2'  *              W := C2**T
 *  *
                DO 70 J = 1, K                 DO 70 J = 1, K
                   CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,                    CALL DCOPY( LASTC, C( M-K+J, 1 ), LDC,
      $                 WORK( 1, J ), 1 )       $                 WORK( 1, J ), 1 )
    70          CONTINUE     70          CONTINUE
 *  *
 *              W := W * V2  *              W := W * V2
 *  *
                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
      $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,       $              LASTC, K, ONE, V( M-K+1, 1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( M.GT.K ) THEN
 *  *
 *                 W := W + C1'*V1  *                 W := W + C1**T*V1
 *  *
                   CALL DGEMM( 'Transpose', 'No transpose',                    CALL DGEMM( 'Transpose', 'No transpose',
      $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,       $                 LASTC, K, M-K, ONE, C, LDC, V, LDV,
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T'  or  W * T  *              W := W * T**T  or  W * T
 *  *
                CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',                 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - V * W'  *              C := C - V * W**T
 *  *
                IF( LASTV.GT.K ) THEN                 IF( M.GT.K ) THEN
 *  *
 *                 C1 := C1 - V1 * W'  *                 C1 := C1 - V1 * W**T
 *  *
                   CALL DGEMM( 'No transpose', 'Transpose',                    CALL DGEMM( 'No transpose', 'Transpose',
      $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,       $                 M-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
      $                 ONE, C, LDC )       $                 ONE, C, LDC )
                END IF                 END IF
 *  *
 *              W := W * V2'  *              W := W * V2**T
 *  *
                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,       $              LASTC, K, ONE, V( M-K+1, 1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
 *  *
 *              C2 := C2 - W'  *              C2 := C2 - W**T
 *  *
                DO 90 J = 1, K                 DO 90 J = 1, K
                   DO 80 I = 1, LASTC                    DO 80 I = 1, LASTC
                      C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)                       C( M-K+J, I ) = C( M-K+J, I ) - WORK(I, J)
    80             CONTINUE     80             CONTINUE
    90          CONTINUE     90          CONTINUE
 *  *
             ELSE IF( LSAME( SIDE, 'R' ) ) THEN              ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 *  *
 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )  *              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 *  *
                LASTV = MAX( K, ILADLR( N, K, V, LDV ) )                 LASTC = ILADLR( M, N, C, LDC )
                LASTC = ILADLR( M, LASTV, C, LDC )  
 *  *
 *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)  *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
 *  *
Line 332 Line 451
 *              W := W * V2  *              W := W * V2
 *  *
                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
      $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,       $              LASTC, K, ONE, V( N-K+1, 1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( N.GT.K ) THEN
 *  *
 *                 W := W + C1 * V1  *                 W := W + C1 * V1
 *  *
                   CALL DGEMM( 'No transpose', 'No transpose',                    CALL DGEMM( 'No transpose', 'No transpose',
      $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,       $                 LASTC, K, N-K, ONE, C, LDC, V, LDV,
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T  or  W * T'  *              W := W * T  or  W * T**T
 *  *
                CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',                 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - W * V'  *              C := C - W * V**T
 *  *
                IF( LASTV.GT.K ) THEN                 IF( N.GT.K ) THEN
 *  *
 *                 C1 := C1 - W * V1'  *                 C1 := C1 - W * V1**T
 *  *
                   CALL DGEMM( 'No transpose', 'Transpose',                    CALL DGEMM( 'No transpose', 'Transpose',
      $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,       $                 LASTC, N-K, K, -ONE, WORK, LDWORK, V, LDV,
      $                 ONE, C, LDC )       $                 ONE, C, LDC )
                END IF                 END IF
 *  *
 *              W := W * V2'  *              W := W * V2**T
 *  *
                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,       $              LASTC, K, ONE, V( N-K+1, 1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
 *  *
 *              C2 := C2 - W  *              C2 := C2 - W
 *  *
                DO 120 J = 1, K                 DO 120 J = 1, K
                   DO 110 I = 1, LASTC                    DO 110 I = 1, LASTC
                      C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)                       C( I, N-K+J ) = C( I, N-K+J ) - WORK(I, J)
   110             CONTINUE    110             CONTINUE
   120          CONTINUE    120          CONTINUE
             END IF              END IF
Line 384 Line 503
 *  *
             IF( LSAME( SIDE, 'L' ) ) THEN              IF( LSAME( SIDE, 'L' ) ) THEN
 *  *
 *              Form  H * C  or  H' * C  where  C = ( C1 )  *              Form  H * C  or  H**T * C  where  C = ( C1 )
 *                                                  ( C2 )  *                                                    ( C2 )
 *  *
                LASTV = MAX( K, ILADLC( K, M, V, LDV ) )                 LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
                LASTC = ILADLC( LASTV, N, C, LDC )                 LASTC = ILADLC( LASTV, N, C, LDC )
 *  *
 *              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)  *              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
 *  *
 *              W := C1'  *              W := C1**T
 *  *
                DO 130 J = 1, K                 DO 130 J = 1, K
                   CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )                    CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
   130          CONTINUE    130          CONTINUE
 *  *
 *              W := W * V1'  *              W := W * V1**T
 *  *
                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V, LDV, WORK, LDWORK )       $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( LASTV.GT.K ) THEN
 *  *
 *                 W := W + C2'*V2'  *                 W := W + C2**T*V2**T
 *  *
                   CALL DGEMM( 'Transpose', 'Transpose',                    CALL DGEMM( 'Transpose', 'Transpose',
      $                 LASTC, K, LASTV-K,       $                 LASTC, K, LASTV-K,
Line 412 Line 531
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T'  or  W * T  *              W := W * T**T  or  W * T
 *  *
                CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',                 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - V' * W'  *              C := C - V**T * W**T
 *  *
                IF( LASTV.GT.K ) THEN                 IF( LASTV.GT.K ) THEN
 *  *
 *                 C2 := C2 - V2' * W'  *                 C2 := C2 - V2**T * W**T
 *  *
                   CALL DGEMM( 'Transpose', 'Transpose',                    CALL DGEMM( 'Transpose', 'Transpose',
      $                 LASTV-K, LASTC, K,       $                 LASTV-K, LASTC, K,
Line 434 Line 553
                CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
      $              LASTC, K, ONE, V, LDV, WORK, LDWORK )       $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 *  *
 *              C1 := C1 - W'  *              C1 := C1 - W**T
 *  *
                DO 150 J = 1, K                 DO 150 J = 1, K
                   DO 140 I = 1, LASTC                    DO 140 I = 1, LASTC
Line 444 Line 563
 *  *
             ELSE IF( LSAME( SIDE, 'R' ) ) THEN              ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 *  *
 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )  *              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 *  *
                LASTV = MAX( K, ILADLC( K, N, V, LDV ) )                 LASTV = MAX( K, ILADLC( K, N, V, LDV ) )
                LASTC = ILADLR( M, LASTV, C, LDC )                 LASTC = ILADLR( M, LASTV, C, LDC )
 *  *
 *              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)  *              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
 *  *
 *              W := C1  *              W := C1
 *  *
Line 457 Line 576
                   CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )                    CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
   160          CONTINUE    160          CONTINUE
 *  *
 *              W := W * V1'  *              W := W * V1**T
 *  *
                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V, LDV, WORK, LDWORK )       $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( LASTV.GT.K ) THEN
 *  *
 *                 W := W + C2 * V2'  *                 W := W + C2 * V2**T
 *  *
                   CALL DGEMM( 'No transpose', 'Transpose',                    CALL DGEMM( 'No transpose', 'Transpose',
      $                 LASTC, K, LASTV-K,       $                 LASTC, K, LASTV-K,
Line 471 Line 590
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T  or  W * T'  *              W := W * T  or  W * T**T
 *  *
                CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',                 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
Line 510 Line 629
 *  *
             IF( LSAME( SIDE, 'L' ) ) THEN              IF( LSAME( SIDE, 'L' ) ) THEN
 *  *
 *              Form  H * C  or  H' * C  where  C = ( C1 )  *              Form  H * C  or  H**T * C  where  C = ( C1 )
 *                                                  ( C2 )  *                                                    ( C2 )
 *  *
                LASTV = MAX( K, ILADLC( K, M, V, LDV ) )                 LASTC = ILADLC( M, N, C, LDC )
                LASTC = ILADLC( LASTV, N, C, LDC )  
 *  *
 *              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)  *              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
 *  *
 *              W := C2'  *              W := C2**T
 *  *
                DO 190 J = 1, K                 DO 190 J = 1, K
                   CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,                    CALL DCOPY( LASTC, C( M-K+J, 1 ), LDC,
      $                 WORK( 1, J ), 1 )       $                 WORK( 1, J ), 1 )
   190          CONTINUE    190          CONTINUE
 *  *
 *              W := W * V2'  *              W := W * V2**T
 *  *
                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,       $              LASTC, K, ONE, V( 1, M-K+1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( M.GT.K ) THEN
 *  *
 *                 W := W + C1'*V1'  *                 W := W + C1**T * V1**T
 *  *
                   CALL DGEMM( 'Transpose', 'Transpose',                    CALL DGEMM( 'Transpose', 'Transpose',
      $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,       $                 LASTC, K, M-K, ONE, C, LDC, V, LDV,
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T'  or  W * T  *              W := W * T**T  or  W * T
 *  *
                CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',                 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - V' * W'  *              C := C - V**T * W**T
 *  *
                IF( LASTV.GT.K ) THEN                 IF( M.GT.K ) THEN
 *  *
 *                 C1 := C1 - V1' * W'  *                 C1 := C1 - V1**T * W**T
 *  *
                   CALL DGEMM( 'Transpose', 'Transpose',                    CALL DGEMM( 'Transpose', 'Transpose',
      $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,       $                 M-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
      $                 ONE, C, LDC )       $                 ONE, C, LDC )
                END IF                 END IF
 *  *
 *              W := W * V2  *              W := W * V2
 *  *
                CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',                 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
      $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,       $              LASTC, K, ONE, V( 1, M-K+1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
 *  *
 *              C2 := C2 - W'  *              C2 := C2 - W**T
 *  *
                DO 210 J = 1, K                 DO 210 J = 1, K
                   DO 200 I = 1, LASTC                    DO 200 I = 1, LASTC
                      C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)                       C( M-K+J, I ) = C( M-K+J, I ) - WORK(I, J)
   200             CONTINUE    200             CONTINUE
   210          CONTINUE    210          CONTINUE
 *  *
             ELSE IF( LSAME( SIDE, 'R' ) ) THEN              ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 *  *
 *              Form  C * H  or  C * H'  where  C = ( C1  C2 )  *              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 *  *
                LASTV = MAX( K, ILADLC( K, N, V, LDV ) )                 LASTC = ILADLR( M, N, C, LDC )
                LASTC = ILADLR( M, LASTV, C, LDC )  
 *  *
 *              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)  *              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
 *  *
 *              W := C2  *              W := C2
 *  *
                DO 220 J = 1, K                 DO 220 J = 1, K
                   CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1,                    CALL DCOPY( LASTC, C( 1, N-K+J ), 1,
      $                 WORK( 1, J ), 1 )       $                 WORK( 1, J ), 1 )
   220          CONTINUE    220          CONTINUE
 *  *
 *              W := W * V2'  *              W := W * V2**T
 *  *
                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',                 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
      $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,       $              LASTC, K, ONE, V( 1, N-K+1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
                IF( LASTV.GT.K ) THEN                 IF( N.GT.K ) THEN
 *  *
 *                 W := W + C1 * V1'  *                 W := W + C1 * V1**T
 *  *
                   CALL DGEMM( 'No transpose', 'Transpose',                    CALL DGEMM( 'No transpose', 'Transpose',
      $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,       $                 LASTC, K, N-K, ONE, C, LDC, V, LDV,
      $                 ONE, WORK, LDWORK )       $                 ONE, WORK, LDWORK )
                END IF                 END IF
 *  *
 *              W := W * T  or  W * T'  *              W := W * T  or  W * T**T
 *  *
                CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',                 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
      $              LASTC, K, ONE, T, LDT, WORK, LDWORK )       $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 *  *
 *              C := C - W * V  *              C := C - W * V
 *  *
                IF( LASTV.GT.K ) THEN                 IF( N.GT.K ) THEN
 *  *
 *                 C1 := C1 - W * V1  *                 C1 := C1 - W * V1
 *  *
                   CALL DGEMM( 'No transpose', 'No transpose',                    CALL DGEMM( 'No transpose', 'No transpose',
      $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,       $                 LASTC, N-K, K, -ONE, WORK, LDWORK, V, LDV,
      $                 ONE, C, LDC )       $                 ONE, C, LDC )
                END IF                 END IF
 *  *
 *              W := W * V2  *              W := W * V2
 *  *
                CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',                 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
      $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,       $              LASTC, K, ONE, V( 1, N-K+1 ), LDV,
      $              WORK, LDWORK )       $              WORK, LDWORK )
 *  *
 *              C1 := C1 - W  *              C1 := C1 - W
 *  *
                DO 240 J = 1, K                 DO 240 J = 1, K
                   DO 230 I = 1, LASTC                    DO 230 I = 1, LASTC
                      C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)                       C( I, N-K+J ) = C( I, N-K+J ) - WORK(I, J)
   230             CONTINUE    230             CONTINUE
   240          CONTINUE    240          CONTINUE
 *  *
Removed from v.1.2  
changed lines
  Added in v.1.12

CVSweb interface <joel.bertrand@systella.fr>