1 SUBROUTINE pcgebdrv( M, N, A, IA, JA, DESCA, D, E, TAUQ, TAUP,
10 INTEGER INFO, IA, JA, M, N
15 COMPLEX A( * ), TAUP( * ), TAUQ( * ), WORK( * )
153 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
154 $ lld_, mb_, m_, nb_, n_, rsrc_
155 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
156 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
157 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
159 parameter( reight = 8.0e+0, rzero = 0.0e+0 )
161 parameter( one = ( 1.0e+0, 0.0e+0 ),
162 $ zero = ( 0.0e+0, 0.0e+0 ) )
165 INTEGER I, IACOL, IAROW, ICTXT, IIA, IL, IPTP, IPTQ,
166 $ ipv, ipw, ipwk, ioff, iv, j, jb, jja, jl, jv,
167 $ k, mn, mp, mycol, myrow, nb, npcol, nprow, nq
172 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), DESCV( DLEN_ ),
181 INTEGER INDXG2P, NUMROC
183 EXTERNAL indxg2p, numroc, pslamch
192 ictxt = desca( ctxt_ )
193 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
197 ioff = mod( ia-1, nb )
198 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
200 mp = numroc( m+ioff, nb, myrow, iarow, nprow )
201 nq = numroc( n+ioff, nb, mycol, iacol, npcol )
211 il =
max( ( (ia+mn-2) / nb )*nb + 1, ia )
212 jl =
max( ( (ja+mn-2) / nb )*nb + 1, ja )
213 iarow = indxg2p( il, nb, myrow, desca( rsrc_ ), nprow )
214 iacol = indxg2p( jl, nb, mycol, desca( csrc_ ), npcol )
215 CALL descset( descv, ia+m-il, nb, nb, nb, iarow, iacol, ictxt,
217 CALL descset( descw, nb, ja+n-jl, nb, nb, iarow, iacol, ictxt,
220 addbnd = reight * pslamch( ictxt,
'eps' )
226 CALL descset( descd, 1, ja+mn-1, 1, desca( nb_ ), myrow,
227 $ desca( csrc_ ), desca( ctxt_ ), 1 )
228 CALL descset( desce, ia+mn-1, 1, desca( mb_ ), 1,
229 $ desca( rsrc_ ), mycol, desca( ctxt_ ),
237 CALL pselget(
' ',
' ', d2, d, 1, ja+j, descd )
238 CALL pcelget(
'Columnwise',
' ', d1, a, ia+j, ja+j, desca )
239 IF( j.LT.(mn-1) )
THEN
240 CALL pselget(
' ',
' ', e2, e, ia+j, 1, desce )
241 CALL pcelget(
'Rowwise',
' ', e1, a, ia+j, ja+j+1,
245 IF( ( abs( d1 -
cmplx( d2 ) ).GT.( abs( d2 )*addbnd ) ) .OR.
246 $ ( abs( e1 -
cmplx( e2 ) ).GT.( abs( e2 )*addbnd ) ) )
250 DO 20 j = jl, ja+nb-ioff, -nb
251 jb =
min( ja+n-j, nb )
257 CALL pclarft(
'Forward',
'Columnwise', m-k+1, jb, a, i, j,
258 $ desca, tauq, work( iptq ), work( ipwk ) )
262 CALL pclacpy(
'Lower', m-k+1, jb, a, i, j, desca,
263 $ work( ipv ), iv, jv, descv )
264 CALL pclaset(
'Upper', m-k+1, jb, zero, one, work( ipv ),
269 CALL pclaset(
'Lower', m-k, jb, zero, zero, a, i+1, j,
274 CALL pclarft(
'Forward',
'Rowwise', n-k, jb, a, i, j+1,
275 $ desca, taup, work( iptp ), work( ipwk ) )
279 CALL pclacpy(
'Upper', jb, n-k, a, i, j+1, desca,
280 $ work( ipw ), iv, jv+1, descw )
281 CALL pclaset(
'Lower', jb, n-k, zero, one, work( ipw ), iv,
286 CALL pclaset(
'Upper', jb, n-k-1, zero, zero, a, i, j+2,
291 CALL pclarfb(
'Left',
'No transpose',
'Forward',
292 $
'Columnwise', m-k+1, n-k+1, jb, work( ipv ),
293 $ iv, jv, descv, work( iptq ), a, i, j, desca,
298 CALL pclarfb(
'Right',
'Conjugate transpose',
'Forward',
299 $
'Rowwise', m-k+1, n-k, jb, work( ipw ), iv,
300 $ jv+1, descw, work( iptp ), a, i, j+1, desca,
303 descv( m_ ) = descv( m_ ) + nb
304 descv( rsrc_ ) = mod( descv( rsrc_ ) + nprow - 1, nprow )
305 descv( csrc_ ) = mod( descv( csrc_ ) + npcol - 1, npcol )
306 descw( n_ ) = descw( n_ ) + nb
307 descw( rsrc_ ) = descv( rsrc_ )
308 descw( csrc_ ) = descv( csrc_ )
314 jb =
min( n, nb - ioff )
320 CALL pclarft(
'Forward',
'Columnwise', m, jb, a, ia, ja, desca,
321 $ tauq, work( iptq ), work( ipwk ) )
325 CALL pclacpy(
'Lower', m, jb, a, ia, ja, desca, work( ipv ),
327 CALL pclaset(
'Upper', m, jb, zero, one, work( ipv ), iv, jv,
332 CALL pclaset(
'Lower', m-1, jb, zero, zero, a, ia+1, ja,
337 CALL pclarft(
'Forward',
'Rowwise', n-1, jb, a, ia, ja+1,
338 $ desca, taup, work( iptp ), work( ipwk ) )
342 CALL pclacpy(
'Upper', jb, n-1, a, ia, ja+1, desca,
343 $ work( ipw ), iv, jv+1, descw )
344 CALL pclaset(
'Lower', jb, n-1, zero, one, work( ipw ), iv,
349 CALL pclaset(
'Upper', jb, n-2, zero, zero, a, ia, ja+2,
354 CALL pclarfb(
'Left',
'No transpose',
'Forward',
'Columnwise',
355 $ m, n, jb, work( ipv ), iv, jv, descv,
356 $ work( iptq ), a, ia, ja, desca, work( ipwk ) )
360 CALL pclarfb(
'Right',
'Conjugate transpose',
'Forward',
361 $
'Rowwise', m, n-1, jb, work( ipw ), iv, jv+1,
362 $ descw, work( iptp ), a, ia, ja+1, desca,
367 CALL descset( descd, ia+mn-1, 1, desca( mb_ ), 1,
368 $ desca( rsrc_ ), mycol, desca( ctxt_ ),
370 CALL descset( desce, 1, ja+mn-2, 1, desca( nb_ ), myrow,
371 $ desca( csrc_ ), desca( ctxt_ ), 1 )
378 CALL pselget(
' ',
' ', d2, d, ia+j, 1, descd )
379 CALL pcelget(
'Rowwise',
' ', d1, a, ia+j, ja+j, desca )
380 IF( j.LT.(mn-1) )
THEN
381 CALL pselget(
' ',
' ', e2, e, 1, ja+j, desce )
382 CALL pcelget(
'Columnwise',
' ', e1, a, ia+j+1, ja+j,
386 IF( ( abs( d1 -
cmplx( d2 ) ).GT.( abs( d2 )*addbnd ) ) .OR.
387 $ ( abs( e1 -
cmplx( e2 ) ).GT.( abs( e2 )*addbnd ) ) )
391 DO 40 i = il, ia+nb-ioff, -nb
392 jb =
min( ia+m-i, nb )
398 CALL pclarft(
'Forward',
'Columnwise', m-k, jb, a, i+1, j,
399 $ desca, tauq, work( iptq ), work( ipwk ) )
403 CALL pclacpy(
'Lower', m-k, jb, a, i+1, j, desca,
404 $ work( ipv ), iv+1, jv, descv )
405 CALL pclaset(
'Upper', m-k, jb, zero, one, work( ipv ),
410 CALL pclaset(
'Lower', m-k-1, jb, zero, zero, a, i+2, j,
415 CALL pclarft(
'Forward',
'Rowwise', n-k+1, jb, a, i, j,
416 $ desca, taup, work( iptp ), work( ipwk ) )
420 CALL pclacpy(
'Upper', jb, n-k+1, a, i, j, desca,
421 $ work( ipw ), iv, jv, descw )
422 CALL pclaset(
'Lower', jb, n-k+1, zero, one, work( ipw ),
427 CALL pclaset(
'Upper', jb, n-k, zero, zero, a, i, j+1,
432 CALL pclarfb(
'Left',
'No transpose',
'Forward',
433 $
'Columnwise', m-k, n-k+1, jb, work( ipv ),
434 $ iv+1, jv, descv, work( iptq ), a, i+1, j,
435 $ desca, work( ipwk ) )
439 CALL pclarfb(
'Right',
'Conjugate transpose',
'Forward',
440 $
'Rowwise', m-k+1, n-k+1, jb, work( ipw ), iv,
441 $ jv, descw, work( iptp ), a, i, j, desca,
444 descv( m_ ) = descv( m_ ) + nb
445 descv( rsrc_ ) = mod( descv( rsrc_ ) + nprow - 1, nprow )
446 descv( csrc_ ) = mod( descv( csrc_ ) + npcol - 1, npcol )
447 descw( n_ ) = descw( n_ ) + nb
448 descw( rsrc_ ) = descv( rsrc_ )
449 descw( csrc_ ) = descv( csrc_ )
455 jb =
min( m, nb - ioff )
461 CALL pclarft(
'Forward',
'Columnwise', m-1, jb, a, ia+1, ja,
462 $ desca, tauq, work( iptq ), work( ipwk ) )
466 CALL pclacpy(
'Lower', m-1, jb, a, ia+1, ja, desca,
467 $ work( ipv ), iv+1, jv, descv )
468 CALL pclaset(
'Upper', m-1, jb, zero, one, work( ipv ), iv+1,
473 CALL pclaset(
'Lower', m-2, jb, zero, zero, a, ia+2, ja,
478 CALL pclarft(
'Forward',
'Rowwise', n, jb, a, ia, ja, desca,
479 $ taup, work( iptp ), work( ipwk ) )
483 CALL pclacpy(
'Upper', jb, n, a, ia, ja, desca, work( ipw ),
485 CALL pclaset(
'Lower', jb, n, zero, one, work( ipw ), iv, jv,
490 CALL pclaset(
'Upper', jb, n-1, zero, zero, a, ia, ja+1,
495 CALL pclarfb(
'Left',
'No transpose',
'Forward',
'Columnwise',
496 $ m-1, n, jb, work( ipv ), iv+1, jv, descv,
497 $ work( iptq ), a, ia+1, ja, desca, work( ipwk ) )
501 CALL pclarfb(
'Right',
'Conjugate transpose',
'Forward',
502 $
'Rowwise', m, n, jb, work( ipw ), iv, jv, descw,
503 $ work( iptp ), a, ia, ja, desca, work( ipwk ) )
506 CALL igsum2d( ictxt,
'All',
' ', 1, 1, info, 1, -1, 0 )