--- rpl/lapack/lapack/dlarfb.f 2010/08/06 15:28:41 1.3
+++ rpl/lapack/lapack/dlarfb.f 2011/11/21 20:42:57 1.9
@@ -1,11 +1,204 @@
+*> \brief \b DLARFB
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DLARFB + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \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 November 2011
+*
+*> \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,
$ T, LDT, C, LDC, WORK, LDWORK )
- IMPLICIT NONE
*
-* -- LAPACK auxiliary routine (version 3.2) --
+* -- LAPACK auxiliary routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
+* November 2011
*
* .. Scalar Arguments ..
CHARACTER DIRECT, SIDE, STOREV, TRANS
@@ -16,78 +209,6 @@
$ 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 ..
@@ -129,15 +250,15 @@
*
IF( LSAME( SIDE, 'L' ) ) THEN
*
-* Form H * C or H' * C where C = ( C1 )
-* ( C2 )
+* Form H * C or H**T * C where C = ( C1 )
+* ( C2 )
*
LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
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
CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
@@ -149,7 +270,7 @@
$ LASTC, K, ONE, V, LDV, WORK, LDWORK )
IF( LASTV.GT.K ) THEN
*
-* W := W + C2'*V2
+* W := W + C2**T *V2
*
CALL DGEMM( 'Transpose', 'No transpose',
$ LASTC, K, LASTV-K,
@@ -157,16 +278,16 @@
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T' or W * T
+* W := W * T**T or W * T
*
CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
*
-* C := C - V * W'
+* C := C - V * W**T
*
IF( LASTV.GT.K ) THEN
*
-* C2 := C2 - V2 * W'
+* C2 := C2 - V2 * W**T
*
CALL DGEMM( 'No transpose', 'Transpose',
$ LASTV-K, LASTC, K,
@@ -174,12 +295,12 @@
$ C( K+1, 1 ), LDC )
END IF
*
-* W := W * V1'
+* W := W * V1**T
*
CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
$ LASTC, K, ONE, V, LDV, WORK, LDWORK )
*
-* C1 := C1 - W'
+* C1 := C1 - W**T
*
DO 30 J = 1, K
DO 20 I = 1, LASTC
@@ -189,7 +310,7 @@
*
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, LASTV, C, LDC )
@@ -216,16 +337,16 @@
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T or W * T'
+* W := W * T or W * T**T
*
CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
*
-* C := C - W * V'
+* C := C - W * V**T
*
IF( LASTV.GT.K ) THEN
*
-* C2 := C2 - W * V2'
+* C2 := C2 - W * V2**T
*
CALL DGEMM( 'No transpose', 'Transpose',
$ LASTC, LASTV-K, K,
@@ -233,7 +354,7 @@
$ C( 1, K+1 ), LDC )
END IF
*
-* W := W * V1'
+* W := W * V1**T
*
CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
$ LASTC, K, ONE, V, LDV, WORK, LDWORK )
@@ -255,15 +376,15 @@
*
IF( LSAME( SIDE, 'L' ) ) THEN
*
-* Form H * C or H' * C where C = ( C1 )
-* ( C2 )
+* Form H * C or H**T * C where C = ( C1 )
+* ( C2 )
*
LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
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
CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
@@ -277,36 +398,36 @@
$ WORK, LDWORK )
IF( LASTV.GT.K ) THEN
*
-* W := W + C1'*V1
+* W := W + C1**T*V1
*
CALL DGEMM( 'Transpose', 'No transpose',
$ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T' or W * T
+* W := W * T**T or W * T
*
CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
*
-* C := C - V * W'
+* C := C - V * W**T
*
IF( LASTV.GT.K ) THEN
*
-* C1 := C1 - V1 * W'
+* C1 := C1 - V1 * W**T
*
CALL DGEMM( 'No transpose', 'Transpose',
$ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
$ ONE, C, LDC )
END IF
*
-* W := W * V2'
+* W := W * V2**T
*
CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
$ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
$ WORK, LDWORK )
*
-* C2 := C2 - W'
+* C2 := C2 - W**T
*
DO 90 J = 1, K
DO 80 I = 1, LASTC
@@ -316,7 +437,7 @@
*
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, LASTV, C, LDC )
@@ -343,23 +464,23 @@
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T or W * T'
+* W := W * T or W * T**T
*
CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
*
-* C := C - W * V'
+* C := C - W * V**T
*
IF( LASTV.GT.K ) THEN
*
-* C1 := C1 - W * V1'
+* C1 := C1 - W * V1**T
*
CALL DGEMM( 'No transpose', 'Transpose',
$ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
$ ONE, C, LDC )
END IF
*
-* W := W * V2'
+* W := W * V2**T
*
CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
$ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
@@ -384,27 +505,27 @@
*
IF( LSAME( SIDE, 'L' ) ) THEN
*
-* Form H * C or H' * C where C = ( C1 )
-* ( C2 )
+* Form H * C or H**T * C where C = ( C1 )
+* ( C2 )
*
LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
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
CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
130 CONTINUE
*
-* W := W * V1'
+* W := W * V1**T
*
CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
$ LASTC, K, ONE, V, LDV, WORK, LDWORK )
IF( LASTV.GT.K ) THEN
*
-* W := W + C2'*V2'
+* W := W + C2**T*V2**T
*
CALL DGEMM( 'Transpose', 'Transpose',
$ LASTC, K, LASTV-K,
@@ -412,16 +533,16 @@
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T' or W * T
+* W := W * T**T or W * T
*
CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
*
-* C := C - V' * W'
+* C := C - V**T * W**T
*
IF( LASTV.GT.K ) THEN
*
-* C2 := C2 - V2' * W'
+* C2 := C2 - V2**T * W**T
*
CALL DGEMM( 'Transpose', 'Transpose',
$ LASTV-K, LASTC, K,
@@ -434,7 +555,7 @@
CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
$ LASTC, K, ONE, V, LDV, WORK, LDWORK )
*
-* C1 := C1 - W'
+* C1 := C1 - W**T
*
DO 150 J = 1, K
DO 140 I = 1, LASTC
@@ -444,12 +565,12 @@
*
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, 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
*
@@ -457,13 +578,13 @@
CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
160 CONTINUE
*
-* W := W * V1'
+* W := W * V1**T
*
CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
$ LASTC, K, ONE, V, LDV, WORK, LDWORK )
IF( LASTV.GT.K ) THEN
*
-* W := W + C2 * V2'
+* W := W + C2 * V2**T
*
CALL DGEMM( 'No transpose', 'Transpose',
$ LASTC, K, LASTV-K,
@@ -471,7 +592,7 @@
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T or W * T'
+* W := W * T or W * T**T
*
CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
@@ -510,45 +631,45 @@
*
IF( LSAME( SIDE, 'L' ) ) THEN
*
-* Form H * C or H' * C where C = ( C1 )
-* ( C2 )
+* Form H * C or H**T * C where C = ( C1 )
+* ( C2 )
*
LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
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
CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
$ WORK( 1, J ), 1 )
190 CONTINUE
*
-* W := W * V2'
+* W := W * V2**T
*
CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
$ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
$ WORK, LDWORK )
IF( LASTV.GT.K ) THEN
*
-* W := W + C1'*V1'
+* W := W + C1**T * V1**T
*
CALL DGEMM( 'Transpose', 'Transpose',
$ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T' or W * T
+* W := W * T**T or W * T
*
CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )
*
-* C := C - V' * W'
+* C := C - V**T * W**T
*
IF( LASTV.GT.K ) THEN
*
-* C1 := C1 - V1' * W'
+* C1 := C1 - V1**T * W**T
*
CALL DGEMM( 'Transpose', 'Transpose',
$ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
@@ -561,7 +682,7 @@
$ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
$ WORK, LDWORK )
*
-* C2 := C2 - W'
+* C2 := C2 - W**T
*
DO 210 J = 1, K
DO 200 I = 1, LASTC
@@ -571,12 +692,12 @@
*
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, 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
*
@@ -585,21 +706,21 @@
$ WORK( 1, J ), 1 )
220 CONTINUE
*
-* W := W * V2'
+* W := W * V2**T
*
CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
$ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
$ WORK, LDWORK )
IF( LASTV.GT.K ) THEN
*
-* W := W + C1 * V1'
+* W := W + C1 * V1**T
*
CALL DGEMM( 'No transpose', 'Transpose',
$ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
$ ONE, WORK, LDWORK )
END IF
*
-* W := W * T or W * T'
+* W := W * T or W * T**T
*
CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
$ LASTC, K, ONE, T, LDT, WORK, LDWORK )