SUBROUTINE PBDTRSM( ICONTXT, MATBLK, SIDE, UPLO, TRANSA, DIAG, M,
     $                    N, NB, ALPHA, A, LDA, B, LDB, IAROW, IACOL,
     $                    IBPOS, ACOMM, ABWORK, 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        ABWORK, ACOMM, DIAG, MATBLK, SIDE, TRANSA,
     $                   UPLO
      INTEGER            IACOL, IAROW, IBPOS, ICONTXT, LDA, LDB, M, N,
     $                   NB
      DOUBLE PRECISION   ALPHA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PBDTRSM is a parallel blocked version of he Level 3 BLAS routine
*  DTRSM.
*  PBDTRSM solves one of the matrix equations  based on block
*  cyclic distribution.
*
*     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
*
*  where alpha is a scalar, X and B are m-by-n matrices, A is a unit, or
*  non-unit,  upper or lower triangular matrix.  op( A ) is one of
*
*     op( A ) = A,  A**T,  or  A**H
*
*  where the size of the matrix op( A ) is M-by-M  if SIDE = 'L', and N-
*  by-N  otherwise.  The M-by-N  matrix B  is a column block  (only one
*  column of processes have B) if SIDE = 'L', and a row block otherwise
*  (only one row of processes have B).  The matrix X is overwritten  on
*  B.
*
*  The first elements  of the matrices A, and B  should  be  located  at
*  the beginnings of their first blocks. (not the middle of the blocks.)
*  When MATBLK = 'M',  B can be moved  or  transposed  to  the  starting
*  column or row processes if necessary.  The communication scheme is
*  predetermined.
*  And when MATBLK = 'B',  A can be  broadcast columnwise or rowwise  if
*  necessary.  The communication scheme can be selected.
*
*  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.
*
*  MATBLK  (input) CHARACTER*1
*          MATBLK specifies whether op( A ) is a (full) block matrix or
*          a single block as follows:
*
*             MATBLK = 'M',  op( A ) is a (full) block matrix,
*             MATBLK = 'B',  op( A ) is a single block.
*
*  SIDE    (input) CHARACTER*1
*          SIDE specifies whether op( A ) appears on the left or right
*          of X as follows:
*
*             SIDE = 'L',  op( A )*X = alpha*B,
*             SIDE = 'R',  X*op( A ) = alpha*B.
*
*  UPLO    (input) CHARACTER*1
*          UPLO specifies whether the matrix A is an upper or lower
*          triangular matrix as follows:
*
*             UPLO = 'U',  A is an upper triangular matrix.
*             UPLO = 'L',  A is a lower triangular matrix.
*
*  TRANSA  (input) CHARACTER*1
*          TRANSA specifies the form of op( A ) to be used in
*          the matrix multiplication as follows:
*
*             TRANSA = 'N',  op( A ) = A.
*             TRANSA = 'T',  op( A ) = A**T.
*             TRANSA = 'C',  op( A ) = A**H.
*
*  DIAG    (input) CHARACTER*1
*          DIAG specifies whether or not A is unit triangular as
*          follows:
*
*             DIAG = 'U'   A is assumed to be unit triangular.
*             DIAG = 'N'   A is not assumed to be unit
*                                 triangular.
*
*  M       (input) INTEGER
*          M specifies the number of rows of B.  M >= 0.
*
*  N       (input) INTEGER
*          N specifies the number of columns of B.  N >= 0.
*
*  NB      (input) INTEGER
*          NB specifies the row and column block size of matrix A.
*          It also specifies the row block size of the matrix B if
*          MATBLK = 'M' and SIDE = 'L', or MATBLK = 'B' and SIDE = 'R';
*          and the column block size of the matrix B if MATBLK = 'M'
*          and SIDE = 'R', or MATBLK = 'B' and SIDE = 'L'.  NB >= 1.
*
*  ALPHA   (input) DOUBLE PRECISION
*          ALPHA specifies the scalar  alpha.  When alpha is zero,
*          A is not referenced and  B need not be set before entry.
*
*  A       (input) DOUBLE PRECISION array of DIMENSION ( LDA, Kq ),
*          where kq is Mq (Kp is Mp) when SIDE = 'L' and is Nq (Kp is
*          Np) when SIDE = 'R'.
*          If SIDE = `L', the M-by-M part of the array A must contain
*          the (global) triangular matrix, such that when UPLO = 'U',
*          the leading M-by-M upper triangular part of the array A must
*          contain the upper triangular part of the (global) matrix and
*          the strictly  lower triangular part of A is not referenced,
*          and when  UPLO = 'L', the leading M-by-M lower triangular
*          part of the array A must  contain the lower triangular part
*          of the (global) matrix and the strictly upper triangular
*          part of A is not referenced.
*          And if SIDE = 'R', the N-by-N part of the (global) array A
*          must contain the (global) matrix, such that when UPLO = 'U',
*          the leading N-by-N upper triangular part of the array A must
*          contain the upper triangular part of the (global) matrix and
*          the strictly lower triangular part of A is not referenced,
*          and when UPLO = 'L', the leading N-by-N lower triangular
*          part of the array A must contain the lower triangular part
*          of the (global) matrix and the strictly upper triangular
*          part of A is not referenced.
*          Note that when DIAG = `U', the diagonal elements of A are
*          not referenced either, but are assumed to be unity.
*
*  LDA     (input) INTEGER
*          LDA specifies the first dimension of A as declared in the
*          calling (sub) program.  LDA >= MAX(1,Mp) if SIDE = 'L', and
*          LDA >= MAX(1,Np) otherwise.
*
*  B       (input/output) DOUBLE PRECISION array of DIMENSION (LDB,Nq)
*          On entry,  the leading Mp-by-Nq part of the array B must
*          contain the matrix B  when SIDE = 'R', or the leading Mp-by-
*          Nq part of the array B  must contain the (local) matrix B
*          otherwise.
*          On exit B is overwritten by the transformed matrix.  Input
*          values of B would be changed after the computation in the
*          processes which don't have the resultant column block or
*          row block of B if MATBLK = 'M'.
*
*  LDB     (input) INTEGER
*          LDB specifies the leading dimension of (local) B as declared
*          in the calling (sub) program.  LDB >= MAX(1,Mp).
*
*  IAROW   (input) INTEGER
*          It specifies a row of process template which has the
*          first block of A.  When MATBLK = 'B', and all rows of
*          processes have their own copies of A, set IAROW =  -1.
*
*  IACOL   (input) INTEGER
*          It specifies a column of process template which has the
*          first block of A.  When MATBLK = 'B', and all columns of
*          processes have their own copies of A, set IACOL = -1.
*
*  IBPOS   (input) INTEGER
*          When MATBLK = 'M', if SIDE = 'L', IBPOS specifies a column of
*          the process template, which holds the column of blocks of B
*          (0 <= IBPOS < NPCOL).  And if SIDE = 'R', it specifies a row
*          of the template, which holds the row of blocks of B (0 <=
*          IBPOS < NPROW).
*          When MATBLK = 'B', if SIDE = 'L', it specifies a column of
*          the template which has the first block of B (0 <= IBPOS
*          < NPCOL), and if SIDE = 'R', it specifies a row of the
*          template, which has the first block of B (0 <=IBPOS <NPROW).
*
*  ACOMM   (input) CHARACTER*1
*          When MATBLK = 'B', ACOMM specifies the communication scheme
*          of a block of A.  It follows topology definition of BLACS.
*          When MATBLK = 'M', the argument is ignored.
*
*  ABWORK  (input) CHARACTER*1
*          When MATBLK = 'M', ABWORK determines whether B is a
*          workspace or not.  If transposition of B is involved with
*          the computation, the argument is ignored.
*
*             ABWORK = 'Y':  B is workspace in other processes.
*                            B is overwitten with temporal B in other
*                            processes. It is assumed that processes
*                            have sufficient space to store temporal
*                            (local) B.
*             ABWORK = 'N':  Data of B in other processes will be
*                            untouched (unchanged).
*
*          And MATBLK = 'B', ABWORK determines whether A is a
*          workspace or not.
*
*             ABWORK = 'Y':  A is workspace in other processes.
*                            A is sent to A position in other processes.
*                            It is assumed that processes have
*                            sufficient space to store a single block A.
*             ABWORK = 'N':  Data of A in other processes will be
*                            untouched (unchanged).
*
*  WORK    (workspace) DOUBLE PRECISION array of dimension Size(WORK).
*          It will store copy of A or B if necessary.
*
*  Communication Scheme
*  ====================
*
*  If MATBLK='M', the communication scheme of the routine is determined
*  by the conditions  and it is independent of  machine characteristics,
*  so that it is not an option of the routine.  Increasing  ring  or
*  Decreasing ring is used depeding on the following input conditions.
*
*   COMM='Increasing ring' when UPLO = 'U', SIDE = 'L', TRANSA = 'T'/'C'
*                            or UPLO = 'U', SIDE = 'R', TRANSA = 'N'
*                            or UPLO = 'L', SIDE = 'L', TRANSA = 'N'
*                            or UPLO = 'L', SIDE = 'R', TRANSA = 'T'/'C'
*
*   COMM='Decreasing ring' when UPLO = 'U', SIDE = 'L', TRANSA = 'N'
*                            or UPLO = 'U', SIDE = 'R', TRANSA = 'T'/'C'
*                            or UPLO = 'L', SIDE = 'L', TRANSA = 'T'/'C'
*                            or UPLO = 'L', SIDE = 'R', TRANSA = 'N'
*
*  Parameters Details
*  ==================
*
*  Lx      It is  a local portion  of L  owned  by  a process,  (L is
*          replaced by M, or N,  and x  is replaced  by  either  p
*          (=NPROW) or q (=NPCOL)).  The value is determined by  L, LB,
*          x, and MI,  where  LB is  a block size  and MI is a  row  or
*          column position in a process template.  Lx is equal to  or
*          less than  Lx0 = CEIL( L, LB*x ) * LB.
*
*  Memory Requirement of WORK
*  ==========================
*
*  Mqb    = CEIL( M, NB*NPCOL )
*  Npb    = CEIL( N, NB*NPROW )
*  Mq0    = NUMROC( M, NB, 0, 0, NPCOL ) ~= Mqb * NB
*  Np0    = NUMROC( N, NB, 0, 0, NPROW ) ~= Npb * NB
*  LCMQ = LCM / NPCOL
*  LCMP = LCM / NPROW
*
*  (1) MATBLK = 'M'
*    (a) SIDE = 'Left'
*      (i)  TRANSA = 'N'
*          Size(WORK) = N*Mp0                       (if ABWORK <> 'Y')
*                     + N*NB*MAX[ 1, CEIL(Q-1,P) ]
*      (ii) TRANSA = 'T'/'C'
*          Size(WORK) = N*Mq0
*                     + MAX[ N*NB*MAX( CEIL(Mpb,LCMP), CEIL(Mqb,LCMQ) ),
*                            N*NB*CEIL(P-1,Q) ]
*
*    (b) SIDE = 'Right'
*      (i)  TRANSA = 'N'
*          Size(WORK) = M*Nq0                       (if ABWORK <> 'Y')
*                     + M*NB*MAX[ 1, CEIL(P-1,Q) ]
*      (ii) TRANSA = 'T'/'C'
*          Size(WORK) = M*Np0
*                     + MAX[ M*NB*MAX( CEIL(Npb,LCMP), CEIL(Nqb,LCMQ) ),
*                            M*NB*CEIL(Q-1,P) ]
*
*  (2) MATBLK = 'B'
*    (a) SIDE = 'Left'
*      Size(WORK) = M * M  (in IAROW; if IACOL <> -1 and ABWORK <> 'Y')
*    (b) SIDE = 'Right'
*      Size(WORK) = N * N  (in IACOL; if IAROW <> -1 and ABWORK <> 'Y')
*
*  Notes
*  -----
*  More precise space can be computed as
*
*  CEIL(Mqb,LCMQ)*NB => NUMROC( NUMROC(M,NB,0,0,NPCOL), NB, 0, 0, LCMQ )
*                    = NUMROC( Mq0, NB, 0, 0, LCMQ )
*  CEIL(Npb,LCMP)*NB => NUMROC( NUMROC(N,NB,0,0,NPROW), NB, 0, 0, LCMP )
*                    = NUMROC( Np0, NB, 0, 0, LCMP )
*
*  =====================================================================
*
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      CHARACTER*1        COMMA
      LOGICAL            ADATA, AMAT, ASPACE, BDATA, BSPACE, LSIDE,
     $                   NOTRAN, UPPER
      INTEGER            ICURCOL, ICURROW, IDEST, II, IIN, IN, INFO,
     $                   IPART, IPT, IRDB, IRPB, J, JB, JJ, JJN, JN, KB,
     $                   KDIST, LCM, LDW, MBTROW, MLFCOL, MP, MQ,
     $                   MRTCOL, MTPROW, MYCOL, MYROW, NCOMM, NDIM,
     $                   NLENG, NPART, NPCOL, NPROW, NQ, NREST, NS,
     $                   NXTCOL, NXTROW
      DOUBLE PRECISION   DUMMY
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL, ILCM, NUMROC
      EXTERNAL           ICEIL, ILCM, LSAME, NUMROC
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, DGEBR2D, DGEBS2D, DGEMM,
     $                   DGERV2D, DGESD2D, DTRBR2D, DTRBS2D, DTRSM,
     $                   PBDMATADD, PBDTRAN, PXERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN, MOD
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 ) RETURN
*
      CALL BLACS_GRIDINFO( ICONTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters.
*
      AMAT   = LSAME( MATBLK, 'M' )
      LSIDE  = LSAME( SIDE,   'L' )
      UPPER  = LSAME( UPLO,   'U' )
      NOTRAN = LSAME( TRANSA, 'N' )
*
      INFO = 0
      IF(      ( .NOT.AMAT                 ).AND.
     $         ( .NOT.LSAME( MATBLK, 'B' ) )           ) THEN
        INFO = 2
      ELSE IF( ( .NOT.LSIDE                ).AND.
     $         ( .NOT.LSAME( SIDE  , 'R' ) )           ) THEN
        INFO = 3
      ELSE IF( ( .NOT.UPPER                ).AND.
     $         ( .NOT.LSAME( UPLO  , 'L' ) )           ) THEN
        INFO = 4
      ELSE IF( ( .NOT.NOTRAN               ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) )           ) THEN
        INFO = 5
      ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
     $         ( .NOT.LSAME( DIAG  , 'N' ) )           ) THEN
        INFO = 6
      ELSE IF( M  .LT. 0                               ) THEN
        INFO = 7
      ELSE IF( N  .LT. 0                               ) THEN
        INFO = 8
      ELSE IF( NB .LT. 1                               ) THEN
        INFO = 9
      END IF
*
   10 CONTINUE
      IF( INFO.NE.0 ) THEN
        CALL PXERBLA( ICONTXT, 'PBDTRSM ', INFO )
        RETURN
      END IF
*
*     Start the operations.
*
* === If A is a matrix ===
*
      IF( AMAT ) THEN
*
*       Initialize parameters
*
        IF( LSIDE ) THEN
          NDIM = M
          NS   = N
        ELSE
          NDIM = N
          NS   = M
        END IF
        MP = NUMROC( NDIM, NB, MYROW, IAROW, NPROW )
        MQ = NUMROC( NDIM, NB, MYCOL, IACOL, NPCOL )
*
        IF( LDA.LT.MAX(1,MP)                    ) THEN
          INFO = 12
        ELSE IF( IAROW.LT.0 .OR. IAROW.GE.NPROW ) THEN
          INFO = 15
        ELSE IF( IACOL.LT.0 .OR. IACOL.GE.NPCOL ) THEN
          INFO = 16
        END IF
*
        BSPACE = LSAME( ABWORK, 'Y' )
        IF( LSIDE ) THEN
          IF( LDB.LT.MAX(1,MP) .AND. ( BSPACE .OR.
     $          IBPOS.EQ.MYCOL .OR. IBPOS.EQ.-1 )       ) THEN
            INFO = 14
          ELSE IF( IBPOS.LT.0  .OR. IBPOS.GE.NPCOL      ) THEN
            INFO = 17
          END IF
          IF( INFO.NE.0 ) GO TO 10
*
          IF( MYCOL.EQ.IBPOS )
     $      CALL PBDMATADD( ICONTXT, 'V', MP, NS, ZERO, DUMMY, 1, ALPHA,
     $                      B, LDB )
*
        ELSE
          IF( LDB.LT.MAX(1,M) .AND. ( BSPACE .OR.
     $          IBPOS.EQ.MYROW .OR. IBPOS.EQ.-1 )      ) THEN
            INFO = 14
          ELSE IF( IBPOS.LT.0  .OR. IBPOS.GE.NPROW     ) THEN
            INFO = 17
          END IF
          IF( INFO.NE.0 ) GO TO 10
*
          IF( MYROW.EQ.IBPOS )
     $      CALL PBDMATADD( ICONTXT, 'V', NS, MQ, ZERO, DUMMY, 1, ALPHA,
     $                      B, LDB )
        END IF
*
*       When alpha = zero, quick return.
*
        IF( ALPHA.EQ.ZERO ) RETURN
*
*                                    MTPROW (Top)
*                                        |
*       MLFCOL <- MYCOL -> MRTCOL,     MYROW
*       (Left)             (Right)       |
*                                   MBTROW (Bottom)
*
        MLFCOL = MOD( NPCOL+MYCOL-1, NPCOL )
        MRTCOL = MOD( MYCOL+1, NPCOL )
        MTPROW = MOD( NPROW+MYROW-1, NPROW )
        MBTROW = MOD( MYROW+1, NPROW )
        LCM    = ILCM( NPROW, NPCOL )
        BDATA  = .FALSE.
*
*       Start the operations.
*
        IF( UPPER ) THEN
          IF( LSIDE.AND.NOTRAN .OR. .NOT.(LSIDE.OR.NOTRAN) ) THEN
*
*           Form  B := Up( A ) \ alpha * B.
*            _       __________              _
*           | |      \_        |            | |
*           | |        \_      |            | |
*           |B|  :=      \_ A  |  \ alpha * |B|
*           | |            \_  |            | |
*           |_|              \_|            |_|
*
            IPT     = MP * NS + 1
            ICURROW = MOD( ICEIL(NDIM,NB)+IAROW-1, NPROW )
            ICURCOL = MOD( ICEIL(NDIM,NB)+IACOL-1, NPCOL )
*
            IF( LSIDE ) THEN
              IF( BSPACE ) THEN
                IF( MYCOL.EQ.IBPOS ) THEN
                  IF( MYCOL.NE.ICURCOL ) THEN
                    CALL DGESD2D( ICONTXT, MP, NS, B, LDB,
     $                            MYROW, ICURCOL )
                    CALL PBDMATADD( ICONTXT, 'G', MP, NS, ZERO, DUMMY,
     $                              1, ZERO, B, LDB )
                  END IF
                ELSE IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGERV2D( ICONTXT, MP, NS, B, LDB, MYROW, IBPOS )
                ELSE
                  CALL PBDMATADD( ICONTXT, 'G', MP, NS, ZERO, DUMMY, 1,
     $                            ZERO, B, LDB )
                END IF
                BDATA = .TRUE.
                IPT = 1
*
              ELSE
                IF( MYCOL.EQ.IBPOS ) THEN
                  IF( MYCOL.EQ.ICURCOL ) THEN
                     CALL PBDMATADD( ICONTXT, 'V', MP, NS, ONE, B, LDB,
     $                               ZERO, WORK, MP )
                  ELSE
                     CALL DGESD2D( ICONTXT, MP, NS, B, LDB,
     $                             MYROW, ICURCOL )
                  END IF
                ELSE IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGERV2D( ICONTXT, MP, NS, WORK, MP,
     $                          MYROW, IBPOS )
                END IF
              END IF
              IDEST = IBPOS
*
            ELSE
              CALL PBDTRAN( ICONTXT, 'Row', TRANSA, NS, NDIM, NB, B,
     $                      LDB, ZERO, WORK, MP, IBPOS, IACOL, IAROW,
     $                      ICURCOL, WORK(IPT) )
              IDEST = ICURCOL
            END IF
            LDW = MAX( 1, MP )
*
            IRPB  = MOD( NPCOL+IDEST-MYCOL-1, NPCOL )
            IRDB  = NB * MOD( IRPB+1, NPCOL )
            IRPB  = NB * IRPB
            NCOMM = NB * (NPCOL-1)
            KB    = MOD( NDIM, NB )
            IF( KB.EQ.0 ) KB = NB
*
            II = MP - NB + 1
            IF( MYROW.EQ.ICURROW ) II = MP - KB + 1
            IN = II
            JJ = MQ - NB + 1
            IF( MYCOL.EQ.ICURCOL ) JJ = MQ - KB + 1
            JB = KB
*
*           If B can be used as a working space,
*
            IF( BDATA ) THEN
              DO 20 J = 1, NDIM, NB
                NLENG  = NDIM - J - KB + 1
                NXTROW = MOD( NPROW+ICURROW-1, NPROW )
                NXTCOL = MOD( NPCOL+ICURCOL-1, NPCOL )
                IF( MYROW.EQ.ICURROW ) IN = II - NB
*
                IF( MYCOL.EQ.ICURCOL ) THEN
*
*                 Receive updated blocks from previous column of
*                 processes
*
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG+JB, NCOMM), NB, ICURROW,
     $                              MYROW, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, WORK(IPT), NPART,
     $                            MYROW, MRTCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NPART, NS, ONE,
     $                              WORK(IPT), NPART, ONE,
     $                              B(II+NB-NPART,1), LDB )
                  END IF
*
*                 B(II,1) <== A(II,JJ) \ B(II,1), ( B(II,1) = WORK(II) )
*                 where A(II,JJ) is a upper triangular matrix
*
                  IF( MYROW.EQ.ICURROW ) THEN
                    CALL DTRSM( 'Left', 'Upper', 'No', DIAG, JB, NS,
     $                          ONE, A(II,JJ), LDA, B(II,1), LDB )
                    CALL PBDMATADD( ICONTXT, 'G', JB, NS, ONE, B(II,1),
     $                              LDB, ZERO, WORK(IPT), JB )
                    CALL DGEBS2D( ICONTXT, 'Col', 'D-ring', JB, NS,
     $                            WORK(IPT), JB )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Col', 'D-ring', JB, NS,
     $                            WORK(IPT), JB, ICURROW, MYCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.0 ) THEN
*
*                   Update the (NPCOL-1) blocks first
*
                    NREST = MIN( NLENG, NCOMM )
                    NPART = NUMROC( NREST, NB, ICURROW, MYROW+1, NPROW )
                    IIN   = IN + NB - NPART
*
                    CALL DGEMM( 'No', 'No', NPART, NS, JB, -ONE,
     $                          A(IIN,JJ), LDA, WORK(IPT), JB, ONE,
     $                          B(IIN,1), LDB )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPCOL.GT.1 )
     $                CALL DGESD2D( ICONTXT, NPART, NS, B(IIN,1), LDB,
     $                              MYROW, MLFCOL)
*
*                   Update the rest of the matrix
*
                    CALL DGEMM( 'No', 'No', IIN-1, NS, JB, -ONE,
     $                          A(1,JJ), LDA, WORK(IPT), JB, ONE, B,
     $                          LDB )
                  END IF
*
*                 Send the solution blocks to destination (IDEST column)
*
                  JJN = J + KB - 1
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, JJN-JB), NB, MYROW,
     $                              ICURROW+1, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, B(II+NB,1), LDB,
     $                            MYROW, MRTCOL)
                  END IF
*
                  IF( NLENG.GT.0 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB, JJN), NB, MYROW, ICURROW,
     $                              NPROW )
                    CALL DGESD2D( ICONTXT, NPART, NS, B(IN+NB,1), LDB,
     $                            MYROW, MLFCOL )
                  END IF
*
                  JJ = JJ - NB
                END IF
*
                II = IN
                JB = NB
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   20         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURCOL = MOD( ICURCOL+1, NPCOL )
              IF( ICURCOL.NE.IDEST ) THEN
                KDIST = MOD( NPCOL+IDEST-ICURCOL, NPCOL )
                NPART = NUMROC( MIN(NDIM, KDIST*NB), NB, MYROW, IAROW,
     $                          NPROW )
*
                IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGESD2D( ICONTXT, NPART, NS, B, LDB,
     $                          MYROW, IDEST )
                ELSE IF( MYCOL.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NPART, NS, B, LDB,
     $                          MYROW, ICURCOL )
                END IF
              END IF
*
*           If B can't be used as a working space,
*
            ELSE
              IF( MYCOL.NE.ICURCOL )
     $          CALL PBDMATADD( ICONTXT, 'G', MP, NS, ZERO, DUMMY, 1,
     $                          ZERO, WORK, MP )
*
              DO 30 J = 1, NDIM, NB
                NLENG  = NDIM - J - KB + 1
                NXTROW = MOD( NPROW+ICURROW-1, NPROW )
                NXTCOL = MOD( NPCOL+ICURCOL-1, NPCOL )
                IF( MYROW.EQ.ICURROW ) IN = II - NB
*
                IF( MYCOL.EQ.ICURCOL ) THEN
*
*                 Receive updated blocks from previous column of
*                 processes
*
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG+JB, NCOMM), NB, ICURROW,
     $                              MYROW, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, WORK(IPT), NPART,
     $                            MYROW, MRTCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NPART, NS, ONE,
     $                              WORK(IPT), NPART, ONE,
     $                              WORK(II+NB-NPART), MP )
                  END IF
*
*                 B(II,1) <== A(II,JJ) \ B(II,1), ( B(II,1) = WORK(II) )
*                 where A(II,JJ) is a upper triangular matrix
*
                  IF( MYROW.EQ.ICURROW ) THEN
                    CALL DTRSM( 'Left', 'Upper', 'No', DIAG, JB, NS,
     $                          ONE, A(II,JJ), LDA, WORK(II), LDW )
                    CALL PBDMATADD( ICONTXT, 'G', JB, NS, ONE, WORK(II),
     $                              MP, ZERO, WORK(IPT), JB )
                    CALL DGEBS2D( ICONTXT, 'Col', 'D-ring', JB, NS,
     $                            WORK(IPT), JB )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Col', 'D-ring', JB, NS,
     $                            WORK(IPT), JB, ICURROW, MYCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.0 ) THEN
*
*                   Update the (NPCOL-1) blocks first
*
                    NREST = MIN( NLENG, NCOMM )
                    NPART = NUMROC( NREST, NB, ICURROW, MYROW+1, NPROW )
                    IIN   = IN + NB - NPART
*
                    CALL DGEMM( 'No', 'No', NPART, NS, JB, -ONE,
     $                          A(IIN,JJ), LDA, WORK(IPT), JB, ONE,
     $                          WORK(IIN), LDW )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPCOL.GT.1 )
     $                CALL DGESD2D( ICONTXT, NPART, NS, WORK(IIN), MP,
     $                              MYROW, MLFCOL)
*
*                   Update the rest of the matrix
*
                    CALL DGEMM( 'No','No', IIN-1, NS, JB, -ONE, A(1,JJ),
     $                          LDA, WORK(IPT),JB, ONE, WORK, LDW )
                  END IF
*
*                 Send the solution blocks to destination (IDEST column)
*
                  JJN = J + KB - 1
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, JJN-JB), NB, MYROW,
     $                              ICURROW+1, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, WORK(II+NB), MP,
     $                            MYROW, MRTCOL )
                  END IF
*
                  IF( NLENG.GT.0 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB, JJN), NB, MYROW, ICURROW,
     $                              NPROW )
                    CALL DGESD2D( ICONTXT, NPART, NS, WORK(IN+NB), MP,
     $                            MYROW, MLFCOL )
                  END IF
*
                  JJ = JJ - NB
                END IF
*
                II = IN
                JB = NB
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   30         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURCOL = MOD( ICURCOL+1, NPCOL )
              IF( ICURCOL.NE.IDEST ) THEN
                KDIST = MOD( NPCOL+IDEST-ICURCOL, NPCOL )
                NPART = NUMROC( MIN(NDIM, KDIST*NB), NB, MYROW, IAROW,
     $                          NPROW )
                IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGESD2D( ICONTXT, NPART, NS, WORK, MP,
     $                          MYROW, IDEST )
                ELSE IF( MYCOL.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NPART, NS, WORK, MP,
     $                          MYROW, ICURCOL )
                END IF
              END IF
            END IF
*
            IF( LSIDE ) THEN
              IF( .NOT.BDATA .AND. MYCOL.EQ.IDEST )
     $          CALL PBDMATADD( ICONTXT, 'V', MP, NS, ONE, WORK, MP,
     $                          ZERO, B, LDB )
            ELSE
              CALL PBDTRAN( ICONTXT, 'Col', TRANSA, NDIM, NS, NB, WORK,
     $                      MP, ZERO, B, LDB, IAROW, IDEST, IBPOS,
     $                      IACOL, WORK(IPT) )
            END IF
*
          ELSE IF( ( LSIDE .AND. .NOT.NOTRAN ) .OR.
     $             ( .NOT.LSIDE .AND. NOTRAN ) ) THEN
*
*
*           Form  B := alpha * B / Up( A ).
*                                                  __________
*                                                  \_        |
*            __________              __________      \_      |
*           |_____B____| := alpha * |_____B____| /     \_ A  |
*                                                        \_  |
*                                                          \_|
*
            IPT     = NS * MQ + 1
            ICURROW = IAROW
            ICURCOL = IACOL
*
            IF( LSIDE ) THEN
              CALL PBDTRAN( ICONTXT, 'Col', TRANSA, NDIM, NS, NB, B,
     $                      LDB, ZERO, WORK, NS, IAROW, IBPOS, ICURROW,
     $                      IACOL, WORK(IPT) )
              IDEST = ICURROW
*
            ELSE
              IF( BSPACE ) THEN
                IF( MYROW.EQ.IBPOS ) THEN
                  IF( MYROW.NE.ICURROW ) THEN
                    CALL DGESD2D( ICONTXT, NS, MQ, B, LDB,
     $                            ICURROW, MYCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ZERO, DUMMY,
     $                              1, ZERO, B, LDB )
                  END IF
                ELSE IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGERV2D( ICONTXT, NS, MQ, B, LDB, IBPOS, MYCOL )
                ELSE
                  CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ZERO, DUMMY, 1,
     $                            ZERO, B, LDB )
                END IF
                BDATA = .TRUE.
                IPT = 1
*
              ELSE
                IF( MYROW.EQ.IBPOS ) THEN
                  IF( MYROW.EQ.ICURROW ) THEN
                    CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ONE, B, LDB,
     $                              ZERO, WORK, NS )
                  ELSE
                    CALL DGESD2D( ICONTXT, NS, MQ, B, LDB,
     $                            ICURROW, MYCOL )
                  END IF
                ELSE IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGERV2D( ICONTXT, NS, MQ, WORK, NS,
     $                          IBPOS, MYCOL )
                END IF
              END IF
              IDEST = IBPOS
            END IF
*
            IRPB  = MOD( NPROW+MYROW-IDEST-1, NPROW )
            IRDB  = NB * MOD( IRPB+1, NPROW )
            IRPB  = NB * IRPB
            NCOMM = NB * (NPROW-1)
*
            II = 1
            JJ = 0
            JN = 0
*
*           If B can be used as a working space,
*
            IF( BDATA ) THEN
              DO 40 J = 1, NDIM, NB
                NLENG  = NDIM - J + 1
                JB     = MIN( NLENG, NB )
                NXTROW = MOD( ICURROW+1, NPROW )
                NXTCOL = MOD( ICURCOL+1, NPCOL )
                IF( MYCOL.EQ.ICURCOL ) JN = JJ + JB
*
                IF( MYROW.EQ.ICURROW ) THEN
*
*                 Receive updated blocks from previous row of processes
*
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG, NCOMM), NB, MYCOL,
     $                              ICURCOL, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART, WORK(IPT), NS,
     $                            MTPROW, MYCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NS, NPART, ONE,
     $                              WORK(IPT), NS, ONE, B(1,JJ+1), LDB )
                  END IF
*
*                 B(1,JJ+1) <== A(II,JJ+1) \ B(1,JJ+1),
*                                          ( B(1,JJ+1) = WORK(JJN) )
*                 where A(II,JJ+1) is a upper triangular matrix
*
                  IF( MYCOL.EQ.ICURCOL ) THEN
                    CALL DTRSM( 'Right', 'Upper', 'No', DIAG, NS, JB,
     $                          ONE, A(II,JJ+1), LDA, B(1,JJ+1), LDB )
                    CALL PBDMATADD( ICONTXT, 'G', NS, JB, ONE,
     $                              B(1,JJ+1), LDB, ZERO, WORK(IPT),
     $                              NS )
                    CALL DGEBS2D( ICONTXT, 'Row', 'I-ring', NS, JB,
     $                            WORK(IPT), NS )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Row', 'I-ring', NS, JB,
     $                            WORK(IPT), NS, MYROW, ICURCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.JB ) THEN
*
*                   Update the (NPROW-1) blocks first
*
                    NREST = MIN( NLENG-JB, NCOMM )
                    NPART = NUMROC( NREST, NB, MYCOL, ICURCOL+1, NPCOL )
*
                    CALL DGEMM( 'No', 'No', NS, NPART, JB, -ONE,
     $                          WORK(IPT), NS, A(II,JN+1), LDA, ONE,
     $                          B(1,JN+1), LDB )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPROW.GT.1 )
     $                CALL DGESD2D( ICONTXT, NS, NPART, B(1,JN+1), LDB,
     $                              MBTROW, MYCOL )
*
*                   Update the rest of the matrix
*
                    IPART = NUMROC( NLENG-JB-NREST, NB, MYCOL+LCM,
     $                              ICURCOL+NPROW, NPCOL )
                    CALL DGEMM( 'No', 'No', NS,IPART, JB, -ONE,
     $                          WORK(IPT), NS, A(II,JN+NPART+1), LDA,
     $                          ONE, B(1,JN+NPART+1), LDB )
                  END IF
*
*                 Send the solution blocks to destination (IDEST row)
*
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, J-1), NB, ICURCOL,
     $                              MYCOL+1, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART, B(1,JJ-NPART+1),
     $                            LDB, MTPROW, MYCOL )
                  END IF
*
                  IF( NLENG.GT.JB .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB, J+JB-1), NB, ICURCOL,
     $                              MYCOL, NPCOL )
                    CALL DGESD2D( ICONTXT, NS, NPART, B(1,JN-NPART+1),
     $                            LDB, MBTROW, MYCOL )
                  END IF
*
                  II = II + JB
                END IF
*
                JJ = JN
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   40         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURROW = MOD( NPROW+ICURROW-1, NPROW )
              IF( ICURROW.NE.IDEST ) THEN
                KDIST = MOD( NPROW+ICURROW-IDEST-1, NPROW )
                IF( ICEIL(NDIM,NB).GT.MOD(IDEST-IAROW+NPROW,NPROW) )
     $              THEN
                  NPART = NUMROC( KDIST*NB+JB, NB, MYCOL+KDIST,
     $                            ICURCOL-1, NPCOL )
                ELSE
                  NPART = NUMROC( NDIM, NB, MYCOL, IACOL, NPCOL )
                END IF
*
                IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGESD2D( ICONTXT, NS, NPART, B(1,JJ-NPART+1),
     $                          LDB, IDEST, MYCOL )
                ELSE IF( MYROW.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NS, NPART, B(1,JJ-NPART+1),
     $                          LDB, ICURROW, MYCOL )
                END IF
              END IF
*
*           If B can't be used as a working space,
*
            ELSE
              IF( MYROW.NE.ICURROW )
     $          CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ZERO, DUMMY, 1,
     $                          ZERO, WORK, NS )
*
              DO 50 J = 1, NDIM, NB
                NLENG  = NDIM - J + 1
                JB     = MIN( NLENG, NB )
                NXTROW = MOD( ICURROW+1, NPROW )
                NXTCOL = MOD( ICURCOL+1, NPCOL )
                IF( MYCOL.EQ.ICURCOL ) JN = JJ + JB
*
                IF( MYROW.EQ.ICURROW ) THEN
*
*                 Receive updated blocks from previous row of processes
*
                  JJN = JJ * NS + 1
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG, NCOMM), NB, MYCOL,
     $                              ICURCOL, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART, WORK(IPT), NS,
     $                            MTPROW, MYCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NS, NPART, ONE,
     $                              WORK(IPT), NS, ONE, WORK(JJN), NS )
                  END IF
*
*                 B(1,JJ+1) <== A(II,JJ+1) \ B(1,JJ+1),
*                                          ( B(1,JJ+1) = WORK(JJN) )
*                 where A(II,JJ+1) is a upper triangular matrix
*
                  IF( MYCOL.EQ.ICURCOL ) THEN
                    CALL DTRSM( 'Right', 'Upper', 'No', DIAG, NS, JB,
     $                          ONE, A(II,JJ+1), LDA, WORK(JJN), NS )
                    CALL PBDMATADD( ICONTXT, 'G', NS, JB, ONE,
     $                              WORK(JJN), NS, ZERO, WORK(IPT), NS )
                    CALL DGEBS2D( ICONTXT, 'Row', 'I-ring', NS, JB,
     $                            WORK(IPT), NS )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Row', 'I-ring', NS, JB,
     $                            WORK(IPT), NS, MYROW, ICURCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.JB ) THEN
*
*                   Update the (NPROW-1) blocks first
*
                    NREST = MIN( NLENG-JB, NCOMM )
                    NPART = NUMROC( NREST, NB, MYCOL, ICURCOL+1, NPCOL )
*
                    CALL DGEMM( 'No', 'No', NS, NPART, JB, -ONE,
     $                          WORK(IPT), NS, A(II,JN+1), LDA, ONE,
     $                          WORK(JN*NS+1), NS )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPROW.GT.1 )
     $                CALL DGESD2D( ICONTXT, NS, NPART, WORK(JN*NS+1),
     $                              NS, MBTROW, MYCOL )
*
*                   Update the rest of the matrix
*
                    IPART = NUMROC( NLENG-JB-NREST, NB, MYCOL+LCM,
     $                              ICURCOL+NPROW, NPCOL )
                    CALL DGEMM( 'No', 'No', NS,IPART, JB, -ONE,
     $                          WORK(IPT), NS, A(II,JN+NPART+1), LDA,
     $                          ONE, WORK((JN+NPART)*NS+1), NS )
                  END IF
*
*                 Send the solution blocks to destination (IDEST row)
*
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, J-1), NB, ICURCOL,
     $                              MYCOL+1, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART,
     $                            WORK((JJ-NPART)*NS+1), NS,
     $                            MTPROW, MYCOL )
                  END IF
*
                  IF( NLENG.GT.JB .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB, J+JB-1), NB, ICURCOL,
     $                              MYCOL, NPCOL )
                    CALL DGESD2D( ICONTXT, NS, NPART,
     $                            WORK((JN-NPART)*NS+1), NS,
     $                            MBTROW, MYCOL )
                  END IF
*
                  II = II + JB
                END IF
*
                JJ = JN
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   50         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURROW = MOD( NPROW+ICURROW-1, NPROW )
              IF( ICURROW.NE.IDEST ) THEN
                KDIST = MOD( NPROW+ICURROW-IDEST-1, NPROW )
                IF( ICEIL(NDIM,NB).GT.MOD(IDEST-IAROW+NPROW,NPROW) )
     $              THEN
                  NPART = NUMROC( KDIST*NB+JB, NB, MYCOL+KDIST,
     $                            ICURCOL-1, NPCOL )
                ELSE
                  NPART = NUMROC( NDIM, NB, MYCOL, IACOL, NPCOL )
                END IF
*
                IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGESD2D( ICONTXT, NS, NPART,
     $                           WORK((JJ-NPART)*NS+1), NS,
     $                           IDEST, MYCOL )
                ELSE IF( MYROW.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NS, NPART,
     $                          WORK((JJ-NPART)*NS+1), NS,
     $                          ICURROW, MYCOL )
                END IF
              END IF
            END IF
*
            IF( LSIDE ) THEN
              CALL PBDTRAN( ICONTXT, 'Row', TRANSA, NS, NDIM, NB, WORK,
     $                      NS, ZERO, B, LDB, IDEST, IACOL, IAROW,
     $                      IBPOS, WORK(IPT) )
            ELSE
              IF( .NOT.BDATA .AND. MYROW.EQ.IDEST )
     $          CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ONE, WORK, NS,
     $                          ZERO, B, LDB )
            END IF
          END IF
*
*       if ( LSAME( UPLO, 'L' ) then
*
        ELSE
*
          IF( LSIDE.AND.NOTRAN .OR. .NOT.(LSIDE.OR.NOTRAN) ) THEN
*
*           Form  B := Lo( A ) \ alpha * B.
*            _        _                      _
*           | |      | \_                   | |
*           | |      |   \_                 | |
*           |B|  :=  |  A  \_     \ alpha * |B|
*           | |      |       \_             | |
*           |_|      |_________|            |_|
*
            IPT     = MP * NS + 1
            ICURROW = IAROW
            ICURCOL = IACOL
*
            IF( LSIDE ) THEN
              IF( BSPACE ) THEN
                IF( MYCOL.EQ.IBPOS ) THEN
                  IF( MYCOL.NE.ICURCOL ) THEN
                    CALL DGESD2D( ICONTXT, MP, NS, B, LDB,
     $                            MYROW, ICURCOL )
                    CALL PBDMATADD( ICONTXT, 'G', MP, NS, ZERO, DUMMY,
     $                              1, ZERO, B, LDB )
                  END IF
                ELSE IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGERV2D( ICONTXT, MP, NS, B, LDB, MYROW, IBPOS )
                ELSE
                  CALL PBDMATADD( ICONTXT, 'G', MP, NS, ZERO, DUMMY, 1,
     $                            ZERO, B, LDB )
                END IF
                BDATA = .TRUE.
                IPT = 1
*
              ELSE
                IF( MYCOL.EQ.IBPOS ) THEN
                  IF( MYCOL.EQ.ICURCOL ) THEN
                    CALL PBDMATADD( ICONTXT, 'V', MP, NS, ONE, B, LDB,
     $                              ZERO, WORK, MP )
                  ELSE
                    CALL DGESD2D( ICONTXT, MP, NS, B, LDB,
     $                            MYROW, ICURCOL )
                  END IF
                ELSE IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGERV2D( ICONTXT, MP, NS, WORK, MP,
     $                          MYROW, IBPOS )
                END IF
              END IF
              IDEST = IBPOS
*
            ELSE
              CALL PBDTRAN( ICONTXT, 'Row', TRANSA, NS, NDIM, NB, B,
     $                      LDB, ZERO, WORK, MP, IBPOS, IACOL, IAROW,
     $                      ICURCOL, WORK(IPT) )
              IDEST = IACOL
            END IF
*
            LDW = MAX( 1, MP )
            IRPB  = MOD( NPCOL+MYCOL-IDEST-1, NPCOL )
            IRDB  = NB * MOD( IRPB+1, NPCOL )
            IRPB  = NB * IRPB
            NCOMM = NB * (NPCOL-1)
*
            II = 1
            IN = 1
            JJ = 1
*
*           If B can be used as a working space,
*
            IF( BDATA ) THEN
              DO 60 J = 1, NDIM, NB
                NLENG  = NDIM - J + 1
                JB     = MIN( NLENG, NB )
                NXTROW = MOD( ICURROW+1, NPROW )
                NXTCOL = MOD( ICURCOL+1, NPCOL )
                IF( MYROW.EQ.ICURROW ) IN = II + JB
*
                IF( MYCOL.EQ.ICURCOL ) THEN
*
*                 Receive updated blocks from previous column of
*                 processes
*
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG, NCOMM), NB, MYROW,
     $                              ICURROW, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, WORK(IPT), NPART,
     $                            MYROW, MLFCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NPART, NS, ONE,
     $                              WORK(IPT), NPART, ONE, B(II,1),
     $                              LDB )
                  END IF
*
*                 B(II,1) <== A(II,JJ) \ B(II,1), ( B(II,1) = WORK(II) )
*                 where A(II,JJ) is a lower triangular matrix
*
                  IF( MYROW.EQ.ICURROW ) THEN
                    CALL DTRSM( 'Left', 'Lower', 'No', DIAG, JB, NS,
     $                          ONE, A(II,JJ), LDA, B(II,1), LDB )
                    CALL PBDMATADD( ICONTXT, 'G', JB, NS, ONE, B(II,1),
     $                              LDB, ZERO, WORK(IPT), JB )
                    CALL DGEBS2D( ICONTXT, 'Col', 'I-ring', JB, NS,
     $                            WORK(IPT), JB )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Col', 'I-ring', JB, NS,
     $                            WORK(IPT), JB, ICURROW, MYCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.JB ) THEN
*
*                   Update the (NPCOL-1) blocks first
*
                    NREST = MIN( NLENG-JB, NCOMM )
                    NPART = NUMROC( NREST, NB, MYROW, ICURROW+1, NPROW )
*
                    CALL DGEMM( 'No', 'No', NPART, NS, JB, -ONE,
     $                          A(IN,JJ), LDA, WORK(IPT), JB, ONE,
     $                          B(IN,1), LDB )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPCOL.GT.1 )
     $                CALL DGESD2D( ICONTXT, NPART, NS, B(IN,1), LDB,
     $                              MYROW, MRTCOL )
*
*                   Update the rest of the matrix
*
                    IPART = NUMROC( NLENG-JB-NREST, NB, MYROW+LCM,
     $                              ICURROW+NPCOL, NPROW )
                    CALL DGEMM( 'No', 'No', IPART, NS, JB, -ONE,
     $                          A(IN+NPART,JJ), LDA, WORK(IPT), JB,
     $                          ONE, B(IN+NPART,1), LDB )
                  END IF
*
*                 Send the solution blocks to destination (IDEST column)
*
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, J-1), NB, ICURROW,
     $                              MYROW+1, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, B(II-NPART,1),
     $                            LDB, MYROW, MLFCOL )
                  END IF
*
                  IF( NLENG.GT.JB .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB,J+JB-1), NB, ICURROW,
     $                              MYROW, NPROW )
                    CALL DGESD2D( ICONTXT, NPART, NS, B(IN-NPART,1),
     $                            LDB, MYROW, MRTCOL )
                  END IF
*
                  JJ = JJ + JB
                END IF
*
                II = IN
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   60         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURCOL = MOD( NPCOL+ICURCOL-1, NPCOL )
              IF( ICURCOL.NE.IDEST ) THEN
                KDIST = MOD( NPCOL+ICURCOL-IDEST-1, NPCOL )
                IF( ICEIL(NDIM,NB).GT.MOD(IDEST-IACOL+NPCOL,NPCOL) )
     $              THEN
                  NPART = NUMROC( KDIST*NB+JB, NB, MYROW+KDIST,
     $                            ICURROW-1, NPROW )
                ELSE
                  NPART = NUMROC( NDIM, NB, MYROW, IAROW, NPROW )
                END IF
*
                IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGESD2D( ICONTXT, NPART, NS, B(II-NPART,1), LDB,
     $                          MYROW, IDEST )
                ELSE IF( MYCOL.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NPART, NS, B(II-NPART,1), LDB,
     $                          MYROW, ICURCOL )
                END IF
              END IF
*
*           If B can't be used as a working space,
*
            ELSE
              IF( MYCOL.NE.ICURCOL )
     $          CALL PBDMATADD( ICONTXT, 'G', MP, NS, ZERO, DUMMY, 1,
     $                          ZERO, WORK, MP )
*
              DO 70 J = 1, NDIM, NB
                NLENG  = NDIM - J + 1
                JB     = MIN( NLENG, NB )
                NXTROW = MOD( ICURROW+1, NPROW )
                NXTCOL = MOD( ICURCOL+1, NPCOL )
                IF( MYROW.EQ.ICURROW ) IN = II + JB
*
                IF( MYCOL.EQ.ICURCOL ) THEN
*
*                 Receive updated blocks from previous column of
*                 processes
*
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG, NCOMM), NB, MYROW,
     $                              ICURROW, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, WORK(IPT), NPART,
     $                            MYROW, MLFCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NPART, NS, ONE,
     $                              WORK(IPT), NPART, ONE, WORK(II),
     $                              MP )
                  END IF
*
*                 B(II,1) <== A(II,JJ) \ B(II,1), ( B(II,1) = WORK(II) )
*                 where A(II,JJ) is a lower triangular matrix
*
                  IF( MYROW.EQ.ICURROW ) THEN
                    CALL DTRSM( 'Left', 'Lower', 'No', DIAG, JB, NS,
     $                          ONE, A(II,JJ), LDA, WORK(II), LDW )
                    CALL PBDMATADD( ICONTXT, 'G', JB, NS, ONE, WORK(II),
     $                              MP, ZERO, WORK(IPT), JB )
                    CALL DGEBS2D( ICONTXT, 'Col', 'I-ring', JB, NS,
     $                            WORK(IPT), JB )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Col', 'I-ring', JB, NS,
     $                            WORK(IPT), JB, ICURROW, MYCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.JB ) THEN
*
*                   Update the (NPCOL-1) blocks first
*
                    NREST = MIN( NLENG-JB, NCOMM )
                    NPART = NUMROC( NREST, NB, MYROW, ICURROW+1, NPROW )
*
                    CALL DGEMM( 'No', 'No', NPART, NS, JB, -ONE,
     $                          A(IN,JJ), LDA, WORK(IPT), JB, ONE,
     $                          WORK(IN), LDW )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPCOL.GT.1 )
     $                CALL DGESD2D( ICONTXT, NPART, NS, WORK(IN), MP,
     $                              MYROW, MRTCOL )
*
*                   Update the rest of the matrix
*
                    IPART = NUMROC( NLENG-JB-NREST, NB, MYROW+LCM,
     $                              ICURROW+NPCOL, NPROW )
                    CALL DGEMM( 'No', 'No', IPART, NS, JB, -ONE,
     $                          A(IN+NPART,JJ), LDA, WORK(IPT), JB,
     $                          ONE, WORK(IN+NPART), LDW )
                  END IF
*
*                 Send the solution blocks to destination (IDEST column)
*
                  IF( J.GT.1 .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, J-1), NB, ICURROW,
     $                              MYROW+1, NPROW )
                    CALL DGERV2D( ICONTXT, NPART, NS, WORK(II-NPART),
     $                            MP, MYROW, MLFCOL )
                  END IF
*
                  IF( NLENG.GT.JB .AND. NPCOL.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB,J+JB-1), NB, ICURROW,
     $                              MYROW, NPROW )
                    CALL DGESD2D( ICONTXT, NPART, NS, WORK(IN-NPART),
     $                            MP, MYROW, MRTCOL )
                  END IF
*
                  JJ = JJ + JB
                END IF
*
                II = IN
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   70         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURCOL = MOD( NPCOL+ICURCOL-1, NPCOL )
              IF( ICURCOL.NE.IDEST ) THEN
                KDIST = MOD( NPCOL+ICURCOL-IDEST-1, NPCOL )
                IF( ICEIL(NDIM,NB).GT.MOD(IDEST-IACOL+NPCOL,NPCOL) )
     $              THEN
                  NPART = NUMROC( KDIST*NB+JB, NB, MYROW+KDIST,
     $                            ICURROW-1, NPROW )
                ELSE
                  NPART = NUMROC( NDIM, NB, MYROW, IAROW, NPROW )
                END IF
*
                IF( MYCOL.EQ.ICURCOL ) THEN
                  CALL DGESD2D( ICONTXT, NPART, NS, WORK(II-NPART), MP,
     $                          MYROW, IDEST )
                ELSE IF( MYCOL.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NPART, NS, WORK(II-NPART), MP,
     $                          MYROW, ICURCOL )
                END IF
              END IF
            END IF
*
            IF( LSIDE ) THEN
              IF( .NOT.BDATA .AND. MYCOL.EQ.IDEST )
     $          CALL PBDMATADD( ICONTXT, 'V', MP, NS, ONE, WORK, MP,
     $                          ZERO, B, LDB )
            ELSE
              CALL PBDTRAN( ICONTXT, 'Col', TRANSA, NDIM, NS, NB, WORK,
     $                      MP, ZERO, B, LDB, IAROW, IDEST, IBPOS,
     $                      IACOL, WORK(IPT) )
            END IF
*
          ELSE IF( ( LSIDE .AND. .NOT.NOTRAN ) .OR.
     $             ( .NOT.LSIDE .AND. NOTRAN ) ) THEN
*
*           Form  B := alpha * B / Lo( A ).
*                                                  _
*                                                 | \_
*           __________              __________    |   \_
*          |_____B____| := alpha * |_____B____| / |  A  \_
*                                                 |       \_
*                                                 |_________|
*
            IPT     = NS * MQ + 1
            ICURROW = MOD( ICEIL(NDIM,NB)+IAROW-1, NPROW )
            ICURCOL = MOD( ICEIL(NDIM,NB)+IACOL-1, NPCOL )
*
            IF( LSIDE ) THEN
              CALL PBDTRAN( ICONTXT, 'Col', TRANSA, NDIM, NS, NB, B,
     $                      LDB, ZERO, WORK, NS, IAROW, IBPOS, ICURROW,
     $                      IACOL, WORK(IPT) )
              IDEST = ICURROW
*
            ELSE
              IF( BSPACE ) THEN
                IF( MYROW.EQ.IBPOS ) THEN
                  IF( MYROW.NE.ICURROW ) THEN
                    CALL DGESD2D( ICONTXT, NS, MQ, B, LDB,
     $                            ICURROW, MYCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ZERO, DUMMY,
     $                              1, ZERO, B, LDB )
                  END IF
                ELSE IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGERV2D( ICONTXT, NS, MQ, B, LDB, IBPOS, MYCOL )
                ELSE
                  CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ZERO, DUMMY, 1,
     $                            ZERO, B, LDB )
                END IF
                BDATA = .TRUE.
                IPT = 1
*
              ELSE
                IF( MYROW.EQ.IBPOS ) THEN
                  IF( MYROW.EQ.ICURROW ) THEN
                    CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ONE, B, LDB,
     $                              ZERO, WORK, NS )
                  ELSE
                    CALL DGESD2D( ICONTXT, NS, MQ, B, LDB,
     $                            ICURROW, MYCOL )
                  END IF
                ELSE IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGERV2D( ICONTXT, NS, MQ, WORK, NS,
     $                          IBPOS, MYCOL )
                END IF
              END IF
              IDEST = IBPOS
            END IF
*
            IRPB  = MOD( NPROW+IDEST-MYROW-1, NPROW )
            IRDB  = NB * MOD( IRPB+1, NPROW )
            IRPB  = NB * IRPB
            NCOMM = NB * (NPROW-1)
            KB    = MOD( NDIM, NB )
            IF( KB.EQ.0 ) KB = NB
*
            II = MP - NB + 1
            IF( MYROW.EQ.ICURROW ) II = MP - KB + 1
            JJ = MQ - NB
            IF( MYCOL.EQ.ICURCOL ) JJ = MQ - KB
            JN = JJ
            JB = KB
*
*           If B can be used as a working space,
*
            IF( BDATA ) THEN
              DO 80 J = 1, NDIM, NB
                NLENG  = NDIM - J - KB + 1
                NXTROW = MOD( NPROW+ICURROW-1, NPROW )
                NXTCOL = MOD( NPCOL+ICURCOL-1, NPCOL )
                IF( MYCOL.EQ.ICURCOL ) JN = JJ - NB
*
                IF( MYROW.EQ.ICURROW ) THEN
*
*                 Receive updated blocks from previous row of processes
*
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG+JB, NCOMM), NB, ICURCOL,
     $                              MYCOL, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART, WORK(IPT), NS,
     $                            MBTROW, MYCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NS, NPART, ONE,
     $                              WORK(IPT), NS, ONE,
     $                              B(1,JJ+NB-NPART+1), LDB )
                  END IF
*
*                 B(1,JJ+1) <== A(II,JJ+1) / B(1,JJ+1)
*                                          ( B(1,JJ+1) = WORK(JJ*NS+1) )
*                 where A(II,JJ+1) is a lower triangular matrix
*
                  IF( MYCOL.EQ.ICURCOL ) THEN
                    CALL DTRSM( 'Right', 'Lower', 'No', DIAG, NS, JB,
     $                          ONE, A(II,JJ+1), LDA, B(1,JJ+1), LDB )
                    CALL PBDMATADD( ICONTXT, 'G', NS, JB, ONE,
     $                             B(1,JJ+1), LDB, ZERO, WORK(IPT), NS )
                    CALL DGEBS2D( ICONTXT, 'Row', 'D-ring', NS, JB,
     $                            WORK(IPT), NS )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Row', 'D-ring', NS, JB,
     $                            WORK(IPT), NS, MYROW, ICURCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.0 ) THEN
*
*                   Update the (NPROW-1) blocks first
*
                    NREST = MIN( NLENG, NCOMM )
                    NPART = NUMROC( NREST, NB, ICURCOL, MYCOL+1, NPCOL )
                    JJN   = JN + NB - NPART
                    CALL DGEMM( 'No', 'No', NS, NPART, JB, -ONE,
     $                          WORK(IPT), NS, A(II,JJN+1), LDA, ONE,
     $                          B(1,JJN+1), LDB )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPROW.GT.1 )
     $                CALL DGESD2D( ICONTXT, NS, NPART, B(1,JJN+1), LDB,
     $                              MTPROW, MYCOL )
*
*                   Update the rest of the matrix
*
                    CALL DGEMM( 'No','No', NS, JJN, JB, -ONE, WORK(IPT),
     $                          NS, A(II,1), LDA, ONE, B, LDB )
                  END IF
*
*                 Send the solution blocks to destination (IDEST row)
*
                  JJN = J + KB - 1
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, JJN-JB), NB, MYCOL,
     $                              ICURCOL+1, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART, B(1,JJ+NB+1), LDB,
     $                            MBTROW, MYCOL )
                  END IF
*
                  IF( NLENG.GT.0 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB, JJN), NB, MYCOL, ICURCOL,
     $                              NPCOL )
                    CALL DGESD2D( ICONTXT, NS, NPART, B(1,JN+NB+1), LDB,
     $                            MTPROW, MYCOL )
                  END IF
*
                  II = II - NB
                END IF
*
                JJ = JN
                JB = NB
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   80         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURROW = MOD( ICURROW+1, NPROW )
              IF( ICURROW.NE.IDEST ) THEN
                KDIST = MOD( NPROW+IDEST-ICURROW, NPROW )
                NPART = NUMROC( MIN(NDIM, KDIST*NB), NB, MYCOL, IACOL,
     $                          NPCOL )
                IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGESD2D( ICONTXT, NS, NPART, B, LDB,
     $                          IDEST, MYCOL )
                ELSE IF( MYROW.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NS, NPART, B, LDB,
     $                          ICURROW, MYCOL )
                END IF
              END IF
*
*           If B can't be used as a working space,
*
            ELSE
              IF( MYROW.NE.ICURROW )
     $          CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ZERO, DUMMY, 1,
     $                          ZERO, WORK, NS )
*
              DO 90 J = 1, NDIM, NB
                NLENG  = NDIM - J - KB + 1
                NXTROW = MOD( NPROW+ICURROW-1, NPROW )
                NXTCOL = MOD( NPCOL+ICURCOL-1, NPCOL )
                IF( MYCOL.EQ.ICURCOL ) JN = JJ - NB
*
                IF( MYROW.EQ.ICURROW ) THEN
*
*                 Receive updated blocks from previous row of processes
*
                  JJN = JJ * NS + 1
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(NLENG+JB, NCOMM), NB, ICURCOL,
     $                              MYCOL, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART, WORK(IPT), NS,
     $                            MBTROW, MYCOL )
                    CALL PBDMATADD( ICONTXT, 'G', NS, NPART, ONE,
     $                              WORK(IPT), NS, ONE,
     $                              WORK((JJ+NB-NPART)*NS+1), NS )
                  END IF
*
*                 B(1,JJ+1) <== A(II,JJ+1) / B(1,JJ+1)
*                                          ( B(1,JJ+1) = WORK(JJ*NS+1) )
*                 where A(II,JJ+1) is a lower triangular matrix
*
                  IF( MYCOL.EQ.ICURCOL ) THEN
                    CALL DTRSM( 'Right', 'Lower', 'No', DIAG, NS, JB,
     $                          ONE, A(II,JJ+1),LDA, WORK(JJ*NS+1),NS )
                    CALL PBDMATADD( ICONTXT, 'G', NS, JB, ONE,
     $                              WORK(JJ*NS+1), NS, ZERO, WORK(IPT),
     $                              NS )
                    CALL DGEBS2D( ICONTXT, 'Row', 'D-ring', NS, JB,
     $                            WORK(IPT), NS )
                  ELSE
                    CALL DGEBR2D( ICONTXT, 'Row', 'D-ring', NS, JB,
     $                            WORK(IPT), NS, MYROW, ICURCOL )
                  END IF
*
*                 Update the rest of data and prepare for the next step
*
                  IF( NLENG.GT.0 ) THEN
*
*                   Update the (NPROW-1) blocks first
*
                    NREST = MIN( NLENG, NCOMM )
                    NPART = NUMROC( NREST, NB, ICURCOL, MYCOL+1, NPCOL )
                    JJN   = JN + NB - NPART
                    CALL DGEMM( 'No', 'No', NS, NPART, JB, -ONE,
     $                          WORK(IPT), NS, A(II,JJN+1), LDA, ONE,
     $                          WORK(JJN*NS+1), NS )
*
*                   Send updated blocks to next column of processes
*
                    IF( NPROW.GT.1 )
     $                CALL DGESD2D( ICONTXT, NS, NPART, WORK(JJN*NS+1),
     $                              NS, MTPROW, MYCOL )
*
*                   Update the rest of the matrix
*
                    CALL DGEMM( 'No','No', NS, JJN, JB, -ONE, WORK(IPT),
     $                          NS, A(II,1), LDA, ONE, WORK, NS )
                  END IF
*
*                 Send the solution blocks to destination (IDEST row)
*
                  JJN = J + KB - 1
                  IF( J.GT.1 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRPB, JJN-JB), NB, MYCOL,
     $                              ICURCOL+1, NPCOL )
                    CALL DGERV2D( ICONTXT, NS, NPART,
     $                            WORK((JJ+NB)*NS+1), NS,
     $                            MBTROW, MYCOL )
                  END IF
*
                  IF( NLENG.GT.0 .AND. NPROW.GT.1 ) THEN
                    NPART = NUMROC( MIN(IRDB, JJN), NB, MYCOL, ICURCOL,
     $                              NPCOL )
                    CALL DGESD2D( ICONTXT, NS, NPART,
     $                            WORK((JN+NB)*NS+1), NS,
     $                            MTPROW, MYCOL )
                  END IF
*
                  II = II - NB
                END IF
*
                JJ = JN
                JB = NB
                ICURROW = NXTROW
                ICURCOL = NXTCOL
   90         CONTINUE
*
*             Uncopied solutions are moved to the first column of procs.
*
              ICURROW = MOD( ICURROW+1, NPROW )
              IF( ICURROW.NE.IDEST ) THEN
                KDIST = MOD( NPROW+IDEST-ICURROW, NPROW )
                NPART = NUMROC( MIN(NDIM, KDIST*NB), NB, MYCOL, IACOL,
     $                          NPCOL )
                IF( MYROW.EQ.ICURROW ) THEN
                  CALL DGESD2D( ICONTXT, NS, NPART, WORK, NS,
     $                          IDEST, MYCOL )
                ELSE IF( MYROW.EQ.IDEST ) THEN
                  CALL DGERV2D( ICONTXT, NS, NPART, WORK, NS,
     $                          ICURROW, MYCOL )
                END IF
              END IF
            END IF
*
            IF( LSIDE ) THEN
              CALL PBDTRAN( ICONTXT, 'Row', TRANSA, NS, NDIM, NB, WORK,
     $                      NS, ZERO, B, LDB, IDEST, IACOL, IAROW,
     $                      IBPOS, WORK(IPT) )
            ELSE
              IF( .NOT.BDATA .AND. MYROW.EQ.IDEST )
     $          CALL PBDMATADD( ICONTXT, 'G', NS, MQ, ONE, WORK, NS,
     $                          ZERO, B, LDB )
            END IF
          END IF
        END IF
*
* === If A is just a block ===
*
      ELSE
        ADATA = .FALSE.
        ASPACE = LSAME( ABWORK, 'Y' )
        COMMA = ACOMM
        IF( LSAME( COMMA, ' ' ) ) COMMA = '1'
*
        IF( LSIDE .AND. MYROW.EQ.IAROW ) THEN
*
*         Form  B := op( A ) \ alpha * B.
*            _____________       _       _____________
*           |______B______|  =  |_|  \  |______B______|
*                              op(A)
*
          IF( IACOL.EQ.-1 )  ADATA = .TRUE.
          NQ = NUMROC( N, NB, MYCOL, IBPOS, NPCOL )
*
          IF( LDA.LT.MAX(1,M) .AND. ( ASPACE .OR.
     $             IACOL.EQ.MYCOL .OR. IACOL.EQ.-1 ) ) THEN
            INFO = 12
          ELSE IF( LDB.LT. MAX(1,M)                  ) THEN
            INFO = 14
          ELSE IF( IAROW.LT. 0 .OR. IAROW.GE.NPROW   ) THEN
            INFO = 15
          ELSE IF( IACOL.LT.-1 .OR. IACOL.GE.NPCOL   ) THEN
            INFO = 16
          ELSE IF( IBPOS.LT. 0 .OR. IBPOS.GE.NPCOL   ) THEN
            INFO = 17
          END IF
          IF( INFO.NE.0 ) GO TO 10
*
*         Broadcast A if necessary
*
          IF( .NOT. ADATA ) THEN
            IF( ASPACE ) THEN
              IF( MYCOL.EQ.IACOL ) THEN
                CALL DTRBS2D( ICONTXT, 'Row', COMMA, UPLO, DIAG, M, M,
     $                        A, LDA )
              ELSE
                CALL DTRBR2D( ICONTXT, 'Row', COMMA, UPLO, DIAG, M, M,
     $                        A, LDA, MYROW, IACOL )
              END IF
              ADATA = .TRUE.
            ELSE
              IF( MYCOL.EQ.IACOL ) THEN
                CALL DTRBS2D( ICONTXT, 'Row', COMMA, UPLO, DIAG, M, M,
     $                        A, LDA )
                CALL PBDMATADD( ICONTXT, UPLO, M, M, ONE, A, LDA, ZERO,
     $                          WORK, M )
              ELSE
                CALL DTRBR2D( ICONTXT, 'Row', COMMA, UPLO, DIAG, M, M,
     $                        WORK, M, MYROW, IACOL )
              END IF
            END IF
          END IF
*
*         Compute DTRSM
*
          IF( ADATA ) THEN
            CALL DTRSM( 'Left', UPLO, TRANSA, DIAG, M, NQ, ALPHA,
     $                  A, LDA, B, LDB )
          ELSE
            CALL DTRSM( 'Left', UPLO, TRANSA, DIAG, M, NQ, ALPHA,
     $                  WORK, M, B, LDB )
          END IF
*
        ELSE IF( LSAME( SIDE, 'R' ) .AND. MYCOL.EQ.IACOL ) THEN
*
*         Form  B := alpha*B / op( A ).
*                _         _
*               | |       | |
*               | |       | |
*               | |       | |        _
*               |B|    =  |B|   /   |_|
*               | |       | |       op(A)
*               | |       | |
*               |_|       |_|
*
          IF( IAROW.EQ.-1 )  ADATA = .TRUE.
          MP = NUMROC( M, NB, MYROW, IBPOS, NPROW )
*
          IF( LDA.LT.MAX(1,N) .AND. ( ASPACE .OR.
     $             IAROW.EQ.MYROW .OR. IAROW.EQ.-1 ) ) THEN
            INFO = 12
          ELSE IF( LDB .LT.MAX(1,MP)                 ) THEN
            INFO = 14
          ELSE IF( IAROW.LT.-1 .OR. IAROW.GE.NPROW   ) THEN
            INFO = 15
          ELSE IF( IACOL.LT. 0 .OR. IACOL.GE.NPCOL   ) THEN
            INFO = 16
          ELSE IF( IBPOS.LT. 0 .OR. IBPOS.GE.NPROW   ) THEN
            INFO = 17
          END IF
          IF( INFO.NE.0 ) GO TO 10
*
*         Broadcast B if necessary
*
          IF( .NOT. ADATA ) THEN
            IF( ASPACE ) THEN
              IF ( MYROW.EQ.IAROW ) THEN
                CALL DTRBS2D( ICONTXT, 'Col', COMMA, UPLO, DIAG, N, N,
     $                        A, LDA )
              ELSE
                CALL DTRBR2D( ICONTXT, 'Col', COMMA, UPLO, DIAG, N, N,
     $                        A, LDA, IAROW, MYCOL )
              END IF
              ADATA = .TRUE.
            ELSE
              IF ( MYROW.EQ.IAROW ) THEN
                CALL DTRBS2D( ICONTXT, 'Col', COMMA, UPLO, DIAG, N, N,
     $                        A, LDA )
                CALL PBDMATADD( ICONTXT, UPLO, N, N, ONE, A, LDA, ZERO,
     $                          WORK, N )
              ELSE
                CALL DTRBR2D( ICONTXT, 'Col', COMMA, UPLO, DIAG, N, N,
     $                        WORK, N, IAROW, MYCOL )
              END IF
            END IF
          END IF
*
*         Compute DTRSM
*
          IF( ADATA ) THEN
            CALL DTRSM( 'Right', UPLO, TRANSA, DIAG, MP, N, ALPHA,
     $                  A, LDA, B, LDB )
          ELSE
            CALL DTRSM( 'Right', UPLO, TRANSA, DIAG, MP, N, ALPHA,
     $                  WORK, N, B, LDB )
          END IF
        END IF
      END IF
*
      RETURN
*
*     End of PBDTRSM
*
      END