SUBROUTINE PBDGEMV( ICONTXT, TRANS, XDIST, YDIST, M, N, MB, NB, $ MZ, NZ, ALPHA, A, LDA, X, INCX, BETA, Y, INCY, $ IAROW, IACOL, IXROW, IXCOL, IYROW, IYCOL, $ XCOMM, XWORK, YWORK, WORK ) * * -- PB-BLAS routine (version 2.1) -- * University of Tennessee, Knoxville, Oak Ridge National Laboratory. * April 28, 1996 * * Jaeyoung Choi, Oak Ridge National Laboratory * Jack Dongarra, University of Tennessee and Oak Ridge National Lab. * David Walker, Oak Ridge National Laboratory * * .. Scalar Arguments .. CHARACTER*1 TRANS, XCOMM, XDIST, XWORK, YDIST, YWORK INTEGER IACOL, IAROW, ICONTXT, INCX, INCY, IXCOL, $ IXROW, IYCOL, IYROW, LDA, M, MB, MZ, N, NB, NZ DOUBLE PRECISION ALPHA, BETA * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ), WORK( * ) * .. * * Purpose * ======= * * PBDGEMV is a parallel blocked version of DGEMV. * PBDGEMV performs one of the matrix-vector operations based on block * cyclic distribution. * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * M-by-N matrix. * * The first element of the matrices A is located in the middle of * the first block ((MZ+1,NZ+1) position). If TRANS = `N', the first * elements of X and Y start from (NZ+1)-th position and (MZ+1)-th * position, respectively, otherwise (MZ+1)-th position and (NZ+1)-th * position, respectively. * * X is broadcast columnwise or rowwise if necessary, and the resultant * Y is collected. * * Parameters * ========== * * ICONTXT (input) INTEGER * ICONTXT is the BLACS mechanism for partitioning communication * space. A defining property of a context is that a message in * a context cannot be sent or received in another context. The * BLACS context includes the definition of a grid, and each * process' coordinates in it. * * TRANS (input) CHARACTER*1 * TRANS specifies the operation to be performed as follows: * * TRANS = 'N', y := alpha*A*x + beta*y. * TRANS = 'T', y := alpha*A'*x + beta*y. * TRANS = 'C', y := alpha*A'*x + beta*y. * * XDIST (input) CHARACTER*1 * XDIST specifies the distribution of vector x as follows: * * XDIST = 'C', x is distributed columnwise * or in a column of processes * XDIST = 'R', x is distributed rowwise * or in a row of processes * * YDIST (input) CHARACTER*1 * YDIST specifies the distribution of vector y as follows: * * YDIST = 'C', y is distributed columnwise * or in a column of processes * YDIST = 'R', y is distributed rowwise * or in a row of processes * * M (input) INTEGER * M specifies the (global) number of rows of the matrix A. * M >= 0. * * N (input) INTEGER * N specifies the (global) number of columns of the matrix A. * N >= 0. * * MB (input) INTEGER * MB specifies the row block size of the matrix A. It also * specifies the block size of the vector Y if TRANS = 'N', or * the vector X if TRANS = 'T'/'C'. MB >= 1. * * NB (input) INTEGER * NB specifies the column block size of the matrix A. It also * specifies the block size of the vector X if TRANS = 'N', or * the vector Y if TRANS = 'T'/'C'. NB >= 1. * * MZ (input) INTEGER * MZ is the row offset to specify the row distance from the * beginning of the block to the first element of A. * It also specifies the offset to the first element of Y if * TRANS = 'N', or to the first element of X if TRANS = 'T'/'C'. * 0 <= MZ < MB. * * NZ (input) INTEGER * NZ is the column offset to specify the column distance from * the beginning of the block to the first element of A. It * also specifies the offset to the first element of X if TRANS * = 'N', or to the first element of Y if TRANS = 'T'/'C'. * 0 <= NZ < NB. * * ALPHA (input) DOUBLE PRECISION * ALPHA specifies the scalar alpha. * * A (input) DOUBLE PRECISION array of DIMENSION ( LDA, Nq ), * The leading Mp-by-Nq part of the array A must contain the * matrix A. * * LDA (input) INTEGER * LDA specifies the first dimension of (local) A as declared * in the calling (sub) program. LDA >= MAX(1,Mp). * * X (input) DOUBLE PRECISION array of DIMENSION at least * (1+(Np-1)*abs(INCX)) if TRANS = 'N' and XDIST = 'C', * (1+(Nq-1)*abs(INCX)) if TRANS = 'N' and XDIST = 'R', * (1+(Mp-1)*abs(INCX)) if TRANS = 'T'/'C' and XDIST = 'C', * (1+(Mq-1)*abs(INCX)) if TRANS = 'T'/'C' and XDIST = 'R'. * The incremented array X must contain the vector X. * * INCX (input) INTEGER * INCX specifies the increment for the elements of X. * INCX <> 0. * * BETA (input) DOUBLE PRECISION * BETA specifies the scalar beta. When BETA is supplied as * zero then Y need not be set on input. * * Y (input/output) DOUBLE PRECISION array of DIMENSION at least * (1+(Mp-1)*abs(INCY)) if TRANS = 'N' and YDIST = 'C', * (1+(Mq-1)*abs(INCY)) if TRANS = 'N' and YDIST = 'R', * (1+(Np-1)*abs(INCY)) if TRANS = 'T'/'C' and YDIST = 'C', * (1+(Nq-1)*abs(INCY)) if TRANS = 'T'/'C' and YDIST = 'R'. * On entry with BETA non-zero, the incremented array Y must * contain the vector Y. * On exit, Y is overwritten by the updated vector Y. * * INCY (input) INTEGER * INCY specifies the increment for the elements of Y. * INCY <> 0. * * IAROW (input) INTEGER * IAROW specifies a row of the process template, which holds * the first block of the matrix A. * * IACOL (input) INTEGER * IACOL specifies a column of the process template, which * holds the first block of the matrix A. * * IXROW (input) INTEGER * IXROW specifies the current row of process template * which has the first element of X. If all row processes * have their own copies of X, which is a row vector, * set IXROW = -1. * * IXCOL (input) INTEGER * IXCOL specifies the current column of process template * which has the first element of X. If all column processes * have their own copies of X, which is a column vector, set * IXCOL = -1. * * IYROW (input) INTEGER * IYROW specifies the current row of process template which * has the first element of Y. * * IYCOL (input) INTEGER * IYCOL specifies the current column of process template * which has the first element of Y. * * XCOMM (input) CHARACTER*1 * XCOMM specifies the communication scheme of row or column of * X. It follows topology definition of BLACS. If vector * transpose of X is involved, the value is ignored and * it is set to '1-tree'. * * XWORK (input) CHARACTER*1 * XWORK determines whether X is a workspace or not. * * XWORK = 'Y': X is workspace in other processes. * X is sent to X position in other processes. * It is assumed that processes have * sufficient space to store (local) X. * XWORK = 'N': Data in X will be untouched (unchanged). * * If transposition of X is involved with the computation of Y, * the argument is ignored. * * YWORK (input) CHARACTER*1 * YWORK determines whether Y is a workspace or not. * * YWORK = 'Y': Y is workspace in other processes. * It is assumed that processes have * sufficient space to store temporary * (local) Y. * YWORK = 'N': Data in Y will be untouched (unchanged) * in other processes. * * If transposition of Y is required during computation, * the argument is ignored. * * WORK (workspace) DOUBLE PRECISION array of dimension Size(WORK). * It will store copy of X and/or copy of Y. (see requirements) * * Parameters Details * ================== * * Lx It is a local portion of L (L is replaced by either M or N, * and x is replaced by either p (=NPROW) or q (=NPCOL)), owned * by a process. The value is determined by L, LB, LZ, x, and * MI, where LB is a block size, LZ is a offset from the * beginning of the block, MI is a row or column position in * process template. Lx is equal to or less than CEIL( L+LZ, * LB*x ) * LB. * * Memory Requirement of WORK * ========================== * * MM = M + MZ * NN = N + NZ * Mpb = CEIL( MM, MB*NPROW ) * Nqb = CEIL( NN, NB*NPCOL ) * Mp0 = NUMROC( MM, MB, 0, 0, NPROW ) ~= Mpb * MB * Nq0 = NUMROC( NN, NB, 0, 0, NPCOL ) ~= Nqb * NB * LCMQ = LCM / NPCOL * LCMP = LCM / NPROW * * (1) TRANS = 'N' * (i) XDIST = 'Col' & YDIST = 'Col' * Size(WORK) = Nq0 * + MAX[ Mp0 (if YWORK <> 'Y'), * CEIL(Nqb,LCMQ)*NB (if IXCOL <> -1), * CEIL(Nqb,LCMQ)*NB*MIN(LCMQ,CEIL(NN,NB)) * (if IXCOL = -1) ] * (ii) XDIST = 'Col' & YDIST = 'Row' * Size(WORK) = Mp0 * + MAX[ CEIL(Mpb,LCMP)*MB, * Nq0 + MAX[ CEIL(Nqb,LCMQ)*NB (if IXCOL <> -1), * CEIL(Nqb,LCMQ)*NB*MIN(LCMQ,CEIL(NN,NB)) * (if IXCOL = -1) ] ] * * (iii) XDIST = 'Row' & YDIST = 'Col' * Size(WORK) = Nq0 (if XWORK <> 'Y' & IXCOL <> -1) * + Mp0 (if YWORK <> 'Y') * * (iv) XDIST = 'Row' & YDIST = 'Row' * Size(WORK) = Mp0 * + MAX[ Nq0 (if XWORK <> 'Y' & IXCOL <> -1), * CEIL(Mqb,LCMQ)*MB ] * * (2) TRANS = 'T'/'C' * (i) XDIST = 'Col' & YDIST = 'Col' * Size(WORK) = Nq0 * + MAX[ CEIL(Npb,LCMP)*NB, * Mp0 (if XWORK <> 'Y' & IXCOL <> -1) ] * * (ii) XDIST = 'Col' & YDIST = 'Row' * Size(WORK) = Mp0 (if XWORK <> 'Y' & IXCOL <> -1) * + Nq0 (if YWORK <> 'Y') * * (iii) XDIST = 'Row' & YDIST = 'Col' * Size(WORK) = Nq0 * + MAX[ CEIL(Nqb,LCMQ)*NB, * Mp0 + MAX[ CEIL(Mpb,LCMP)*MB (if IXROW <> -1), * CEIL(Mpb,LCMP)*MB*MIN(LCMP,CEIL(MM,MB)) * (if IXROW = -1) ] ] * * (iv) XDIST = 'Row' & YDIST = 'Row' * Size(WORK) = Mp0 * + MAX[ Nq0 (if YWORK <> 'Y'), * CEIL(Mpb,LCMP)*MB (if IXROW <> -1), * CEIL(Mpb,LCMP)*MB*MIN(LCMP,CEIL(MM,MB)) * (if IXROW = -1) ] * * Notes * ----- * More precise space can be computed as * * CEIL(Mpb,LCMP)*MB => NUMROC( NUMROC(MM,MB,0,0,NPROW), MB, 0, 0, LCMP) * = NUMROC( Mp0, MB, 0, 0, LCMP ) * CEIL(Nqb,LCMQ)*NB => NUMROC( NUMROC(NN,NB,0,0,NPCOL), NB, 0, 0, LCMQ) * = NUMROC( Nq0, NB, 0, 0, LCMQ ) * * ===================================================================== * * .. * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. CHARACTER*1 COMMX LOGICAL NOTRAN, XDATA, XCOL, YCOL INTEGER INFO, IPX, IPY, MP, MYCOL, MYROW, NPCOL, NPROW, $ NQ DOUBLE PRECISION DUMMY, TBETA * .. * .. External Functions .. LOGICAL LSAME INTEGER NUMROC EXTERNAL LSAME, NUMROC * .. * .. External Subroutines .. EXTERNAL BLACS_GRIDINFO, DGEBR2D, DGEBS2D, DGEMV, $ DGSUM2D, PBDTRNV, PBDVECADD, PXERBLA * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * * Quick return if possible. * IF( M.EQ.0 .OR. N.EQ.0 .OR. ( ALPHA.EQ.ZERO .AND. BETA.EQ.ONE ) ) $ RETURN * CALL BLACS_GRIDINFO( ICONTXT, NPROW, NPCOL, MYROW, MYCOL ) * NOTRAN = LSAME( TRANS, 'N' ) XCOL = LSAME( XDIST, 'C' ) YCOL = LSAME( YDIST, 'C' ) * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTRAN ).AND. $ ( .NOT.LSAME( TRANS, 'C' ) ).AND. $ ( .NOT.LSAME( TRANS, 'T' ) ) ) THEN INFO = 2 ELSE IF( ( .NOT.XCOL ).AND. $ ( .NOT.LSAME( XDIST, 'R' ) ) ) THEN INFO = 3 ELSE IF( ( .NOT.YCOL ).AND. $ ( .NOT.LSAME( YDIST, 'R' ) ) ) THEN INFO = 4 ELSE IF( M .LT.0 ) THEN INFO = 5 ELSE IF( N .LT.0 ) THEN INFO = 6 ELSE IF( MB .LT.1 ) THEN INFO = 7 ELSE IF( NB .LT.1 ) THEN INFO = 8 ELSE IF( MZ .LT.0 .OR. MZ .GE. MB ) THEN INFO = 9 ELSE IF( NZ .LT.0 .OR. NZ .GE. NB ) THEN INFO = 10 ELSE IF( INCX.EQ.0 ) THEN INFO = 15 ELSE IF( INCY.EQ.0 ) THEN INFO = 18 ELSE IF( IAROW.LT.0 .OR. IAROW.GE.NPROW ) THEN INFO = 19 ELSE IF( IACOL.LT.0 .OR. IACOL.GE.NPCOL ) THEN INFO = 20 END IF * 10 CONTINUE IF( INFO.NE.0 ) THEN CALL PXERBLA( ICONTXT, 'PBDGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Initialize parameters * MP = NUMROC( M+MZ, MB, MYROW, IAROW, NPROW ) IF( MYROW.EQ.IAROW ) MP = MP - MZ NQ = NUMROC( N+NZ, NB, MYCOL, IACOL, NPCOL ) IF( MYCOL.EQ.IACOL ) NQ = NQ - NZ IF( LDA.LT.MAX(1,MP) ) INFO = 13 COMMX = XCOMM IF( LSAME( COMMX, ' ' ) ) COMMX = '1' * * Start the operations. * * If A is not transposed, * IF( NOTRAN ) THEN IF( XCOL ) THEN IF( YCOL ) THEN * * Form y := alpha*A*x + beta*y, * where x and y are distributed columnwise. * __________ * || | | || * || | | || || * || | | || || * (y) = alpha * | A | * (x) + beta *(y) * || | | || || * || | | || || * || |__________| || * IF( IXROW.LT.0 .OR. IXROW.GE.NPROW ) THEN INFO = 21 ELSE IF( IXCOL.LT.-1 .OR. IXCOL.GE.NPCOL ) THEN INFO = 22 ELSE IF( IYROW.NE.IAROW ) THEN INFO = 23 ELSE IF( IYCOL.LT.0 .OR. IYCOL.GE.NPCOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * * Transpose X to WORK * IPY = NQ + 1 CALL PBDTRNV( ICONTXT, 'Col', 'T', N, NB, NZ, X, INCX, ZERO, $ WORK, 1, IXROW, IXCOL, -1, IACOL, WORK(IPY) ) * * Compute Y if Y is distributed columnwise * IF( LSAME( YWORK, 'Y' ) ) THEN TBETA = ZERO IF( MYCOL.EQ.IYCOL ) TBETA = BETA * IF( NQ.GT.0 ) THEN CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, WORK, 1, TBETA, $ Y, INCY ) ELSE IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL $ .OR. MYCOL.EQ.IYCOL ) THEN CALL PBDVECADD( ICONTXT, 'V', MP, ZERO, DUMMY, 1, TBETA, $ Y, INCY ) END IF * * Add Y rowwise * IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, Y, INCY, $ MYROW, IYCOL ) * ELSE IF( NQ.GT.0 ) THEN CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, WORK, 1, ZERO, $ WORK(IPY), 1 ) ELSE CALL PBDVECADD( ICONTXT, 'G', MP, ZERO, DUMMY, 1, ZERO, $ WORK(IPY), 1 ) END IF * * Add WORK(IPY) rowwise * IF( MYCOL.EQ.IYCOL ) THEN CALL PBDVECADD( ICONTXT, 'G', MP, ONE, WORK(IPY), 1, $ BETA, Y, INCY ) IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, Y, $ INCY, MYROW, IYCOL ) ELSE IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, $ WORK(IPY), 1, MYROW, IYCOL ) END IF END IF * ELSE * * Form y := alpha*A*x + beta*y, * where x is distributed columnwise & y is distributed rowwise * __________ * | | * | | || * | | || * =====(y)===== = a * | A | * (x) + b * =====(y)===== * | | || * | | || * |__________| * IF( IXROW.LT.0 .OR. IXROW.GE.NPROW ) THEN INFO = 21 ELSE IF( IXCOL.LT.-1 .OR. IXCOL.GE.NPCOL ) THEN INFO = 22 ELSE IF( IYROW.LT.0 .OR. IYROW.GE.NPROW ) THEN INFO = 23 ELSE IF( IYCOL.LT.0 .OR. IYCOL.GE.NPCOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * * Transpose X to WORK * IPX = MP + 1 CALL PBDTRNV( ICONTXT, 'Col', 'T', N, NB, NZ, X, INCX, ZERO, $ WORK(IPX), 1, IXROW, IXCOL, -1, IACOL, $ WORK(NQ+IPX) ) * IF( NQ.GT.0 ) THEN CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, WORK(IPX), 1, $ ZERO, WORK, 1 ) ELSE CALL PBDVECADD( ICONTXT, 'G', MP, ZERO, DUMMY, 1, ZERO, $ WORK, 1 ) END IF * * Add Y (=WORK) rowwise and transpose it. * IF( N+NZ.GT.NB ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, WORK, 1, $ MYROW, IACOL ) CALL PBDTRNV( ICONTXT, 'Col', 'T', M, MB, MZ, WORK, 1, BETA, $ Y, INCY, IAROW, IACOL, IYROW, IYCOL, $ WORK(IPX) ) END IF * ELSE IF( YCOL ) THEN * * Form y := alpha*A*x + beta*y, * where x is distributed rowwise & y is distributed columnwise * __________ * || | | || * || | | || * || | | || * (y) = alpha*| A |* ====(x)=== + beta*(y) * || | | || * || | | || * || |__________| || * IF( IXROW.LT.-1 .OR. IXROW.GE.NPROW ) THEN INFO = 21 ELSE IF( IXCOL.NE.IACOL ) THEN INFO = 22 ELSE IF( IYROW.NE.IAROW ) THEN INFO = 23 ELSE IF( IYCOL.LT.0 .OR. IYCOL.GE.NPCOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * XDATA = .FALSE. IF( IXROW.EQ.-1 ) XDATA = .TRUE. IPY = 1 * * Broadcast X to X or WORK if necessary ( IXCOL <> -1 ) * IF( .NOT.XDATA ) THEN IF( LSAME( XWORK, 'Y' ) ) THEN IF( MYROW.EQ.IXROW ) THEN CALL DGEBS2D( ICONTXT, 'Col', COMMX, 1, NQ, X, INCX ) ELSE CALL DGEBR2D( ICONTXT, 'Col', COMMX, 1, NQ, X, INCX, $ IXROW, MYCOL ) END IF XDATA = .TRUE. ELSE IF( MYROW.EQ.IXROW ) THEN CALL PBDVECADD( ICONTXT, 'V', NQ, ONE, X, INCX, ZERO, $ WORK, 1 ) CALL DGEBS2D( ICONTXT, 'Col', COMMX, 1, NQ, WORK, 1 ) ELSE CALL DGEBR2D( ICONTXT, 'Col', COMMX, 1, NQ, WORK, 1, $ IXROW, MYCOL ) END IF IPY = NQ + 1 END IF END IF * * Compute Y * IF( LSAME( YWORK, 'Y' ) ) THEN TBETA = ZERO IF( MYCOL.EQ.IYCOL ) TBETA = BETA * IF( NQ.GT.0 ) THEN IF( XDATA ) THEN CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, X, INCX, $ TBETA, Y, INCY ) ELSE CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, WORK, 1, $ TBETA, Y, INCY ) END IF ELSE IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL $ .OR. MYCOL.EQ.IYCOL ) THEN CALL PBDVECADD( ICONTXT, 'V', MP, ZERO, DUMMY, 1, TBETA, $ Y, INCY ) END IF * * Add Y rowwise * IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, Y, INCY, $ MYROW, IYCOL ) ELSE IF( NQ.GT.0 ) THEN IF( XDATA ) THEN CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, X, INCX, $ ZERO, WORK(IPY), 1 ) ELSE CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, WORK, 1, $ ZERO, WORK(IPY), 1 ) END IF ELSE CALL PBDVECADD( ICONTXT, 'G', MP, ZERO, DUMMY, 1, ZERO, $ WORK(IPY), 1 ) END IF * * Add Y rowwise * IF( MYCOL.EQ.IYCOL ) THEN CALL PBDVECADD( ICONTXT, 'G', MP, ONE, WORK(IPY), 1, $ BETA, Y, INCY ) IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, Y, $ INCY, MYROW, IYCOL ) ELSE IF( N+NZ.GT.NB .OR. IACOL.NE.IYCOL ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, $ WORK(IPY), 1, MYROW, IYCOL ) END IF END IF * ELSE * * Form y := alpha*A*x + beta*y, * where x and y are distributed rowwise. * __________ * | | * | | * | | * ======(y)===== = a*| A |* ====(x)=== + b* ======(y)===== * | | * | | * |__________| * IF( IXROW.LT.0 .OR. IXROW.GE.NPROW ) THEN INFO = 21 ELSE IF( IXCOL.NE.IACOL ) THEN INFO = 22 ELSE IF( IYROW.LT.0 .OR. IYROW.GE.NPROW ) THEN INFO = 23 ELSE IF( IYCOL.LT.0 .OR. IYCOL.GE.NPCOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * XDATA = .FALSE. IF( IXROW.EQ.-1 ) XDATA = .TRUE. IPX = MP + 1 * * Broadcast X to X or WORK(IPX) if necessary ( IXCOL <> -1 ) * IF( .NOT.XDATA ) THEN IF( LSAME( XWORK, 'Y' ) ) THEN IF( MYROW.EQ.IXROW ) THEN CALL DGEBS2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ X, INCX ) ELSE CALL DGEBR2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ X, INCX, IXROW, MYCOL ) END IF XDATA = .TRUE. ELSE IF( MYROW.EQ.IXROW ) THEN CALL PBDVECADD( ICONTXT, 'V', NQ, ONE, X, INCX, ZERO, $ WORK(IPX), 1 ) CALL DGEBS2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ WORK(IPX), 1 ) ELSE CALL DGEBR2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ WORK(IPX), 1, IXROW, MYCOL ) END IF END IF END IF * * Compute Y * IF( NQ.GT.0 ) THEN IF( XDATA ) THEN CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, X, INCX, $ ZERO, WORK, 1 ) ELSE CALL DGEMV( 'No', MP, NQ, ALPHA, A, LDA, WORK(IPX), 1, $ ZERO, WORK, 1 ) END IF ELSE CALL PBDVECADD( ICONTXT, 'G', MP, ZERO, DUMMY, 1, ZERO, $ WORK, 1 ) END IF * * Add Y rowwise and transpose it. * IF( N+NZ.GT.NB ) $ CALL DGSUM2D( ICONTXT, 'Row', '1-tree', 1, MP, WORK, 1, $ MYROW, IACOL ) CALL PBDTRNV( ICONTXT, 'Col', 'T', M, MB, MZ, WORK, 1, BETA, $ Y, INCY, IAROW, IACOL, IYROW, IYCOL, $ WORK(IPX) ) END IF END IF * ELSE IF( XCOL ) THEN IF( YCOL ) THEN * * Form Y := alpha*(A')*X + beta*Y. * where X and Y are distributed columnwise. * __________ * | | || * || | | || || * || | | || || * (y) = alpha * | (A') | * (x) + beta *(y) * || | | || || * || | | || || * |__________| || * * IF( IXROW.NE.IAROW ) THEN INFO = 21 ELSE IF( IXCOL.LT.-1 .OR. IXCOL.GE.NPCOL ) THEN INFO = 22 ELSE IF( IYROW.LT.0 .OR. IYROW.GE.NPROW ) THEN INFO = 23 ELSE IF( IYCOL.LT.0 .OR. IYCOL.GE.NPCOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * XDATA = .FALSE. IF( IXCOL.EQ.-1 ) XDATA = .TRUE. IPX = NQ + 1 * * Broadcast X to X or WORK(IPX) if necessary ( IXCOL <> -1 ) * IF( .NOT.XDATA ) THEN IF( LSAME( XWORK, 'Y' ) ) THEN IF( MYCOL.EQ.IXCOL ) THEN CALL DGEBS2D( ICONTXT, 'Row', '1-tree', 1, MP, $ X, INCX ) ELSE CALL DGEBR2D( ICONTXT, 'Row', '1-tree', 1, MP, $ X, INCX, MYROW, IXCOL ) END IF XDATA = .TRUE. ELSE IF( MYCOL.EQ.IXCOL ) THEN CALL PBDVECADD( ICONTXT, 'V', MP, ONE, X, INCX, ZERO, $ WORK(IPX), 1 ) CALL DGEBS2D( ICONTXT, 'Row', '1-tree', 1, MP, $ WORK(IPX), 1 ) ELSE CALL DGEBR2D( ICONTXT, 'Row', '1-tree', 1, MP, $ WORK(IPX), 1, MYROW, IXCOL ) END IF END IF END IF * * Compute Y' (<= X' * A ) * IF( MP.GT.0 ) THEN IF( XDATA ) THEN CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, X, INCX, ZERO, $ WORK, 1 ) ELSE CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, WORK(IPX), 1, $ ZERO, WORK, 1 ) END IF ELSE CALL PBDVECADD( ICONTXT, 'G', NQ, ZERO, DUMMY, 1, ZERO, $ WORK, 1 ) END IF * * Transpose Y (= WORK) and transpose it. * IF( M+MZ.GT.MB ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, WORK, 1, $ IAROW, MYCOL ) CALL PBDTRNV( ICONTXT, 'Row', 'T', N, NB, NZ, WORK, 1, BETA, $ Y, INCY, IAROW, IACOL, IYROW, IYCOL, $ WORK(IPX) ) * ELSE * * Form y := alpha*(A')*x + beta*y. * where x is distributed columnwise & y is distributed rowwise * __________ * | | || * | | || * | | || * ====(y)==== = alpha * | (A') | * (x) + beta * ====(y)==== * | | || * | | || * |__________| || * * IF( IXROW.LT.IAROW ) THEN INFO = 21 ELSE IF( IXCOL.LT.-1 .OR. IXCOL.GE.NPCOL ) THEN INFO = 22 ELSE IF( IYROW.LT.0 .OR. IYROW.GE.NPROW ) THEN INFO = 23 ELSE IF( IYCOL.NE.IACOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * XDATA = .FALSE. IF( IXCOL.EQ.-1 ) XDATA = .TRUE. IPY = 1 * * Broadcast X to X or WORK(IPX) if necessary ( IXCOL <> -1 ) * IF( .NOT.XDATA ) THEN IF( LSAME( XWORK, 'Y' ) ) THEN IF( MYCOL.EQ.IXCOL ) THEN CALL DGEBS2D( ICONTXT, 'Row', COMMX, 1, MP, X, INCX ) ELSE CALL DGEBR2D( ICONTXT, 'Row', COMMX, 1, MP, X, INCX, $ MYROW, IXCOL ) END IF XDATA = .TRUE. ELSE IF( MYCOL.EQ.IXCOL ) THEN CALL PBDVECADD( ICONTXT, 'V', MP, ONE, X, INCX, ZERO, $ WORK, 1 ) CALL DGEBS2D( ICONTXT, 'Row', COMMX, 1, MP, WORK, 1 ) ELSE CALL DGEBR2D( ICONTXT, 'Row', COMMX, 1, MP, WORK, 1, $ MYROW, IXCOL ) END IF IPY = MP + 1 END IF END IF * * Compute Y' (<= X' * A ) * IF( LSAME( YWORK, 'Y' ) ) THEN TBETA = ZERO IF( MYROW.EQ.IYROW ) TBETA = BETA * IF( MP.GT.0 ) THEN IF( XDATA ) THEN CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, X, INCX, $ TBETA, Y, INCY ) ELSE CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, WORK, 1, $ TBETA, Y, INCY ) END IF ELSE IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW .OR. $ MYROW.EQ.IYROW ) THEN CALL PBDVECADD( ICONTXT, 'V', NQ, ZERO, DUMMY, 1, TBETA, $ Y, INCY ) END IF * * Add Y columnwise * IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, Y, INCY, $ IYROW, MYCOL ) * ELSE IF( MP.GT.0 ) THEN IF( XDATA ) THEN CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, X, INCX, $ ZERO, WORK(IPY), 1 ) ELSE CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, WORK, 1, $ ZERO, WORK(IPY), 1 ) END IF ELSE CALL PBDVECADD( ICONTXT, 'G', NQ, ZERO, DUMMY, 1, ZERO, $ WORK(IPY), 1 ) END IF * * Add Y (= WORK(IPY)) columnwise * IF( MYROW.EQ.IYROW ) THEN CALL PBDVECADD( ICONTXT, 'G', NQ, ONE, WORK(IPY), 1, $ BETA, Y, INCY ) IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ Y, INCY, IYROW, MYCOL ) ELSE IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ WORK(IPY), 1, IYROW, MYCOL ) END IF END IF END IF * ELSE IF( YCOL ) THEN * * Form y := alpha*A*x + beta*y, * where x is distributed rowwise & y is distributed columnwise * __________ * | | * || | | || * || | | || * (y) = alpha*| (A') |* =====(x)===== + beta*(y) * || | | || * || | | || * |__________| * IF( IXROW.LT.-1 .OR. IXROW.GE.NPROW ) THEN INFO = 21 ELSE IF( IXCOL.LT.0 .OR. IXCOL.GE.NPCOL ) THEN INFO = 22 ELSE IF( IYROW.LT.0 .OR. IYROW.GE.NPROW ) THEN INFO = 23 ELSE IF( IYCOL.LT.0 .OR. IYCOL.GE.NPCOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * * Transpose X to WORK * IPX = NQ + 1 CALL PBDTRNV( ICONTXT, 'Row', 'T', M, MB, MZ, X, INCX, ZERO, $ WORK(IPX), 1, IXROW, IXCOL, IAROW, -1, $ WORK(MP+IPX) ) * * Compute Y * IF( MP.GT.0 ) THEN CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, WORK(IPX), 1, $ ZERO, WORK, 1 ) ELSE CALL PBDVECADD( ICONTXT, 'G', NQ, ZERO, DUMMY, 1, ZERO, $ WORK, 1 ) END IF * * Add Y columnwise and transpose it. * IF( M+MZ.GT.MB ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, WORK, 1, $ IAROW, MYCOL ) CALL PBDTRNV( ICONTXT, 'Row', 'T', N, NB, NZ, WORK, 1, BETA, $ Y, INCY, IAROW, IACOL, IYROW, IYCOL, $ WORK(IPX) ) * ELSE * * Form y := alpha*(A')*x + beta*y. * where x and y are distributed rowwise. * __________ * | | * | | * | | * =====(y)==== = a*| (A') |* ======(x)===== + b* =====(y)==== * | | * | | * |__________| * IF( IXROW.LT.-1 .OR. IXROW.GE.NPROW ) THEN INFO = 21 ELSE IF( IXCOL.LT.0 .OR. IXCOL.GE.NPCOL ) THEN INFO = 22 ELSE IF( IYROW.LT.0 .OR. IYROW.GE.NPROW ) THEN INFO = 23 ELSE IF( IYCOL.NE.IACOL ) THEN INFO = 24 END IF IF( INFO.NE.0 ) GO TO 10 * * Transpose X to WORK * IPY = MP + 1 CALL PBDTRNV( ICONTXT, 'Row', 'T', M, MB, MZ, X, INCX, ZERO, $ WORK, 1, IXROW, IXCOL, IAROW, -1, WORK(IPY) ) * * Compute Y * IF( LSAME( YWORK, 'Y' ) ) THEN TBETA = ZERO IF( MYROW.EQ.IYROW ) TBETA = BETA * IF( MP.GT.0 ) THEN CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, WORK, 1, $ TBETA, Y, INCY ) ELSE IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW $ .OR. MYROW.EQ.IYROW ) THEN CALL PBDVECADD( ICONTXT, 'V', NQ, ZERO, DUMMY, 1, TBETA, $ Y, INCY ) END IF * * Add Y columnwise if necessary * IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, Y, INCY, $ IYROW, MYCOL ) * ELSE IF( MP.GT.0 ) THEN CALL DGEMV( TRANS, MP, NQ, ALPHA, A, LDA, WORK, 1, ZERO, $ WORK(IPY), 1 ) ELSE CALL PBDVECADD( ICONTXT, 'G', NQ, ZERO, DUMMY, 1, ZERO, $ WORK(IPY), 1 ) END IF * * Add Y columnwise if necessary * IF( MYROW.EQ.IYROW ) THEN CALL PBDVECADD( ICONTXT, 'G', NQ, ONE, WORK(IPY), 1, $ BETA, Y, INCY ) IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ Y, INCY, IYROW, MYCOL ) ELSE IF( M+MZ.GT.MB .OR. IAROW.NE.IYROW ) $ CALL DGSUM2D( ICONTXT, 'Col', '1-tree', 1, NQ, $ WORK(IPY), 1, IYROW, MYCOL ) END IF END IF END IF END IF END IF * RETURN * * End of PBDGEMV * END