1 SUBROUTINE pzgebdrv( M, N, A, IA, JA, DESCA, D, E, TAUQ, TAUP,
10 INTEGER INFO, IA, JA, M, N
14 DOUBLE PRECISION D( * ), E( * )
15 COMPLEX*16 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 )
158 DOUBLE PRECISION REIGHT, RZERO
159 parameter( reight = 8.0d+0, rzero = 0.0d+0 )
161 parameter( one = ( 1.0d+0, 0.0d+0 ),
162 $ zero = ( 0.0d+0, 0.0d+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
168 DOUBLE PRECISION ADDBND, D2, E2
172 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), DESCV( DLEN_ ),
181 INTEGER INDXG2P, NUMROC
182 DOUBLE PRECISION PDLAMCH
183 EXTERNAL indxg2p, numroc, pdlamch
186 INTRINSIC abs, dcmplx,
max,
min, mod
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 * pdlamch( 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 pdelget(
' ',
' ', d2, d, 1, ja+j, descd )
238 CALL pzelget(
'Columnwise',
' ', d1, a, ia+j, ja+j, desca )
239 IF( j.LT.(mn-1) )
THEN
240 CALL pdelget(
' ',
' ', e2, e, ia+j, 1, desce )
241 CALL pzelget(
'Rowwise',
' ', e1, a, ia+j, ja+j+1,
245 IF( ( abs( d1-dcmplx( d2 ) ).GT.( abs( d2 )*addbnd ) ) .OR.
246 $ ( abs( e1-dcmplx( e2 ) ).GT.( abs( e2 )*addbnd ) ) )
250 DO 20 j = jl, ja+nb-ioff, -nb
251 jb =
min( ja+n-j, nb )
257 CALL pzlarft(
'Forward',
'Columnwise', m-k+1, jb, a, i, j,
258 $ desca, tauq, work( iptq ), work( ipwk ) )
262 CALL pzlacpy(
'Lower', m-k+1, jb, a, i, j, desca,
263 $ work( ipv ), iv, jv, descv )
264 CALL pzlaset(
'Upper', m-k+1, jb, zero, one, work( ipv ),
269 CALL pzlaset(
'Lower', m-k, jb, zero, zero, a, i+1, j,
274 CALL pzlarft(
'Forward',
'Rowwise', n-k, jb, a, i, j+1,
275 $ desca, taup, work( iptp ), work( ipwk ) )
279 CALL pzlacpy(
'Upper', jb, n-k, a, i, j+1, desca,
280 $ work( ipw ), iv, jv+1, descw )
281 CALL pzlaset(
'Lower', jb, n-k, zero, one, work( ipw ), iv,
286 CALL pzlaset(
'Upper', jb, n-k-1, zero, zero, a, i, j+2,
291 CALL pzlarfb(
'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 pzlarfb(
'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 pzlarft(
'Forward',
'Columnwise', m, jb, a, ia, ja, desca,
321 $ tauq, work( iptq ), work( ipwk ) )
325 CALL pzlacpy(
'Lower', m, jb, a, ia, ja, desca, work( ipv ),
327 CALL pzlaset(
'Upper', m, jb, zero, one, work( ipv ), iv, jv,
332 CALL pzlaset(
'Lower', m-1, jb, zero, zero, a, ia+1, ja,
337 CALL pzlarft(
'Forward',
'Rowwise', n-1, jb, a, ia, ja+1,
338 $ desca, taup, work( iptp ), work( ipwk ) )
342 CALL pzlacpy(
'Upper', jb, n-1, a, ia, ja+1, desca,
343 $ work( ipw ), iv, jv+1, descw )
344 CALL pzlaset(
'Lower', jb, n-1, zero, one, work( ipw ), iv,
349 CALL pzlaset(
'Upper', jb, n-2, zero, zero, a, ia, ja+2,
354 CALL pzlarfb(
'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 pzlarfb(
'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 pdelget(
' ',
' ', d2, d, ia+j, 1, descd )
379 CALL pzelget(
'Rowwise',
' ', d1, a, ia+j, ja+j, desca )
380 IF( j.LT.(mn-1) )
THEN
381 CALL pdelget(
' ',
' ', e2, e, 1, ja+j, desce )
382 CALL pzelget(
'Columnwise',
' ', e1, a, ia+j+1, ja+j,
386 IF( ( abs( d1-dcmplx( d2 ) ).GT.( abs( d2 )*addbnd ) ) .OR.
387 $ ( abs( e1-dcmplx( e2 ) ).GT.( abs( e2 )*addbnd ) ) )
391 DO 40 i = il, ia+nb-ioff, -nb
392 jb =
min( ia+m-i, nb )
398 CALL pzlarft(
'Forward',
'Columnwise', m-k, jb, a, i+1, j,
399 $ desca, tauq, work( iptq ), work( ipwk ) )
403 CALL pzlacpy(
'Lower', m-k, jb, a, i+1, j, desca,
404 $ work( ipv ), iv+1, jv, descv )
405 CALL pzlaset(
'Upper', m-k, jb, zero, one, work( ipv ),
410 CALL pzlaset(
'Lower', m-k-1, jb, zero, zero, a, i+2, j,
415 CALL pzlarft(
'Forward',
'Rowwise', n-k+1, jb, a, i, j,
416 $ desca, taup, work( iptp ), work( ipwk ) )
420 CALL pzlacpy(
'Upper', jb, n-k+1, a, i, j, desca,
421 $ work( ipw ), iv, jv, descw )
422 CALL pzlaset(
'Lower', jb, n-k+1, zero, one, work( ipw ),
427 CALL pzlaset(
'Upper', jb, n-k, zero, zero, a, i, j+1,
432 CALL pzlarfb(
'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 pzlarfb(
'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 pzlarft(
'Forward',
'Columnwise', m-1, jb, a, ia+1, ja,
462 $ desca, tauq, work( iptq ), work( ipwk ) )
466 CALL pzlacpy(
'Lower', m-1, jb, a, ia+1, ja, desca,
467 $ work( ipv ), iv+1, jv, descv )
468 CALL pzlaset(
'Upper', m-1, jb, zero, one, work( ipv ), iv+1,
473 CALL pzlaset(
'Lower', m-2, jb, zero, zero, a, ia+2, ja,
478 CALL pzlarft(
'Forward',
'Rowwise', n, jb, a, ia, ja, desca,
479 $ taup, work( iptp ), work( ipwk ) )
483 CALL pzlacpy(
'Upper', jb, n, a, ia, ja, desca, work( ipw ),
485 CALL pzlaset(
'Lower', jb, n, zero, one, work( ipw ), iv, jv,
490 CALL pzlaset(
'Upper', jb, n-1, zero, zero, a, ia, ja+1,
495 CALL pzlarfb(
'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 pzlarfb(
'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 )