#include "ESMF_LapackBlas.inc" !> \brief \b DLARFT ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DLARFT + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarft.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) ! ! .. Scalar Arguments .. ! CHARACTER DIRECT, STOREV ! INTEGER K, LDT, LDV, N ! .. ! .. Array Arguments .. ! DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> DLARFT forms the triangular factor T of a real block reflector H !> of order n, which is defined as a product of k elementary reflectors. !> !> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; !> !> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. !> !> If STOREV = 'C', the vector which defines the elementary reflector !> H(i) is stored in the i-th column of the array V, and !> !> H = I - V * T * V**T !> !> If STOREV = 'R', the vector which defines the elementary reflector !> H(i) is stored in the i-th row of the array V, and !> !> H = I - V**T * T * V !> \endverbatim ! ! Arguments: ! ========== ! !> \param[in] DIRECT !> \verbatim !> DIRECT is CHARACTER*1 !> Specifies the order in which the elementary reflectors are !> multiplied to form the block reflector: !> = '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 !> Specifies how the vectors which define the elementary !> reflectors are stored (see also Further Details): !> = 'C': columnwise !> = 'R': rowwise !> \endverbatim !> !> \param[in] N !> \verbatim !> N is INTEGER !> The order of the block reflector H. N >= 0. !> \endverbatim !> !> \param[in] K !> \verbatim !> K is INTEGER !> The order of the triangular factor T (= the number of !> elementary reflectors). K >= 1. !> \endverbatim !> !> \param[in,out] V !> \verbatim !> V is DOUBLE PRECISION array, dimension !> (LDV,K) if STOREV = 'C' !> (LDV,N) if STOREV = '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', LDV >= max(1,N); if STOREV = 'R', LDV >= K. !> \endverbatim !> !> \param[in] TAU !> \verbatim !> TAU is DOUBLE PRECISION array, dimension (K) !> TAU(i) must contain the scalar factor of the elementary !> reflector H(i). !> \endverbatim !> !> \param[out] T !> \verbatim !> T is DOUBLE PRECISION array, dimension (LDT,K) !> The k by k triangular factor T of the block reflector. !> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is !> lower triangular. The rest of the array is not used. !> \endverbatim !> !> \param[in] LDT !> \verbatim !> LDT is INTEGER !> The leading dimension of the array T. LDT >= K. !> \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 DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) ! ! -- 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 2011 ! ! .. Scalar Arguments .. CHARACTER DIRECT, STOREV INTEGER K, LDT, LDV, N ! .. ! .. Array Arguments .. DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) ! .. ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Local Scalars .. INTEGER I, J, PREVLASTV, LASTV DOUBLE PRECISION VII ! .. ! .. External Subroutines .. EXTERNAL DGEMV, DTRMV ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. Executable Statements .. ! ! Quick return if possible ! IF( N.EQ.0 ) & & RETURN ! IF( LSAME( DIRECT, 'F' ) ) THEN PREVLASTV = N DO 20 I = 1, K PREVLASTV = MAX( I, PREVLASTV ) IF( TAU( I ).EQ.ZERO ) THEN ! ! H(i) = I ! DO 10 J = 1, I T( J, I ) = ZERO 10 CONTINUE ELSE ! ! general case ! VII = V( I, I ) V( I, I ) = ONE IF( LSAME( STOREV, 'C' ) ) THEN ! Skip any trailing zeros. DO LASTV = N, I+1, -1 IF( V( LASTV, I ).NE.ZERO ) EXIT END DO J = MIN( LASTV, PREVLASTV ) ! ! T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) ! CALL DGEMV( 'Transpose', J-I+1, I-1, -TAU( I ), & & V( I, 1 ), LDV, V( I, I ), 1, ZERO, & & T( 1, I ), 1 ) ELSE ! Skip any trailing zeros. DO LASTV = N, I+1, -1 IF( V( I, LASTV ).NE.ZERO ) EXIT END DO J = MIN( LASTV, PREVLASTV ) ! ! T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T ! CALL DGEMV( 'No transpose', I-1, J-I+1, -TAU( I ), & & V( 1, I ), LDV, V( I, I ), LDV, ZERO, & & T( 1, I ), 1 ) END IF V( I, I ) = VII ! ! T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) ! CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, & & LDT, T( 1, I ), 1 ) T( I, I ) = TAU( I ) IF( I.GT.1 ) THEN PREVLASTV = MAX( PREVLASTV, LASTV ) ELSE PREVLASTV = LASTV END IF END IF 20 CONTINUE ELSE PREVLASTV = 1 DO 40 I = K, 1, -1 IF( TAU( I ).EQ.ZERO ) THEN ! ! H(i) = I ! DO 30 J = I, K T( J, I ) = ZERO 30 CONTINUE ELSE ! ! general case ! IF( I.LT.K ) THEN IF( LSAME( STOREV, 'C' ) ) THEN VII = V( N-K+I, I ) V( N-K+I, I ) = ONE ! Skip any leading zeros. DO LASTV = 1, I-1 IF( V( LASTV, I ).NE.ZERO ) EXIT END DO J = MAX( LASTV, PREVLASTV ) ! ! T(i+1:k,i) := ! - tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) ! CALL DGEMV( 'Transpose', N-K+I-J+1, K-I, -TAU( I ), & & V( J, I+1 ), LDV, V( J, I ), 1, ZERO, & & T( I+1, I ), 1 ) V( N-K+I, I ) = VII ELSE VII = V( I, N-K+I ) V( I, N-K+I ) = ONE ! Skip any leading zeros. DO LASTV = 1, I-1 IF( V( I, LASTV ).NE.ZERO ) EXIT END DO J = MAX( LASTV, PREVLASTV ) ! ! T(i+1:k,i) := ! - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T ! CALL DGEMV( 'No transpose', K-I, N-K+I-J+1, & & -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, & & ZERO, T( I+1, I ), 1 ) V( I, N-K+I ) = VII END IF ! ! T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) ! CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, & & T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) IF( I.GT.1 ) THEN PREVLASTV = MIN( PREVLASTV, LASTV ) ELSE PREVLASTV = LASTV END IF END IF T( I, I ) = TAU( I ) END IF 40 CONTINUE END IF RETURN ! ! End of DLARFT ! END