SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) ! ! -- LAPACK auxiliary routine (version 3.1) -- ! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. ! November 2006 ! ! .. Scalar Arguments .. CHARACTER DIRECT, STOREV INTEGER K, LDT, LDV, N ! .. ! .. Array Arguments .. DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) ! .. ! ! Purpose ! ======= ! ! 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' ! ! 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 * V ! ! Arguments ! ========= ! ! DIRECT (input) 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) ! ! STOREV (input) CHARACTER*1 ! Specifies how the vectors which define the elementary ! reflectors are stored (see also Further Details): ! = 'C': columnwise ! = 'R': rowwise ! ! N (input) INTEGER ! The order of the block reflector H. N >= 0. ! ! K (input) INTEGER ! The order of the triangular factor T (= the number of ! elementary reflectors). K >= 1. ! ! V (input/output) DOUBLE PRECISION array, dimension ! (LDV,K) if STOREV = 'C' ! (LDV,N) if STOREV = 'R' ! The matrix V. See further details. ! ! LDV (input) INTEGER ! The leading dimension of the array V. ! If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. ! ! TAU (input) DOUBLE PRECISION array, dimension (K) ! TAU(i) must contain the scalar factor of the elementary ! reflector H(i). ! ! T (output) 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. ! ! LDT (input) INTEGER ! The leading dimension of the array T. LDT >= K. ! ! Further Details ! =============== ! ! 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 ) ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Local Scalars .. INTEGER I, J 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 DO 20 I = 1, K 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 ! ! T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) ! CALL DGEMV( 'Transpose', N-I+1, I-1, -TAU( I ), & V( I, 1 ), LDV, V( I, I ), 1, ZERO, & T( 1, I ), 1 ) ELSE ! ! T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' ! CALL DGEMV( 'No transpose', I-1, N-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 ) END IF 20 CONTINUE ELSE 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 ! ! T(i+1:k,i) := ! - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) ! CALL DGEMV( 'Transpose', N-K+I, K-I, -TAU( I ), & V( 1, I+1 ), LDV, V( 1, 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 ! ! T(i+1:k,i) := ! - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' ! CALL DGEMV( 'No transpose', K-I, N-K+I, -TAU( I ), & V( I+1, 1 ), LDV, V( I, 1 ), 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 ) END IF T( I, I ) = TAU( I ) END IF 40 CONTINUE END IF RETURN ! ! End of DLARFT ! END SUBROUTINE DLARFT