162 SUBROUTINE dlarft( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
169 CHARACTER DIRECT, STOREV
170 INTEGER K, LDT, LDV, N
173 DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * )
179 DOUBLE PRECISION ONE, ZERO
180 parameter( one = 1.0d+0, zero = 0.0d+0 )
183 INTEGER I, J, PREVLASTV, LASTV
199 IF( lsame( direct,
'F' ) )
THEN
202 prevlastv = max( i, prevlastv )
203 IF( tau( i ).EQ.zero )
THEN
214 IF( lsame( storev,
'C' ) )
THEN
216 DO lastv = n, i+1, -1
217 IF( v( lastv, i ).NE.zero )
EXIT
220 t( j, i ) = -tau( i ) * v( i , j )
222 j = min( lastv, prevlastv )
226 CALL dgemv(
'Transpose', j-i, i-1, -tau( i ),
227 $ v( i+1, 1 ), ldv, v( i+1, i ), 1, one,
231 DO lastv = n, i+1, -1
232 IF( v( i, lastv ).NE.zero )
EXIT
235 t( j, i ) = -tau( i ) * v( j , i )
237 j = min( lastv, prevlastv )
241 CALL dgemv(
'No transpose', i-1, j-i, -tau( i ),
242 $ v( 1, i+1 ), ldv, v( i, i+1 ), ldv, one,
248 CALL dtrmv(
'Upper',
'No transpose',
'Non-unit', i-1, t,
249 $ ldt, t( 1, i ), 1 )
252 prevlastv = max( prevlastv, lastv )
261 IF( tau( i ).EQ.zero )
THEN
273 IF( lsame( storev,
'C' ) )
THEN
276 IF( v( lastv, i ).NE.zero )
EXIT
279 t( j, i ) = -tau( i ) * v( n-k+i , j )
281 j = max( lastv, prevlastv )
285 CALL dgemv(
'Transpose', n-k+i-j, k-i, -tau( i ),
286 $ v( j, i+1 ), ldv, v( j, i ), 1, one,
291 IF( v( i, lastv ).NE.zero )
EXIT
294 t( j, i ) = -tau( i ) * v( j, n-k+i )
296 j = max( lastv, prevlastv )
300 CALL dgemv(
'No transpose', k-i, n-k+i-j,
301 $ -tau( i ), v( i+1, j ), ldv, v( i, j ), ldv,
302 $ one, t( i+1, i ), 1 )
307 CALL dtrmv(
'Lower',
'No transpose',
'Non-unit', k-i,
308 $ t( i+1, i+1 ), ldt, t( i+1, i ), 1 )
310 prevlastv = min( prevlastv, lastv )
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH