1 SUBROUTINE pbdtrnv( ICONTXT, XDIST, TRANS, N, NB, NZ, X, INCX,
2 $ BETA, Y, INCY, IXROW, IXCOL, IYROW, IYCOL,
14 CHARACTER*1 TRANS, XDIST
15 INTEGER ICONTXT, INCX, INCY, IXCOL, IXROW, IYCOL,
20 DOUBLE PRECISION WORK( * ), X( * ), Y( * )
169 DOUBLE PRECISION ONE, ZERO
170 PARAMETER ( ONE = 1.0d+0, zero = 0.0d+0 )
173 LOGICAL COLFORM, ROWFORM
174 INTEGER I, IDEX, IGD, INFO, JDEX, JYCOL, JYROW, JZ, KZ,
175 $ lcm, lcmp, lcmq, mccol, mcrow, mrcol, mrrow,
176 $ mycol, myrow, nn, np, np0, np1, npcol, nprow,
178 DOUBLE PRECISION TBETA
182 INTEGER ILCM, ICEIL, NUMROC
183 EXTERNAL lsame, ilcm, iceil, numroc
186 EXTERNAL blacs_gridinfo, dgebr2d, dgebs2d, dgerv2d,
199 CALL blacs_gridinfo( icontxt, nprow, npcol, myrow, mycol )
201 colform = lsame( xdist,
'C' )
202 rowform = lsame( xdist,
'R' )
207 IF( ( .NOT.colform ) .AND. ( .NOT.rowform ) )
THEN
209 ELSE IF( n .LT.0 )
THEN
211 ELSE IF( nb .LT.1 )
THEN
213 ELSE IF( nz .LT.0 .OR. nz.GE.nb )
THEN
215 ELSE IF( incx.EQ.0 )
THEN
217 ELSE IF( incy.EQ.0 )
THEN
219 ELSE IF( ixrow.LT.-1 .OR. ixrow.GE.nprow .OR.
220 $ ( ixrow.EQ.-1 .AND. colform ) )
THEN
222 ELSE IF( ixcol.LT.-1 .OR. ixcol.GE.npcol .OR.
223 $ ( ixcol.EQ.-1 .AND. rowform ) )
THEN
225 ELSE IF( iyrow.LT.-1 .OR. iyrow.GE.nprow .OR.
226 $ ( iyrow.EQ.-1 .AND. rowform ) )
THEN
228 ELSE IF( iycol.LT.-1 .OR. iycol.GE.npcol .OR.
229 $ ( iycol.EQ.-1 .AND. colform ) )
THEN
235 CALL pxerbla( icontxt,
'PBDTRNV ', info )
243 lcm = ilcm( nprow, npcol )
263 IF( ixrow.LT.0 .OR. ixrow.GE.nprow )
THEN
265 ELSE IF( ixcol.LT.-1 .OR. ixcol.GE.npcol )
THEN
267 ELSE IF( iyrow.LT.-1 .OR. iyrow.GE.nprow )
THEN
269 ELSE IF( iycol.LT.0 .OR. iycol.GE.npcol )
THEN
272 IF( info.NE.0 )
GO TO 10
277 mrrow = mod( nprow+myrow-ixrow, nprow )
278 mrcol = mod( npcol+mycol-iycol, npcol )
280 IF( iyrow.EQ.-1 ) jyrow = ixrow
282 np = numroc( nn, nb, myrow, ixrow, nprow )
283 IF( mrrow.EQ.0 ) np = np - nz
284 nq = numroc( nn, nb, mycol, iycol, npcol )
285 IF( mrcol.EQ.0 ) nq = nq - nz
286 nq0 = numroc( numroc(nn, nb, 0, 0, npcol), nb, 0, 0, lcmq )
290 IF( ixcol .GE. 0 )
THEN
292 IF( myrow.EQ.jyrow ) tbeta = beta
295 DO 20 i = 0,
min( lcm, iceil(nn,nb) ) - 1
296 mcrow = mod( mod(i, nprow) + ixrow, nprow )
297 mccol = mod( mod(i, npcol) + iycol, npcol )
298 IF( lcmq.EQ.1 ) nq0 = numroc( nn, nb, i, 0, npcol )
299 jdex = (i/npcol) * nb
300 IF( mrcol.EQ.0 ) jdex =
max(0, jdex-nz)
304 IF( myrow.EQ.mcrow .AND. mycol.EQ.ixcol )
THEN
308 idex = (i/nprow) * nb
309 IF( mrrow.EQ.0 ) idex =
max( 0, idex-nz )
310 IF( myrow.EQ.jyrow .AND. mycol.EQ.mccol )
THEN
311 CALL pbdtr2b1( icontxt, trans, np-idex, nb, kz,
312 $ x(idex*incx+1), incx, tbeta,
313 $ y(jdex*incy+1), incy, lcmp, lcmq )
318 CALL pbdtr2b1( icontxt, trans, np-idex, nb, kz,
319 $ x(idex*incx+1), incx, zero, work, 1,
321 CALL dgesd2d( icontxt, 1, nq0-kz, work, 1,
327 ELSE IF( myrow.EQ.jyrow .AND. mycol.EQ.mccol )
THEN
328 IF( lcmq.EQ.1 .AND. tbeta.EQ.zero )
THEN
329 CALL dgerv2d( icontxt, 1, nq0-kz, y, incy,
332 CALL dgerv2d( icontxt, 1, nq0-kz, work, 1,
334 CALL pbdtr2a1( icontxt, nq-jdex, nb, kz, work, 1, tbeta,
335 $ y(jdex*incy+1), incy, lcmq*nb )
343 IF( iyrow.EQ.-1 )
THEN
344 IF( myrow.EQ.jyrow )
THEN
345 CALL dgebs2d( icontxt,
'Col',
'1-tree', 1, nq, y, incy )
347 CALL dgebr2d( icontxt,
'Col',
'1-tree', 1, nq, y, incy,
355 IF( lcmq.EQ.1 ) nq0 = nq
361 IF( mrrow.EQ.0 ) kz = nz
363 IF( mrrow.EQ.0 .AND. mycol.EQ.iycol ) jz = nz
365 DO 30 i = 0, lcmp - 1
366 IF( mrcol.EQ.mod(nprow*i+mrrow, npcol) )
THEN
367 idex =
max( 0, i*nb-kz )
368 IF( lcmq.EQ.1 .AND. (iyrow.EQ.-1.OR.iyrow.EQ.myrow) )
THEN
369 CALL pbdtr2b1( icontxt, trans, np-idex, nb, jz,
370 $ x(idex*incx+1), incx, beta, y, incy,
373 CALL pbdtr2b1( icontxt, trans, np-idex, nb, jz,
374 $ x(idex*incx+1), incx, zero, work, 1,
382 mcrow = mod( mod(mrcol, nprow) + ixrow, nprow )
384 mccol = mod( npcol+mycol-iycol, npcol )
385 CALL pbdtrget( icontxt,
'Row', 1, nq0, iceil( nn, nb ),
386 $ work, 1, mcrow, mccol, igd, myrow, mycol,
392 IF( iyrow.EQ.-1 )
THEN
393 IF( myrow.EQ.mcrow )
THEN
396 IF( mycol.EQ.iycol ) kz = nz
397 CALL pbdtrst1( icontxt,
'Row', nq, nb, kz, work, 1,
398 $ beta, y, incy, lcmp, lcmq, nq0 )
400 CALL dgebs2d( icontxt,
'Col',
'1-tree', 1, nq, y, incy )
402 CALL dgebr2d( icontxt,
'Col',
'1-tree', 1, nq, y, incy,
410 IF( myrow.EQ.mcrow )
THEN
412 $
CALL dgesd2d( icontxt, 1, nq0, work, 1, iyrow, mycol )
413 ELSE IF( myrow.EQ.iyrow )
THEN
414 IF( beta.EQ.zero )
THEN
415 CALL dgerv2d( icontxt, 1, nq0, y, incy, mcrow, mycol )
417 CALL dgerv2d( icontxt, 1, nq0, work, 1, mcrow, mycol )
418 CALL pbdvecadd( icontxt,
'G', nq0, one, work, 1,
424 nq1 = nq0 *
min( lcmq,
max( 0, iceil(nn,nb)-mccol ) )
425 IF( myrow.EQ.mcrow )
THEN
427 $
CALL dgesd2d( icontxt, 1, nq1, work, 1, iyrow, mycol )
428 ELSE IF( myrow.EQ.iyrow )
THEN
429 CALL dgerv2d( icontxt, 1, nq1, work, 1, mcrow, mycol )
432 IF( myrow.EQ.iyrow )
THEN
434 IF( mycol.EQ.iycol ) kz = nz
435 CALL pbdtrst1( icontxt,
'Row', nq, nb, kz, work, 1,
436 $ beta, y, incy, lcmp, lcmq, nq0 )
456 IF( ixrow.LT.-1 .OR. ixrow.GE.nprow )
THEN
458 ELSE IF( ixcol.LT.0 .OR. ixcol.GE.npcol )
THEN
460 ELSE IF( iyrow.LT.0 .OR. iyrow.GE.nprow )
THEN
462 ELSE IF( iycol.LT.-1 .OR. iycol.GE.npcol )
THEN
465 IF( info.NE.0 )
GO TO 10
470 mrrow = mod( nprow+myrow-iyrow, nprow )
471 mrcol = mod( npcol+mycol-ixcol, npcol )
473 IF( iycol.EQ.-1 ) jycol = ixcol
475 np = numroc( nn, nb, myrow, iyrow, nprow )
476 IF( mrrow.EQ.0 ) np = np - nz
477 nq = numroc( nn, nb, mycol, ixcol, npcol )
478 IF( mrcol.EQ.0 ) nq = nq - nz
479 np0 = numroc( numroc(nn, nb, 0, 0, nprow), nb, 0, 0, lcmp )
483 IF( ixrow .GE. 0 )
THEN
485 IF( mycol.EQ.jycol ) tbeta = beta
488 DO 40 i = 0,
min( lcm, iceil(nn,nb) ) - 1
489 mcrow = mod( mod(i, nprow) + iyrow, nprow )
490 mccol = mod( mod(i, npcol) + ixcol, npcol )
491 IF( lcmp.EQ.1 ) np0 = numroc( nn, nb, i, 0, nprow )
492 jdex = (i/nprow) * nb
493 IF( mrrow.EQ.0 ) jdex =
max(0, jdex-nz)
497 IF( myrow.EQ.ixrow .AND. mycol.EQ.mccol )
THEN
501 idex = (i/npcol) * nb
502 IF( mrcol.EQ.0 ) idex =
max( 0, idex-nz )
503 IF( myrow.EQ.mcrow .AND. mycol.EQ.jycol )
THEN
504 CALL pbdtr2b1( icontxt, trans, nq-idex, nb, kz,
505 $ x(idex*incx+1), incx, tbeta,
506 $ y(jdex*incy+1), incy, lcmq, lcmp )
511 CALL pbdtr2b1( icontxt, trans, nq-idex, nb, kz,
512 $ x(idex*incx+1), incx, zero, work, 1,
514 CALL dgesd2d( icontxt, 1, np0-kz, work, 1,
520 ELSE IF( myrow.EQ.mcrow .AND. mycol.EQ.jycol )
THEN
521 IF( lcmp.EQ.1 .AND. tbeta.EQ.zero )
THEN
522 CALL dgerv2d( icontxt, 1, np0-kz, y, incy,
525 CALL dgerv2d( icontxt, 1, np0-kz, work, 1,
527 CALL pbdtr2a1( icontxt, np-jdex, nb, kz, work, 1, tbeta,
528 $ y(jdex*incy+1), incy, lcmp*nb )
536 IF( iycol.EQ.-1 )
THEN
537 IF( mycol.EQ.jycol )
THEN
538 CALL dgebs2d( icontxt,
'Row',
'1-tree', 1, np, y, incy )
540 CALL dgebr2d( icontxt,
'Row',
'1-tree', 1, np, y, incy,
548 IF( lcmp.EQ.1 ) np0 = np
554 IF( mrcol.EQ.0 ) kz = nz
556 IF( mrcol.EQ.0 .AND. myrow.EQ.iyrow ) jz = nz
559 IF( mrrow.EQ.mod(npcol*i+mrcol, nprow) )
THEN
560 idex =
max( 0, i*nb-kz )
561 IF( lcmp.EQ.1 .AND. (iycol.EQ.-1.OR.iycol.EQ.mycol) )
THEN
562 CALL pbdtr2b1( icontxt, trans, nq-idex, nb, jz,
563 $ x(idex*incx+1), incx, beta, y, incy,
566 CALL pbdtr2b1( icontxt, trans, nq-idex, nb, jz,
567 $ x(idex*incx+1), incx, zero, work, 1,
575 mccol = mod( mod(mrrow, npcol) + ixcol, npcol )
577 mcrow = mod( nprow+myrow-iyrow, nprow )
578 CALL pbdtrget( icontxt,
'Col', 1, np0, iceil( nn, nb ),
579 $ work, 1, mcrow, mccol, igd, myrow, mycol,
585 IF( iycol.EQ.-1 )
THEN
586 IF( mycol.EQ.mccol )
THEN
589 IF( myrow.EQ.iyrow ) kz = nz
590 CALL pbdtrst1( icontxt,
'Col', np, nb, kz, work, 1,
591 $ beta, y, incy, lcmp, lcmq, np0 )
593 CALL dgebs2d( icontxt,
'Row',
'1-tree', 1, np, y, incy )
595 CALL dgebr2d( icontxt,
'Row',
'1-tree', 1, np, y, incy,
603 IF( mycol.EQ.mccol )
THEN
605 $
CALL dgesd2d( icontxt, 1, np, work, 1, myrow, iycol )
606 ELSE IF( mycol.EQ.iycol )
THEN
607 IF( beta.EQ.zero )
THEN
608 CALL dgerv2d( icontxt, 1, np, y, incy, myrow, mccol )
610 CALL dgerv2d( icontxt, 1, np, work, 1, myrow, mccol )
611 CALL pbdvecadd( icontxt,
'G', np, one, work, 1, beta,
617 np1 = np0 *
min( lcmp,
max( 0, iceil(nn,nb)-mcrow ) )
618 IF( mycol.EQ.mccol )
THEN
620 $
CALL dgesd2d( icontxt, 1, np1, work, 1, myrow, iycol )
621 ELSE IF( mycol.EQ.iycol )
THEN
622 CALL dgerv2d( icontxt, 1, np1, work, 1, myrow, mccol )
625 IF( mycol.EQ.iycol )
THEN
627 IF( myrow.EQ.iyrow ) kz = nz
628 CALL pbdtrst1( icontxt,
'Col', np, nb, kz, work, 1,
629 $ beta, y, incy, lcmp, lcmq, np0 )
646 SUBROUTINE pbdtr2a1( ICONTXT, N, NB, NZ, X, INCX, BETA, Y, INCY,
654 INTEGER ICONTXT, N, NB, NZ, INCX, INCY, INTV
655 DOUBLE PRECISION BETA
658 DOUBLE PRECISION X( * ), Y( * )
680 PARAMETER ( ONE = 1.0d+0 )
683 INTEGER IX, IY, JZ, K, ITER
688 iter = iceil( n+nz, intv )
691 CALL pbdvecadd( icontxt,
'G', nb-jz, one, x(ix*incx+1), incx,
692 $ beta, y(iy*incy+1), incy )
698 CALL pbdvecadd( icontxt,
'G', nb, one, x(ix*incx+1), incx,
699 $ beta, y(iy*incy+1), incy )
706 $ x(ix*incx+1), incx, beta, y(iy*incy+1), incy )
718 SUBROUTINE pbdtr2b1( ICONTXT, TRANS, N, NB, NZ, X, INCX, BETA, Y,
727 INTEGER ICONTXT, N, NB, NZ, INCX, INCY, JINX, JINY
728 DOUBLE PRECISION BETA
731 DOUBLE PRECISION X( * ), Y( * )
753 parameter( one = 1.0d+0 )
756 INTEGER IX, IY, JZ, K, ITER, LENX, LENY
758 IF( jinx.EQ.1 .AND. jiny.EQ.1 )
THEN
759 CALL pbdvecadd( icontxt, trans, n, one, x, incx, beta,
768 iter = iceil( n+nz, lenx )
771 CALL pbdvecadd( icontxt, trans, nb-jz, one, x(ix*incx+1),
772 $ incx, beta, y(iy*incy+1), incy )
778 CALL pbdvecadd( icontxt, trans, nb, one, x(ix*incx+1),
779 $ incx, beta, y(iy*incy+1), incy )
785 CALL pbdvecadd( icontxt, trans,
min( n-ix, nb-jz ), one,
786 $ x(ix*incx+1), incx, beta, y(iy*incy+1), incy )