1 SUBROUTINE pdlarft( DIRECT, STOREV, N, K, V, IV, JV, DESCV, TAU,
10 CHARACTER DIRECT, STOREV
15 DOUBLE PRECISION TAU( * ), T( * ), V( * ), WORK( * )
173 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
174 $ lld_, mb_, m_, nb_, n_, rsrc_
175 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
176 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
177 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
178 DOUBLE PRECISION ONE, ZERO
179 parameter( one = 1.0d+0, zero = 0.0d+0 )
183 INTEGER ICOFF, ICTXT, II, IIV, IROFF, IVCOL, IVROW,
184 $ itmp0, itmp1, iw, jj, jjv, ldv, micol, mirow,
185 $ mycol, myrow, np, npcol, nprow, nq
189 EXTERNAL blacs_gridinfo, dcopy, dgemv, dgsum2d,
194 INTEGER INDXG2P, NUMROC
195 EXTERNAL indxg2p, lsame, numroc
204 IF( n.LE.0 .OR. k.LE.0 )
207 ictxt = descv( ctxt_ )
208 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
210 forward = lsame( direct,
'F' )
211 CALL infog2l( iv, jv, descv, nprow, npcol, myrow, mycol,
212 $ iiv, jjv, ivrow, ivcol )
214 IF( lsame( storev,
'C' ) .AND. mycol.EQ.ivcol )
THEN
218 iroff = mod( iv-1, descv( mb_ ) )
224 np = numroc( n+iroff, descv( mb_ ), myrow, ivrow, nprow )
225 IF( myrow.EQ.ivrow )
THEN
231 IF( iroff+1.EQ.descv( mb_ ) )
THEN
232 mirow = mod( ivrow+1, nprow )
238 DO 10 jj = jjv+1, jjv+k-1
240 IF( myrow.EQ.mirow )
THEN
241 vii = v( ii+(jj-1)*ldv )
242 v( ii+(jj-1)*ldv ) = one
249 IF( np-ii+iiv.GT.0 )
THEN
250 CALL dgemv(
'Transpose', np-ii+iiv, itmp0,
251 $ -tau( jj ), v( ii+(jjv-1)*ldv ), ldv,
252 $ v( ii+(jj-1)*ldv ), 1, zero,
255 CALL dlaset(
'All', itmp0, 1, zero, zero, work( iw ),
260 IF( myrow.EQ.mirow )
THEN
261 v( ii+(jj-1)*ldv ) = vii
265 IF( mod( iv+itmp0, descv( mb_ ) ).EQ.0 )
266 $ mirow = mod( mirow+1, nprow )
270 CALL dgsum2d( ictxt,
'Columnwise',
' ', iw-1, 1, work, iw-1,
273 IF( myrow.EQ.ivrow )
THEN
279 t( itmp1 ) = tau( jjv )
281 DO 20 jj = jjv+1, jjv+k-1
286 itmp1 = itmp1 + descv( nb_ )
287 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
290 CALL dtrmv(
'Upper',
'No transpose',
'Non-unit',
291 $ itmp0, t, descv( nb_ ), t( itmp1 ), 1 )
292 t(itmp1+itmp0) = tau( jj )
302 np = numroc( n+iroff-1, descv( mb_ ), myrow, ivrow, nprow )
305 mirow = indxg2p( iv+n-2, descv( mb_ ), myrow,
306 $ descv( rsrc_ ), nprow )
310 DO 30 jj = jjv+k-2, jjv, -1
312 IF( myrow.EQ.mirow )
THEN
313 vii = v( ii+(jj-1)*ldv )
314 v( ii+(jj-1)*ldv ) = one
321 IF( ii-iiv+1.GT.0 )
THEN
322 CALL dgemv(
'Transpose', ii-iiv+1, itmp0, -tau( jj ),
323 $ v( iiv+jj*ldv ), ldv,
324 $ v( iiv+(jj-1)*ldv ), 1, zero,
327 CALL dlaset(
'All', itmp0, 1, zero, zero, work( iw ),
332 IF( myrow.EQ.mirow )
THEN
333 v( ii+(jj-1)*ldv ) = vii
337 IF( mod( iv+n-itmp0-2, descv(mb_) ).EQ.0 )
338 $ mirow = mod( mirow+nprow-1, nprow )
342 CALL dgsum2d( ictxt,
'Columnwise',
' ', iw-1, 1, work, iw-1,
345 IF( myrow.EQ.ivrow )
THEN
349 itmp1 = k + 1 + (k-1) * descv( nb_ )
351 t( itmp1-1 ) = tau( jjv+k-1 )
353 DO 40 jj = jjv+k-2, jjv, -1
358 itmp1 = itmp1 - descv( nb_ ) - 1
359 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
362 CALL dtrmv(
'Lower',
'No transpose',
'Non-unit',
363 $ itmp0, t( itmp1+descv( nb_ ) ),
364 $ descv( nb_ ), t( itmp1 ), 1 )
365 t( itmp1-1 ) = tau( jj )
373 ELSE IF( lsame( storev,
'R' ) .AND. myrow.EQ.ivrow )
THEN
377 icoff = mod( jv-1, descv( nb_ ) )
383 nq = numroc( n+icoff, descv( nb_ ), mycol, ivcol, npcol )
384 IF( mycol.EQ.ivcol )
THEN
390 IF( icoff+1.EQ.descv( nb_ ) )
THEN
391 micol = mod( ivcol+1, npcol )
397 DO 50 ii = iiv+1, iiv+k-1
399 IF( mycol.EQ.micol )
THEN
400 vii = v( ii+(jj-1)*ldv )
401 v( ii+(jj-1)*ldv ) = one
408 IF( nq-jj+jjv.GT.0 )
THEN
409 CALL dgemv(
'No transpose', itmp0, nq-jj+jjv,
410 $ -tau(ii), v( iiv+(jj-1)*ldv ), ldv,
411 $ v( ii+(jj-1)*ldv ), ldv, zero,
414 CALL dlaset(
'All', itmp0, 1, zero, zero,
415 $ work( iw ), itmp0 )
419 IF( mycol.EQ.micol )
THEN
420 v( ii+(jj-1)*ldv ) = vii
424 IF( mod( jv+itmp0, descv( nb_ ) ).EQ.0 )
425 $ micol = mod( micol+1, npcol )
429 CALL dgsum2d( ictxt,
'Rowwise',
' ', iw-1, 1, work, iw-1,
432 IF( mycol.EQ.ivcol )
THEN
438 t( itmp1 ) = tau( iiv )
440 DO 60 ii = iiv+1, iiv+k-1
445 itmp1 = itmp1 + descv( mb_ )
446 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
449 CALL dtrmv(
'Upper',
'No transpose',
'Non-unit',
450 $ itmp0, t, descv( mb_ ), t( itmp1 ), 1 )
451 t( itmp1+itmp0 ) = tau( ii )
461 nq = numroc( n+icoff-1, descv( nb_ ), mycol, ivcol, npcol )
464 micol = indxg2p( jv+n-2, descv( nb_ ), mycol,
465 $ descv( csrc_ ), npcol )
469 DO 70 ii = iiv+k-2, iiv, -1
471 IF( mycol.EQ.micol )
THEN
472 vii = v( ii+(jj-1)*ldv )
473 v( ii+(jj-1)*ldv ) = one
480 IF( jj-jjv+1.GT.0 )
THEN
481 CALL dgemv(
'No transpose', itmp0, jj-jjv+1,
482 $ -tau( ii ), v( ii+1+(jjv-1)*ldv ), ldv,
483 $ v( ii+(jjv-1)*ldv ), ldv, zero,
486 CALL dlaset(
'All', itmp0, 1, zero, zero,
487 $ work( iw ), itmp0 )
491 IF( mycol.EQ.micol )
THEN
492 v( ii+(jj-1)*ldv ) = vii
496 IF( mod( jv+n-itmp0-2, descv( nb_ ) ).EQ.0 )
497 $ micol = mod( micol+npcol-1, npcol )
501 CALL dgsum2d( ictxt,
'Rowwise',
' ', iw-1, 1, work, iw-1,
504 IF( mycol.EQ.ivcol )
THEN
508 itmp1 = k + 1 + (k-1) * descv( mb_ )
510 t( itmp1-1 ) = tau( iiv+k-1 )
512 DO 80 ii = iiv+k-2, iiv, -1
517 itmp1 = itmp1 - descv( mb_ ) - 1
518 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
521 CALL dtrmv(
'Lower',
'No transpose',
'Non-unit',
522 $ itmp0, t( itmp1+descv( mb_ ) ),
523 $ descv( mb_ ), t( itmp1 ), 1 )
524 t( itmp1-1 ) = tau( ii )