pdblastst.f

Go to the documentation of this file.

      SUBROUTINE pdoptee( ICTXT, NOUT, SUBPTR, SCODE, SNAME )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)      SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL           subptr
*     ..
*
*  Purpose
*  =======
*
*  PDOPTEE  tests  whether  the  PBLAS respond correctly to a bad option
*  argument.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  Calling sequence encodings
*  ==========================
*
*  code Formal argument list                                Examples
*
*  11   (n,      v1,v2)                                     _SWAP, _COPY
*  12   (n,s1,   v1   )                                     _SCAL, _SCAL
*  13   (n,s1,   v1,v2)                                     _AXPY, _DOT_
*  14   (n,s1,i1,v1   )                                     _AMAX
*  15   (n,u1,   v1   )                                     _ASUM, _NRM2
*
*  21   (     trans,     m,n,s1,m1,v1,s2,v2)                _GEMV
*  22   (uplo,             n,s1,m1,v1,s2,v2)                _SYMV, _HEMV
*  23   (uplo,trans,diag,  n,   m1,v1      )                _TRMV, _TRSV
*  24   (                m,n,s1,v1,v2,m1)                   _GER_
*  25   (uplo,             n,s1,v1,   m1)                   _SYR
*  26   (uplo,             n,u1,v1,   m1)                   _HER
*  27   (uplo,             n,s1,v1,v2,m1)                   _SYR2, _HER2
*
*  31   (          transa,transb,     m,n,k,s1,m1,m2,s2,m3) _GEMM
*  32   (side,uplo,                   m,n,  s1,m1,m2,s2,m3) _SYMM, _HEMM
*  33   (     uplo,trans,               n,k,s1,m1,   s2,m3) _SYRK
*  34   (     uplo,trans,               n,k,u1,m1,   u2,m3) _HERK
*  35   (     uplo,trans,               n,k,s1,m1,m2,s2,m3) _SYR2K
*  36   (     uplo,trans,               n,k,s1,m1,m2,u2,m3) _HER2K
*  37   (                             m,n,  s1,m1,   s2,m3) _TRAN_
*  38   (side,uplo,transa,       diag,m,n,  s1,m1,m2      ) _TRMM, _TRSM
*  39   (          trans,             m,n,  s1,m1,   s2,m3) _GEADD
*  40   (     uplo,trans,             m,n,  s1,m1,   s2,m3) _TRADD
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER             APOS
*     ..
*     .. External Subroutines ..
      EXTERNAL            pdchkopt
*     ..
*     .. Executable Statements ..
*
*     Level 2 PBLAS
*
      IF( scode.EQ.21 ) THEN
*
*        Check 1st (and only) option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'A', apos )
*
      ELSE IF( scode.EQ.22 .OR. scode.EQ.25 .OR. scode.EQ.26 .OR.
     $         scode.EQ.27 ) THEN
*
*        Check 1st (and only) option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'U', apos )
*
      ELSE IF( scode.EQ.23 ) THEN
*
*        Check 1st option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'U', apos )
*
*        Check 2nd option
*
         apos = 2
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 3rd option
*
         apos = 3
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'D', apos )
*
*     Level 3 PBLAS
*
      ELSE IF( scode.EQ.31 ) THEN
*
*        Check 1st option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2'nd option
*
         apos = 2
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'B', apos )
*
      ELSE IF( scode.EQ.32 ) THEN
*
*        Check 1st option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'S', apos )
*
*        Check 2nd option
*
         apos = 2
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'U', apos )
*
      ELSE IF( scode.EQ.33 .OR. scode.EQ.34 .OR. scode.EQ.35 .OR.
     $         scode.EQ.36 .OR. scode.EQ.40 ) THEN
*
*        Check 1st option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'U', apos )
*
*        Check 2'nd option
*
         apos = 2
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'A', apos )
*
      ELSE IF( scode.EQ.38 ) THEN
*
*        Check 1st option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'S', apos )
*
*        Check 2nd option
*
         apos = 2
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'U', apos )
*
*        Check 3rd option
*
         apos = 3
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 4th option
*
         apos = 4
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'D', apos )
*
*
      ELSE IF( scode.EQ.39 ) THEN
*
*        Check 1st option
*
         apos = 1
         CALL pdchkopt( ictxt, nout, subptr, scode, sname, 'A', apos )
*
      END IF
*
      RETURN
*
*     End of PDOPTEE
*
      END
      SUBROUTINE pdchkopt( ICTXT, NOUT, SUBPTR, SCODE, SNAME, ARGNAM,
     $                     ARGPOS )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1         ARGNAM
      INTEGER             ARGPOS, ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)       SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL            subptr
*     ..
*
*  Purpose
*  =======
*
*  PDCHKOPT tests the option ARGNAM in any PBLAS routine.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  ARGNAM  (global input) CHARACTER*(*)
*          On entry,  ARGNAM  specifies  the  name  of  the option to be
*          checked. ARGNAM can either be 'D', 'S', 'A', 'B', or 'U'.
*
*  ARGPOS  (global input) INTEGER
*          On entry, ARGPOS indicates the position of the option ARGNAM
*          to be tested.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            INFOT
*     ..
*     .. External Subroutines ..
      EXTERNAL           pchkpbe, pdcallsub, pdsetpblas
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. Common Blocks ..
      CHARACTER          DIAG, SIDE, TRANSA, TRANSB, UPLO
      COMMON             /pblasc/diag, side, transa, transb, uplo
*     ..
*     .. Executable Statements ..
*
*     Reiniatilize the dummy arguments to correct values
*
      CALL pdsetpblas( ictxt )
*
      IF( lsame( argnam, 'D' ) ) THEN
*
*        Generate bad DIAG option
*
         diag = '/'
*
      ELSE IF( lsame( argnam, 'S' ) ) THEN
*
*        Generate bad SIDE option
*
         side = '/'
*
      ELSE IF( lsame( argnam, 'A' ) ) THEN
*
*        Generate bad TRANSA option
*
         transa = '/'
*
      ELSE IF( lsame( argnam, 'B' ) ) THEN
*
*        Generate bad TRANSB option
*
         transb = '/'
*
      ELSE IF( lsame( argnam, 'U' ) ) THEN
*
*        Generate bad UPLO option
*
         uplo = '/'
*
      END IF
*
*     Set INFOT to the position of the bad dimension argument
*
      infot = argpos
*
*     Call the PBLAS routine
*
      CALL pdcallsub( subptr, scode )
      CALL pchkpbe( ictxt, nout, sname, infot )
*
      RETURN
*
*     End of PDCHKOPT
*
      END
      SUBROUTINE pddimee( ICTXT, NOUT, SUBPTR, SCODE, SNAME )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)      SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL           subptr
*     ..
*
*  Purpose
*  =======
*
*  PDDIMEE  tests whether the PBLAS respond correctly to a bad dimension
*  argument.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  Calling sequence encodings
*  ==========================
*
*  code Formal argument list                                Examples
*
*  11   (n,      v1,v2)                                     _SWAP, _COPY
*  12   (n,s1,   v1   )                                     _SCAL, _SCAL
*  13   (n,s1,   v1,v2)                                     _AXPY, _DOT_
*  14   (n,s1,i1,v1   )                                     _AMAX
*  15   (n,u1,   v1   )                                     _ASUM, _NRM2
*
*  21   (     trans,     m,n,s1,m1,v1,s2,v2)                _GEMV
*  22   (uplo,             n,s1,m1,v1,s2,v2)                _SYMV, _HEMV
*  23   (uplo,trans,diag,  n,   m1,v1      )                _TRMV, _TRSV
*  24   (                m,n,s1,v1,v2,m1)                   _GER_
*  25   (uplo,             n,s1,v1,   m1)                   _SYR
*  26   (uplo,             n,u1,v1,   m1)                   _HER
*  27   (uplo,             n,s1,v1,v2,m1)                   _SYR2, _HER2
*
*  31   (          transa,transb,     m,n,k,s1,m1,m2,s2,m3) _GEMM
*  32   (side,uplo,                   m,n,  s1,m1,m2,s2,m3) _SYMM, _HEMM
*  33   (     uplo,trans,               n,k,s1,m1,   s2,m3) _SYRK
*  34   (     uplo,trans,               n,k,u1,m1,   u2,m3) _HERK
*  35   (     uplo,trans,               n,k,s1,m1,m2,s2,m3) _SYR2K
*  36   (     uplo,trans,               n,k,s1,m1,m2,u2,m3) _HER2K
*  37   (                             m,n,  s1,m1,   s2,m3) _TRAN_
*  38   (side,uplo,transa,       diag,m,n,  s1,m1,m2      ) _TRMM, _TRSM
*  39   (          trans,             m,n,  s1,m1,   s2,m3) _GEADD
*  40   (     uplo,trans,             m,n,  s1,m1,   s2,m3) _TRADD
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER             APOS
*     ..
*     .. External Subroutines ..
      EXTERNAL            pdchkdim
*     ..
*     .. Executable Statements ..
*
*     Level 1 PBLAS
*
      IF( scode.EQ.11 .OR. scode.EQ.12 .OR. scode.EQ.13 .OR.
     $    scode.EQ.14 .OR. scode.EQ.15 ) THEN
*
*        Check 1st (and only) dimension
*
         apos = 1
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
*     Level 2 PBLAS
*
      ELSE IF( scode.EQ.21 ) THEN
*
*        Check 1st dimension
*
         apos = 2
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 3
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.22 .OR. scode.EQ.25 .OR. scode.EQ.26 .OR.
     $         scode.EQ.27 ) THEN
*
*        Check 1st (and only) dimension
*
         apos = 2
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.23 ) THEN
*
*        Check 1st (and only) dimension
*
         apos = 4
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.24 ) THEN
*
*        Check 1st dimension
*
         apos = 1
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 2
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
*     Level 3 PBLAS
*
      ELSE IF( scode.EQ.31 ) THEN
*
*        Check 1st dimension
*
         apos = 3
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 4
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
*        Check 3rd dimension
*
         apos = 5
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'K', apos )
*
      ELSE IF( scode.EQ.32 ) THEN
*
*        Check 1st dimension
*
         apos = 3
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 4
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.33 .OR. scode.EQ.34 .OR. scode.EQ.35 .OR.
     $         scode.EQ.36 ) THEN
*
*        Check 1st dimension
*
         apos = 3
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
*        Check 2nd dimension
*
         apos = 4
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'K', apos )
*
      ELSE IF( scode.EQ.37 ) THEN
*
*        Check 1st dimension
*
         apos = 1
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 2
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.38 ) THEN
*
*        Check 1st dimension
*
         apos = 5
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 6
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.39 ) THEN
*
*        Check 1st dimension
*
         apos = 2
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 3
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      ELSE IF( scode.EQ.40 ) THEN
*
*        Check 1st dimension
*
         apos = 3
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'M', apos )
*
*        Check 2nd dimension
*
         apos = 4
         CALL pdchkdim( ictxt, nout, subptr, scode, sname, 'N', apos )
*
      END IF
*
      RETURN
*
*     End of PDDIMEE
*
      END
      SUBROUTINE pdchkdim( ICTXT, NOUT, SUBPTR, SCODE, SNAME, ARGNAM,
     $                     ARGPOS )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1         ARGNAM
      INTEGER             ARGPOS, ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)       SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL            subptr
*     ..
*
*  Purpose
*  =======
*
*  PDCHKDIM tests the dimension ARGNAM in any PBLAS routine.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  ARGNAM  (global input) CHARACTER*(*)
*          On entry,  ARGNAM  specifies  the name of the dimension to be
*          checked. ARGNAM can either be 'M', 'N' or 'K'.
*
*  ARGPOS  (global input) INTEGER
*          On entry, ARGPOS indicates the position of the option ARGNAM
*          to be tested.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            INFOT
*     ..
*     .. External Subroutines ..
      EXTERNAL           pchkpbe, pdcallsub, pdsetpblas
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. Common Blocks ..
      INTEGER            KDIM, MDIM, NDIM
      COMMON             /PBLASN/KDIM, MDIM, NDIM
*     ..
*     .. Executable Statements ..
*
*     Reiniatilize the dummy arguments to correct values
*
      CALL pdsetpblas( ictxt )
*
      IF( lsame( argnam, 'M' ) ) THEN
*
*        Generate bad MDIM
*
         mdim = -1
*
      ELSE IF( lsame( argnam, 'N' ) ) THEN
*
*        Generate bad NDIM
*
         ndim = -1
*
      ELSE
*
*        Generate bad KDIM
*
         kdim = -1
*
      END IF
*
*     Set INFOT to the position of the bad dimension argument
*
      infot = argpos
*
*     Call the PBLAS routine
*
      CALL pdcallsub( subptr, scode )
      CALL pchkpbe( ictxt, nout, sname, infot )
*
      RETURN
*
*     End of PDCHKDIM
*
      END
      SUBROUTINE pdvecee( ICTXT, NOUT, SUBPTR, SCODE, SNAME )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER             ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*7         SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL            subptr
*     ..
*
*  Purpose
*  =======
*
*  PDVECEE  tests  whether  the  PBLAS respond correctly to a bad vector
*  argument.  Each  vector <vec> is described by: <vec>, I<vec>, J<vec>,
*  DESC<vec>,  INC<vec>.   Out   of  all  these,  only  I<vec>,  J<vec>,
*  DESC<vec>, and INC<vec> can be tested.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  Calling sequence encodings
*  ==========================
*
*  code Formal argument list                                Examples
*
*  11   (n,      v1,v2)                                     _SWAP, _COPY
*  12   (n,s1,   v1   )                                     _SCAL, _SCAL
*  13   (n,s1,   v1,v2)                                     _AXPY, _DOT_
*  14   (n,s1,i1,v1   )                                     _AMAX
*  15   (n,u1,   v1   )                                     _ASUM, _NRM2
*
*  21   (     trans,     m,n,s1,m1,v1,s2,v2)                _GEMV
*  22   (uplo,             n,s1,m1,v1,s2,v2)                _SYMV, _HEMV
*  23   (uplo,trans,diag,  n,   m1,v1      )                _TRMV, _TRSV
*  24   (                m,n,s1,v1,v2,m1)                   _GER_
*  25   (uplo,             n,s1,v1,   m1)                   _SYR
*  26   (uplo,             n,u1,v1,   m1)                   _HER
*  27   (uplo,             n,s1,v1,v2,m1)                   _SYR2, _HER2
*
*  31   (          transa,transb,     m,n,k,s1,m1,m2,s2,m3) _GEMM
*  32   (side,uplo,                   m,n,  s1,m1,m2,s2,m3) _SYMM, _HEMM
*  33   (     uplo,trans,               n,k,s1,m1,   s2,m3) _SYRK
*  34   (     uplo,trans,               n,k,u1,m1,   u2,m3) _HERK
*  35   (     uplo,trans,               n,k,s1,m1,m2,s2,m3) _SYR2K
*  36   (     uplo,trans,               n,k,s1,m1,m2,u2,m3) _HER2K
*  37   (                             m,n,  s1,m1,   s2,m3) _TRAN_
*  38   (side,uplo,transa,       diag,m,n,  s1,m1,m2      ) _TRMM, _TRSM
*  39   (          trans,             m,n,  s1,m1,   s2,m3) _GEADD
*  40   (     uplo,trans,             m,n,  s1,m1,   s2,m3) _TRADD
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER             APOS
*     ..
*     .. External Subroutines ..
      EXTERNAL            pdchkmat
*     ..
*     .. Executable Statements ..
*
*     Level 1 PBLAS
*
      IF( scode.EQ.11 ) THEN
*
*        Check 1st vector
*
         apos = 2
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
*        Check 2nd vector
*
         apos = 7
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'Y', apos )
*
      ELSE IF( scode.EQ.12 .OR. scode.EQ.15 ) THEN
*
*        Check 1st (and only) vector
*
         apos = 3
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
      ELSE IF( scode.EQ.13 ) THEN
*
*        Check 1st vector
*
         apos = 3
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
*        Check 2nd vector
*
         apos = 8
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'Y', apos )
*
      ELSE IF( scode.EQ.14 ) THEN
*
*        Check 1st (and only) vector
*
         apos = 4
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
*     Level 2 PBLAS
*
      ELSE IF( scode.EQ.21 ) THEN
*
*        Check 1st vector
*
         apos = 9
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
*        Check 2nd vector
*
         apos = 15
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'Y', apos )
*
      ELSE IF( scode.EQ.22 ) THEN
*
*        Check 1st vector
*
         apos = 8
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
*        Check 2nd vector
*
         apos = 14
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'Y', apos )
*
      ELSE IF( scode.EQ.23 ) THEN
*
*        Check 1st (and only) vector
*
         apos = 9
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
      ELSE IF( scode.EQ.24 .OR. scode.EQ.27 ) THEN
*
*        Check 1st vector
*
         apos = 4
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
*        Check 2nd vector
*
         apos = 9
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'Y', apos )
*
      ELSE IF( scode.EQ.26 .OR. scode.EQ.27 ) THEN
*
*        Check 1'st (and only) vector
*
         apos = 4
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'X', apos )
*
      END IF
*
      RETURN
*
*     End of PDVECEE
*
      END
      SUBROUTINE pdmatee( ICTXT, NOUT, SUBPTR, SCODE, SNAME )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER             ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*7         SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL            subptr
*     ..
*
*  Purpose
*  =======
*
*  PDMATEE  tests  whether  the  PBLAS respond correctly to a bad matrix
*  argument.  Each  matrix <mat> is described by: <mat>, I<mat>, J<mat>,
*  and DESC<mat>.  Out  of  all these, only I<vec>, J<vec> and DESC<mat>
*  can be tested.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  Calling sequence encodings
*  ==========================
*
*  code Formal argument list                                Examples
*
*  11   (n,      v1,v2)                                     _SWAP, _COPY
*  12   (n,s1,   v1   )                                     _SCAL, _SCAL
*  13   (n,s1,   v1,v2)                                     _AXPY, _DOT_
*  14   (n,s1,i1,v1   )                                     _AMAX
*  15   (n,u1,   v1   )                                     _ASUM, _NRM2
*
*  21   (     trans,     m,n,s1,m1,v1,s2,v2)                _GEMV
*  22   (uplo,             n,s1,m1,v1,s2,v2)                _SYMV, _HEMV
*  23   (uplo,trans,diag,  n,   m1,v1      )                _TRMV, _TRSV
*  24   (                m,n,s1,v1,v2,m1)                   _GER_
*  25   (uplo,             n,s1,v1,   m1)                   _SYR
*  26   (uplo,             n,u1,v1,   m1)                   _HER
*  27   (uplo,             n,s1,v1,v2,m1)                   _SYR2, _HER2
*
*  31   (          transa,transb,     m,n,k,s1,m1,m2,s2,m3) _GEMM
*  32   (side,uplo,                   m,n,  s1,m1,m2,s2,m3) _SYMM, _HEMM
*  33   (     uplo,trans,               n,k,s1,m1,   s2,m3) _SYRK
*  34   (     uplo,trans,               n,k,u1,m1,   u2,m3) _HERK
*  35   (     uplo,trans,               n,k,s1,m1,m2,s2,m3) _SYR2K
*  36   (     uplo,trans,               n,k,s1,m1,m2,u2,m3) _HER2K
*  37   (                             m,n,  s1,m1,   s2,m3) _TRAN_
*  38   (side,uplo,transa,       diag,m,n,  s1,m1,m2      ) _TRMM, _TRSM
*  39   (          trans,             m,n,  s1,m1,   s2,m3) _GEADD
*  40   (     uplo,trans,             m,n,  s1,m1,   s2,m3) _TRADD
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER             APOS
*     ..
*     .. External Subroutines ..
      EXTERNAL            pdchkmat
*     ..
*     .. Executable Statements ..
*
*     Level 2 PBLAS
*
      IF( scode.EQ.21 .OR. scode.EQ.23 ) THEN
*
*        Check 1st (and only) matrix
*
         apos = 5
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
      ELSE IF( scode.EQ.22 ) THEN
*
*        Check 1st (and only) matrix
*
         apos = 4
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
      ELSE IF( scode.EQ.24 .OR. scode.EQ.27 ) THEN
*
*        Check 1st (and only) matrix
*
         apos = 14
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
      ELSE IF( scode.EQ.25 .OR. scode.EQ.26 ) THEN
*
*        Check 1st (and only) matrix
*
         apos = 9
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*     Level 3 PBLAS
*
      ELSE IF( scode.EQ.31 ) THEN
*
*        Check 1st matrix
*
         apos = 7
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 11
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'B', apos )
*
*        Check 3nd matrix
*
         apos = 16
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'C', apos )
*
      ELSE IF( scode.EQ.32 .OR. scode.EQ.35 .OR. scode.EQ.36 ) THEN
*
*        Check 1st matrix
*
         apos = 6
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 10
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'B', apos )
*
*        Check 3nd matrix
*
         apos = 15
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'C', apos )
*
      ELSE IF( scode.EQ.33 .OR. scode.EQ.34 ) THEN
*
*        Check 1st matrix
*
         apos = 6
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 11
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'C', apos )
*
      ELSE IF( scode.EQ.37 ) THEN
*
*        Check 1st matrix
*
         apos = 4
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 9
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'C', apos )
*
      ELSE IF( scode.EQ.38 ) THEN
*
*        Check 1st matrix
*
         apos = 8
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 12
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'B', apos )
*
      ELSE IF( scode.EQ.39 ) THEN
*
*        Check 1st matrix
*
         apos = 5
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 10
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'C', apos )
*
      ELSE IF( scode.EQ.40 ) THEN
*
*        Check 1st matrix
*
         apos = 6
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'A', apos )
*
*        Check 2nd matrix
*
         apos = 11
         CALL pdchkmat( ictxt, nout, subptr, scode, sname, 'C', apos )
*
      END IF
*
      RETURN
*
*     End of PDMATEE
*
      END
      SUBROUTINE pdsetpblas( ICTXT )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            ICTXT
*     ..
*
*  Purpose
*  =======
*
*  PDSETPBLAS initializes *all* the dummy arguments to correct values.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   rsrc_
      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0d+0 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           pb_descset2
*     ..
*     .. Common Blocks ..
      CHARACTER*1        DIAG, SIDE, TRANSA, TRANSB, UPLO
      INTEGER            IA, IB, IC, INCX, INCY, ISCLR, IX, IY, JA, JB,
     $                   jc, jx, jy, kdim, mdim, ndim
      DOUBLE PRECISION   USCLR, SCLR
      INTEGER            DESCA( DLEN_ ), DESCB( DLEN_ ), DESCC( DLEN_ ),
     $                   descx( dlen_ ), descy( dlen_ )
      DOUBLE PRECISION   A( 2, 2 ), B( 2, 2 ), C( 2, 2 ), X( 2 ), Y( 2 )
      COMMON             /PBLASC/DIAG, SIDE, TRANSA, TRANSB, UPLO
      COMMON             /pblasd/desca, descb, descc, descx, descy
      COMMON             /pblasi/ia, ib, ic, incx, incy, isclr, ix, iy,
     $                   ja, jb, jc, jx, jy
      COMMON             /pblasm/a, b, c
      COMMON             /pblasn/kdim, mdim, ndim
      COMMON             /pblass/sclr, usclr
      COMMON             /pblasv/x, y
*     ..
*     .. Executable Statements ..
*
*     Set default values for options
*
      diag   = 'N'
      side   = 'L'
      transa = 'N'
      transb = 'N'
      uplo   = 'U'
*
*     Set default values for scalars
*
      kdim   = 1
      mdim   = 1
      ndim   = 1
      isclr  = 1
      sclr   = one
      usclr  = one
*
*     Set default values for distributed matrix A
*
      a( 1, 1 ) = one
      a( 2, 1 ) = one
      a( 1, 2 ) = one
      a( 2, 2 ) = one
      ia = 1
      ja = 1
      CALL pb_descset2( desca, 2, 2, 1, 1, 1, 1, 0, 0, ictxt, 2 )
*
*     Set default values for distributed matrix B
*
      b( 1, 1 ) = one
      b( 2, 1 ) = one
      b( 1, 2 ) = one
      b( 2, 2 ) = one
      ib = 1
      jb = 1
      CALL pb_descset2( descb, 2, 2, 1, 1, 1, 1, 0, 0, ictxt, 2 )
*
*     Set default values for distributed matrix C
*
      c( 1, 1 ) = one
      c( 2, 1 ) = one
      c( 1, 2 ) = one
      c( 2, 2 ) = one
      ic = 1
      jc = 1
      CALL pb_descset2( descc, 2, 2, 1, 1, 1, 1, 0, 0, ictxt, 2 )
*
*     Set default values for distributed matrix X
*
      x( 1 ) = one
      x( 2 ) = one
      ix = 1
      jx = 1
      CALL pb_descset2( descx, 2, 1, 1, 1, 1, 1, 0, 0, ictxt, 2 )
      incx = 1
*
*     Set default values for distributed matrix Y
*
      y( 1 ) = one
      y( 2 ) = one
      iy = 1
      jy = 1
      CALL pb_descset2( descy, 2, 1, 1, 1, 1, 1, 0, 0, ictxt, 2 )
      incy = 1
*
      RETURN
*
*     End of PDSETPBLAS
*
      END
      SUBROUTINE pdchkmat( ICTXT, NOUT, SUBPTR, SCODE, SNAME, ARGNAM,
     $                     ARGPOS )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1         ARGNAM
      INTEGER             ARGPOS, ICTXT, NOUT, SCODE
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)       SNAME
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL            subptr
*     ..
*
*  Purpose
*  =======
*
*  PDCHKMAT tests the matrix (or vector) ARGNAM in any PBLAS routine.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  SNAME   (global input) CHARACTER*(*)
*          On entry,  SNAME  specifies  the subroutine name calling this
*          subprogram.
*
*  ARGNAM  (global input) CHARACTER*(*)
*          On entry,  ARGNAM  specifies the name of the matrix or vector
*          to be checked.  ARGNAM can either be 'A', 'B' or 'C' when one
*          wants to check a matrix, and 'X' or 'Y' for a vector.
*
*  ARGPOS  (global input) INTEGER
*          On entry, ARGPOS indicates the position of the first argument
*          of the matrix (or vector) ARGNAM.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      INTEGER             DESCMULT
      PARAMETER           ( DESCMULT = 100 )
*     ..
*     .. Local Scalars ..
      INTEGER             I, INFOT, NPROW, NPCOL, MYROW, MYCOL
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, pchkpbe, pdcallsub, pdsetpblas
*     ..
*     .. External Functions ..
      LOGICAL             LSAME
      EXTERNAL            LSAME
*     ..
*     .. Common Blocks ..
      INTEGER            IA, IB, IC, INCX, INCY, ISCLR, IX, IY, JA, JB,
     $                   JC, JX, JY
      INTEGER            DESCA( DLEN_ ), DESCB( DLEN_ ), DESCC( DLEN_ ),
     $                   descx( dlen_ ), descy( dlen_ )
      COMMON             /pblasd/desca, descb, descc, descx, descy
      COMMON             /pblasi/ia, ib, ic, incx, incy, isclr, ix, iy,
     $                   ja, jb, jc, jx, jy
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      IF( lsame( argnam, 'A' ) ) THEN
*
*        Check IA. Set all other OK, bad IA
*
         CALL pdsetpblas( ictxt )
         ia    = -1
         infot = argpos + 1
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check JA. Set all other OK, bad JA
*
         CALL pdsetpblas( ictxt )
         ja    = -1
         infot = argpos + 2
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check DESCA. Set all other OK, bad DESCA
*
         DO 10 i = 1, dlen_
*
*           Set I'th entry of DESCA to incorrect value, rest ok.
*
            CALL pdsetpblas( ictxt )
            desca( i ) =  -2
            infot = ( ( argpos + 3 ) * descmult ) + i
            CALL pdcallsub( subptr, scode )
            CALL pchkpbe( ictxt, nout, sname, infot )
*
*           Extra tests for RSRCA, CSRCA, LDA
*
            IF( ( i.EQ.rsrc_ ) .OR. ( i.EQ.csrc_ ) .OR.
     $          ( i.EQ.lld_ ) ) THEN
*
               CALL pdsetpblas( ictxt )
*
*              Test RSRCA >= NPROW
*
               IF( i.EQ.rsrc_ )
     $            desca( i ) =  nprow
*
*              Test CSRCA >= NPCOL
*
               IF( i.EQ.csrc_ )
     $            desca( i ) =  npcol
*
*              Test LDA >= MAX(1, PB_NUMROC(...)). Set to 1 as mat 2x2.
*
               IF( i.EQ.lld_ ) THEN
                  IF( myrow.EQ.0 .AND.mycol.EQ.0 ) THEN
                     desca( i ) = 1
                  ELSE
                     desca( i ) = 0
                  END IF
               END IF
*
               infot = ( ( argpos + 3 ) * descmult ) + i
               CALL pdcallsub( subptr, scode )
               CALL pchkpbe( ictxt, nout, sname, infot )
*
            END IF
*
   10    CONTINUE
*
      ELSE IF( lsame( argnam, 'B' ) ) THEN
*
*        Check IB. Set all other OK, bad IB
*
         CALL pdsetpblas( ictxt )
         ib    = -1
         infot = argpos + 1
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check JB. Set all other OK, bad JB
*
         CALL pdsetpblas( ictxt )
         jb    = -1
         infot = argpos + 2
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check DESCB. Set all other OK, bad DESCB
*
         DO 20 i = 1, dlen_
*
*           Set I'th entry of DESCB to incorrect value, rest ok.
*
            CALL pdsetpblas( ictxt )
            descb( i ) =  -2
            infot = ( ( argpos + 3 ) * descmult ) + i
            CALL pdcallsub( subptr, scode )
            CALL pchkpbe( ictxt, nout, sname, infot )
*
*           Extra tests for RSRCB, CSRCB, LDB
*
            IF( ( i.EQ.rsrc_ ) .OR. ( i.EQ.csrc_ ) .OR.
     $          ( i.EQ.lld_ ) ) THEN
*
               CALL pdsetpblas( ictxt )
*
*              Test RSRCB >= NPROW
*
               IF( i.EQ.rsrc_ )
     $            descb( i ) =  nprow
*
*              Test CSRCB >= NPCOL
*
               IF( i.EQ.csrc_ )
     $            descb( i ) =  npcol
*
*              Test LDB >= MAX(1, PB_NUMROC(...)). Set to 1 as mat 2x2.
*
               IF( i.EQ.lld_ ) THEN
                  IF( myrow.EQ.0 .AND.mycol.EQ.0 ) THEN
                     descb( i ) = 1
                  ELSE
                     descb( i ) = 0
                  END IF
               END IF
*
               infot = ( ( argpos + 3 ) * descmult ) + i
               CALL pdcallsub( subptr, scode )
               CALL pchkpbe( ictxt, nout, sname, infot )
*
            END IF
*
   20    CONTINUE
*
      ELSE IF( lsame( argnam, 'C' ) ) THEN
*
*        Check IC. Set all other OK, bad IC
*
         CALL pdsetpblas( ictxt )
         ic    = -1
         infot = argpos + 1
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check JC. Set all other OK, bad JC
*
         CALL pdsetpblas( ictxt )
         jc    = -1
         infot = argpos + 2
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check DESCC. Set all other OK, bad DESCC
*
         DO 30 i = 1, dlen_
*
*           Set I'th entry of DESCC to incorrect value, rest ok.
*
            CALL pdsetpblas( ictxt )
            descc( i ) =  -2
            infot = ( ( argpos + 3 ) * descmult ) + i
            CALL pdcallsub( subptr, scode )
            CALL pchkpbe( ictxt, nout, sname, infot )
*
*           Extra tests for RSRCC, CSRCC, LDC
*
            IF( ( i.EQ.rsrc_ ) .OR. ( i.EQ.csrc_ ) .OR.
     $          ( i.EQ.lld_ ) ) THEN
*
               CALL pdsetpblas( ictxt )
*
*              Test RSRCC >= NPROW
*
               IF( i.EQ.rsrc_ )
     $            descc( i ) =  nprow
*
*              Test CSRCC >= NPCOL
*
               IF( i.EQ.csrc_ )
     $            descc( i ) =  npcol
*
*              Test LDC >= MAX(1, PB_NUMROC(...)). Set to 1 as mat 2x2.
*
               IF( i.EQ.lld_ ) THEN
                  IF( myrow.EQ.0 .AND.mycol.EQ.0 ) THEN
                     descc( i ) = 1
                  ELSE
                     descc( i ) = 0
                  END IF
               END IF
*
               infot = ( ( argpos + 3 ) * descmult ) + i
               CALL pdcallsub( subptr, scode )
               CALL pchkpbe( ictxt, nout, sname, infot )
*
            END IF
*
   30    CONTINUE
*
      ELSE IF( lsame( argnam, 'X' ) ) THEN
*
*        Check IX. Set all other OK, bad IX
*
         CALL pdsetpblas( ictxt )
         ix    = -1
         infot = argpos + 1
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check JX. Set all other OK, bad JX
*
         CALL pdsetpblas( ictxt )
         jx    = -1
         infot = argpos + 2
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check DESCX. Set all other OK, bad DESCX
*
         DO 40 i = 1, dlen_
*
*           Set I'th entry of DESCX to incorrect value, rest ok.
*
            CALL pdsetpblas( ictxt )
            descx( i ) =  -2
            infot = ( ( argpos + 3 ) * descmult ) + i
            CALL pdcallsub( subptr, scode )
            CALL pchkpbe( ictxt, nout, sname, infot )
*
*           Extra tests for RSRCX, CSRCX, LDX
*
            IF( ( i.EQ.rsrc_ ) .OR. ( i.EQ.csrc_ ) .OR.
     $          ( i.EQ.lld_ ) ) THEN
*
               CALL pdsetpblas( ictxt )
*
*              Test RSRCX >= NPROW
*
               IF( i.EQ.rsrc_ )
     $            descx( i ) =  nprow
*
*              Test CSRCX >= NPCOL
*
               IF( i.EQ.csrc_ )
     $            descx( i ) =  npcol
*
*              Test LDX >= MAX(1, PB_NUMROC(...)). Set to 1 as mat 2x2.
*
               IF( i.EQ.lld_ ) THEN
                  IF( myrow.EQ.0 .AND.mycol.EQ.0 ) THEN
                     descx( i ) = 1
                  ELSE
                     descx( i ) = 0
                  END IF
               END IF
*
               infot = ( ( argpos + 3 ) * descmult ) + i
               CALL pdcallsub( subptr, scode )
               CALL pchkpbe( ictxt, nout, sname, infot )
*
            END IF
*
   40    CONTINUE
*
*        Check INCX. Set all other OK, bad INCX
*
         CALL pdsetpblas( ictxt )
         incx  =  -1
         infot = argpos + 4
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
      ELSE
*
*        Check IY. Set all other OK, bad IY
*
         CALL pdsetpblas( ictxt )
         iy    = -1
         infot = argpos + 1
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check JY. Set all other OK, bad JY
*
         CALL pdsetpblas( ictxt )
         jy    = -1
         infot = argpos + 2
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
*        Check DESCY. Set all other OK, bad DESCY
*
         DO 50 i = 1, dlen_
*
*           Set I'th entry of DESCY to incorrect value, rest ok.
*
            CALL pdsetpblas( ictxt )
            descy( i ) =  -2
            infot = ( ( argpos + 3 ) * descmult ) + i
            CALL pdcallsub( subptr, scode )
            CALL pchkpbe( ictxt, nout, sname, infot )
*
*           Extra tests for RSRCY, CSRCY, LDY
*
            IF( ( i.EQ.rsrc_ ) .OR. ( i.EQ.csrc_ ) .OR.
     $          ( i.EQ.lld_ ) ) THEN
*
               CALL pdsetpblas( ictxt )
*
*              Test RSRCY >= NPROW
*
               IF( i.EQ.rsrc_ )
     $            descy( i ) = nprow
*
*              Test CSRCY >= NPCOL
*
               IF( i.EQ.csrc_ )
     $            descy( i ) = npcol
*
*              Test LDY >= MAX(1, PB_NUMROC(...)). Set to 1 as mat 2x2.
*
               IF( i.EQ.lld_ ) THEN
                  IF( myrow.EQ.0 .AND.mycol.EQ.0 ) THEN
                     descy( i ) = 1
                  ELSE
                     descy( i ) = 0
                  END IF
               END IF
*
               infot = ( ( argpos + 3 ) * descmult ) + i
               CALL pdcallsub( subptr, scode )
               CALL pchkpbe( ictxt, nout, sname, infot )
*
            END IF
*
   50    CONTINUE
*
*        Check INCY. Set all other OK, bad INCY
*
         CALL pdsetpblas( ictxt )
         incy =  -1
         infot = argpos + 4
         CALL pdcallsub( subptr, scode )
         CALL pchkpbe( ictxt, nout, sname, infot )
*
      END IF
*
      RETURN
*
*     End of PDCHKMAT
*
      END
      SUBROUTINE pdcallsub( SUBPTR, SCODE )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER             SCODE
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL            subptr
*     ..
*
*  Purpose
*  =======
*
*  PDCALLSUB calls the subroutine SUBPTR with the calling sequence iden-
*  tified by SCODE.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  SUBPTR  (global input) SUBROUTINE
*          On entry,  SUBPTR  is  a  subroutine. SUBPTR must be declared
*          EXTERNAL in the calling subroutine.
*
*  SCODE   (global input) INTEGER
*          On entry, SCODE specifies the calling sequence code.
*
*  Calling sequence encodings
*  ==========================
*
*  code Formal argument list                                Examples
*
*  11   (n,      v1,v2)                                     _SWAP, _COPY
*  12   (n,s1,   v1   )                                     _SCAL, _SCAL
*  13   (n,s1,   v1,v2)                                     _AXPY, _DOT_
*  14   (n,s1,i1,v1   )                                     _AMAX
*  15   (n,u1,   v1   )                                     _ASUM, _NRM2
*
*  21   (     trans,     m,n,s1,m1,v1,s2,v2)                _GEMV
*  22   (uplo,             n,s1,m1,v1,s2,v2)                _SYMV, _HEMV
*  23   (uplo,trans,diag,  n,   m1,v1      )                _TRMV, _TRSV
*  24   (                m,n,s1,v1,v2,m1)                   _GER_
*  25   (uplo,             n,s1,v1,   m1)                   _SYR
*  26   (uplo,             n,u1,v1,   m1)                   _HER
*  27   (uplo,             n,s1,v1,v2,m1)                   _SYR2, _HER2
*
*  31   (          transa,transb,     m,n,k,s1,m1,m2,s2,m3) _GEMM
*  32   (side,uplo,                   m,n,  s1,m1,m2,s2,m3) _SYMM, _HEMM
*  33   (     uplo,trans,               n,k,s1,m1,   s2,m3) _SYRK
*  34   (     uplo,trans,               n,k,u1,m1,   u2,m3) _HERK
*  35   (     uplo,trans,               n,k,s1,m1,m2,s2,m3) _SYR2K
*  36   (     uplo,trans,               n,k,s1,m1,m2,u2,m3) _HER2K
*  37   (                             m,n,  s1,m1,   s2,m3) _TRAN_
*  38   (side,uplo,transa,       diag,m,n,  s1,m1,m2      ) _TRMM, _TRSM
*  39   (          trans,             m,n,  s1,m1,   s2,m3) _GEADD
*  40   (     uplo,trans,             m,n,  s1,m1,   s2,m3) _TRADD
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
*     ..
*     .. Common Blocks ..
      CHARACTER*1        DIAG, SIDE, TRANSA, TRANSB, UPLO
      INTEGER            IA, IB, IC, INCX, INCY, ISCLR, IX, IY, JA, JB,
     $                   JC, JX, JY, KDIM, MDIM, NDIM
      DOUBLE PRECISION   USCLR, SCLR
      INTEGER            DESCA( DLEN_ ), DESCB( DLEN_ ), DESCC( DLEN_ ),
     $                   DESCX( DLEN_ ), DESCY( DLEN_ )
      DOUBLE PRECISION   A( 2, 2 ), B( 2, 2 ), C( 2, 2 ), X( 2 ), Y( 2 )
      COMMON             /pblasc/diag, side, transa, transb, uplo
      COMMON             /pblasd/desca, descb, descc, descx, descy
      COMMON             /pblasi/ia, ib, ic, incx, incy, isclr, ix, iy,
     $                   ja, jb, jc, jx, jy
      COMMON             /pblasm/a, b, c
      COMMON             /pblasn/kdim, mdim, ndim
      COMMON             /pblass/sclr, usclr
      COMMON             /pblasv/x, y
*     ..
*     .. Executable Statements ..
*
*     Level 1 PBLAS
*
      IF( scode.EQ.11 ) THEN
*
         CALL subptr( ndim, x, ix, jx, descx, incx, y, iy, jy, descy,
     $                incy )
*
      ELSE IF( scode.EQ.12 ) THEN
*
         CALL subptr( ndim, sclr, x, ix, jx, descx, incx )
*
      ELSE IF( scode.EQ.13 ) THEN
*
         CALL subptr( ndim, sclr, x, ix, jx, descx, incx, y, iy, jy,
     $                descy, incy )
*
      ELSE IF( scode.EQ.14 ) THEN
*
         CALL subptr( ndim, sclr, isclr, x, ix, jx, descx, incx )
*
      ELSE IF( scode.EQ.15 ) THEN
*
         CALL subptr( ndim, usclr, x, ix, jx, descx, incx )
*
*     Level 2 PBLAS
*
      ELSE IF( scode.EQ.21 ) THEN
*
         CALL subptr( transa, mdim, ndim, sclr, a, ia, ja, desca, x, ix,
     $                jx, descx, incx, sclr, y, iy, jy, descy, incy )
*
      ELSE IF( scode.EQ.22 ) THEN
*
         CALL subptr( uplo, ndim, sclr, a, ia, ja, desca, x, ix, jx,
     $                descx, incx, sclr, y, iy, jy, descy, incy )
*
      ELSE IF( scode.EQ.23 ) THEN
*
         CALL subptr( uplo, transa, diag, ndim, a, ia, ja, desca, x, ix,
     $                jx, descx, incx )
*
      ELSE IF( scode.EQ.24 ) THEN
*
         CALL subptr( mdim, ndim, sclr, x, ix, jx, descx, incx, y, iy,
     $                jy, descy, incy, a, ia, ja, desca )
*
      ELSE IF( scode.EQ.25 ) THEN
*
         CALL subptr( uplo, ndim, sclr, x, ix, jx, descx, incx, a, ia,
     $                ja, desca )
*
      ELSE IF( scode.EQ.26 ) THEN
*
         CALL subptr( uplo, ndim, usclr, x, ix, jx, descx, incx, a, ia,
     $                ja, desca )
*
      ELSE IF( scode.EQ.27 ) THEN
*
         CALL subptr( uplo, ndim, sclr, x, ix, jx, descx, incx, y, iy,
     $                jy, descy, incy, a, ia, ja, desca )
*
*     Level 3 PBLAS
*
      ELSE IF( scode.EQ.31 ) THEN
*
         CALL subptr( transa, transb, mdim, ndim, kdim, sclr, a, ia, ja,
     $                desca, b, ib, jb, descb, sclr, c, ic, jc, descc )
*
      ELSE IF( scode.EQ.32 ) THEN
*
         CALL subptr( side, uplo, mdim, ndim, sclr, a, ia, ja, desca, b,
     $                ib, jb, descb, sclr, c, ic, jc, descc )
*
      ELSE IF( scode.EQ.33 ) THEN
*
         CALL subptr( uplo, transa, ndim, kdim, sclr, a, ia, ja, desca,
     $                sclr, c, ic, jc, descc )
*
      ELSE IF( scode.EQ.34 ) THEN
*
         CALL subptr( uplo, transa, ndim, kdim, usclr, a, ia, ja, desca,
     $                usclr, c, ic, jc, descc )
*
      ELSE IF( scode.EQ.35 ) THEN
*
         CALL subptr( uplo, transa, ndim, kdim, sclr, a, ia, ja, desca,
     $                b, ib, jb, descb, sclr, c, ic, jc, descc )
*
      ELSE IF( scode.EQ.36 ) THEN
*
         CALL subptr( uplo, transa, ndim, kdim, sclr, a, ia, ja, desca,
     $                b, ib, jb, descb, usclr, c, ic, jc, descc )
*
      ELSE IF( scode.EQ.37 ) THEN
*
         CALL subptr( mdim, ndim, sclr, a, ia, ja, desca, sclr, c, ic,
     $                jc, descc )
*
      ELSE IF( scode.EQ.38 ) THEN
*
         CALL subptr( side, uplo, transa, diag, mdim, ndim, sclr, a, ia,
     $                ja, desca, b, ib, jb, descb )
*
      ELSE IF( scode.EQ.39 ) THEN
*
         CALL subptr( transa, mdim, ndim, sclr, a, ia, ja, desca, sclr,
     $                c, ic, jc, descc )
*
      ELSE IF( scode.EQ.40 ) THEN
*
         CALL subptr( uplo, transa, mdim, ndim, sclr, a, ia, ja, desca,
     $                sclr, c, ic, jc, descc )
*
      END IF
*
      RETURN
*
*     End of PDCALLSUB
*
      END
      SUBROUTINE pderrset( ERR, ERRMAX, XTRUE, X )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   ERR, ERRMAX, X, XTRUE
*     ..
*
*  Purpose
*  =======
*
*  PDERRSET  computes the absolute difference ERR = |XTRUE - X| and com-
*  pares it with zero. ERRMAX accumulates the absolute error difference.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ERR     (local output) DOUBLE PRECISION
*          On exit, ERR specifies the absolute difference |XTRUE - X|.
*
*  ERRMAX  (local input/local output) DOUBLE PRECISION
*          On entry,  ERRMAX  specifies  a previously computed error. On
*          exit ERRMAX is the accumulated error MAX( ERRMAX, ERR ).
*
*  XTRUE   (local input) DOUBLE PRECISION
*          On entry, XTRUE specifies the true value.
*
*  X       (local input) DOUBLE PRECISION
*          On entry, X specifies the value to be compared to XTRUE.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. External Functions ..
      DOUBLE PRECISION   PDDIFF
      EXTERNAL           PDDIFF
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
      err = abs( pddiff( xtrue, x ) )
*
      errmax = max( errmax, err )
*
      RETURN
*
*     End of PDERRSET
*
      END
      SUBROUTINE pdchkvin( ERRMAX, N, X, PX, IX, JX, DESCX, INCX,
     $                     INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            INCX, INFO, IX, JX, N
      DOUBLE PRECISION   ERRMAX
*     ..
*     .. Array Arguments ..
      INTEGER            DESCX( * )
      DOUBLE PRECISION   PX( * ), X( * )
*     ..
*
*  Purpose
*  =======
*
*  PDCHKVIN  checks that the submatrix sub( PX ) remained unchanged. The
*  local  array  entries are compared element by element, and their dif-
*  ference  is tested against 0.0 as well as the epsilon machine. Notice
*  that  this difference should be numerically exactly the zero machine,
*  but  because of the possible fluctuation of some of the data we flag-
*  ged differently a difference less than twice the epsilon machine. The
*  largest error is also returned.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ERRMAX  (global output) DOUBLE PRECISION
*          On exit,  ERRMAX  specifies the largest absolute element-wise
*          difference between sub( X ) and sub( PX ).
*
*  N       (global input) INTEGER
*          On entry,  N  specifies  the  length of the subvector operand
*          sub( X ). N must be at least zero.
*
*  X       (local input) DOUBLE PRECISION array
*          On entry, X is an array of  dimension  (DESCX( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PX.
*
*  PX      (local input) DOUBLE PRECISION array
*          On entry, PX is an array of dimension (DESCX( LLD_ ),*). This
*          array contains the local entries of the matrix PX.
*
*  IX      (global input) INTEGER
*          On entry, IX  specifies X's global row index, which points to
*          the beginning of the submatrix sub( X ).
*
*  JX      (global input) INTEGER
*          On entry, JX  specifies X's global column index, which points
*          to the beginning of the submatrix sub( X ).
*
*  DESCX   (global and local input) INTEGER array
*          On entry, DESCX  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix X.
*
*  INCX    (global input) INTEGER
*          On entry,  INCX   specifies  the  global  increment  for  the
*          elements of  X.  Only two values of  INCX   are  supported in
*          this version, namely 1 and M_X. INCX  must not be zero.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO = 0, no error has been found,
*          If INFO > 0, the maximum abolute error found is in (0,eps],
*          If INFO < 0, the maximum abolute error found is in (eps,+oo).
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, ROWREP
      INTEGER            I, IB, ICTXT, ICURCOL, ICURROW, IIX, IN, IXCOL,
     $                   IXROW, J, JB, JJX, JN, KK, LDPX, LDX, LL,
     $                   MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION   ERR, EPS
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, pb_infog2l, pderrset
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           pdlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod
*     ..
*     .. Executable Statements ..
*
      info = 0
      errmax = zero
*
*     Quick return if possible
*
      IF( n.LE.0 )
     $   RETURN
*
      ictxt = descx( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      CALL pb_infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix,
     $                 jjx, ixrow, ixcol )
*
      ldx    = descx( m_ )
      ldpx   = descx( lld_ )
      rowrep = ( ixrow.EQ.-1 )
      colrep = ( ixcol.EQ.-1 )
*
      IF( n.EQ.1 ) THEN
*
         IF( ( myrow.EQ.ixrow .OR. rowrep ) .AND.
     $       ( mycol.EQ.ixcol .OR. colrep ) )
     $      CALL pderrset( err, errmax, x( ix+(jx-1)*ldx ),
     $                     px( iix+(jjx-1)*ldpx ) )
*
      ELSE IF( incx.EQ.descx( m_ ) ) THEN
*
*        sub( X ) is a row vector
*
         jb = descx( inb_ ) - jx + 1
         IF( jb.LE.0 )
     $      jb = ( ( -jb ) / descx( nb_ ) + 1 ) * descx( nb_ ) + jb
         jb = min( jb, n )
         jn = jx + jb - 1
*
         IF( myrow.EQ.ixrow .OR. rowrep ) THEN
*
            icurcol = ixcol
            IF( mycol.EQ.icurcol .OR. colrep ) THEN
               DO 10 j = jx, jn
                  CALL pderrset( err, errmax, x( ix+(j-1)*ldx ),
     $                           px( iix+(jjx-1)*ldpx ) )
                  jjx = jjx + 1
   10          CONTINUE
            END IF
            icurcol = mod( icurcol+1, npcol )
*
            DO 30 j = jn+1, jx+n-1, descx( nb_ )
               jb = min( jx+n-j, descx( nb_ ) )
*
               IF( mycol.EQ.icurcol .OR. colrep ) THEN
*
                  DO 20 kk = 0, jb-1
                     CALL pderrset( err, errmax, x( ix+(j+kk-1)*ldx ),
     $                              px( iix+(jjx+kk-1)*ldpx ) )
   20             CONTINUE
*
                  jjx = jjx + jb
*
               END IF
*
               icurcol = mod( icurcol+1, npcol )
*
   30       CONTINUE
*
         END IF
*
      ELSE
*
*        sub( X ) is a column vector
*
         ib = descx( imb_ ) - ix + 1
         IF( ib.LE.0 )
     $      ib = ( ( -ib ) / descx( mb_ ) + 1 ) * descx( mb_ ) + ib
         ib = min( ib, n )
         in = ix + ib - 1
*
         IF( mycol.EQ.ixcol .OR. colrep ) THEN
*
            icurrow = ixrow
            IF( myrow.EQ.icurrow .OR. rowrep ) THEN
               DO 40 i = ix, in
                  CALL pderrset( err, errmax, x( i+(jx-1)*ldx ),
     $                           px( iix+(jjx-1)*ldpx ) )
                  iix = iix + 1
   40          CONTINUE
            END IF
            icurrow = mod( icurrow+1, nprow )
*
            DO 60 i = in+1, ix+n-1, descx( mb_ )
               ib = min( ix+n-i, descx( mb_ ) )
*
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
*
                  DO 50 kk = 0, ib-1
                     CALL pderrset( err, errmax, x( i+kk+(jx-1)*ldx ),
     $                              px( iix+kk+(jjx-1)*ldpx ) )
   50             CONTINUE
*
                  iix = iix + ib
*
               END IF
*
               icurrow = mod( icurrow+1, nprow )
*
   60       CONTINUE
*
         END IF
*
      END IF
*
      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, errmax, 1, kk, ll, -1,
     $              -1, -1 )
*
      IF( errmax.GT.zero .AND. errmax.LE.eps ) THEN
         info = 1
      ELSE IF( errmax.GT.eps ) THEN
         info = -1
      END IF
*
      RETURN
*
*     End of PDCHKVIN
*
      END
      SUBROUTINE pdchkvout( N, X, PX, IX, JX, DESCX, INCX, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            INCX, INFO, IX, JX, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCX( * )
      DOUBLE PRECISION   PX( * ), X( * )
*     ..
*
*  Purpose
*  =======
*
*  PDCHKVOUT  checks  that the matrix PX \ sub( PX ) remained unchanged.
*  The  local array  entries  are compared element by element, and their
*  difference  is tested against 0.0 as well as the epsilon machine. No-
*  tice that this  difference should be numerically exactly the zero ma-
*  chine, but because  of  the  possible movement of some of the data we
*  flagged differently a difference less than twice the epsilon machine.
*  The largest error is reported.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          On entry,  N  specifies  the  length of the subvector operand
*          sub( X ). N must be at least zero.
*
*  X       (local input) DOUBLE PRECISION array
*          On entry, X is an array of  dimension  (DESCX( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PX.
*
*  PX      (local input) DOUBLE PRECISION array
*          On entry, PX is an array of dimension (DESCX( LLD_ ),*). This
*          array contains the local entries of the matrix PX.
*
*  IX      (global input) INTEGER
*          On entry, IX  specifies X's global row index, which points to
*          the beginning of the submatrix sub( X ).
*
*  JX      (global input) INTEGER
*          On entry, JX  specifies X's global column index, which points
*          to the beginning of the submatrix sub( X ).
*
*  DESCX   (global and local input) INTEGER array
*          On entry, DESCX  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix X.
*
*  INCX    (global input) INTEGER
*          On entry,  INCX   specifies  the  global  increment  for  the
*          elements of  X.  Only two values of  INCX   are  supported in
*          this version, namely 1 and M_X. INCX  must not be zero.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO = 0, no error has been found,
*          If INFO > 0, the maximum abolute error found is in (0,eps],
*          If INFO < 0, the maximum abolute error found is in (eps,+oo).
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, ROWREP
      INTEGER            I, IB, ICTXT, ICURCOL, ICURROW, II, IMBX, INBX,
     $                   J, JB, JJ, KK, LDPX, LDX, LL, MBX, MPALL,
     $                   MYCOL, MYCOLDIST, MYROW, MYROWDIST, NBX, NPCOL,
     $                   nprow, nqall
      DOUBLE PRECISION   EPS, ERR, ERRMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, DGAMX2D, PDERRSET
*     ..
*     .. External Functions ..
      INTEGER            PB_NUMROC
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           PDLAMCH, PB_NUMROC
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod
*     ..
*     .. Executable Statements ..
*
      info = 0
      errmax = zero
*
*     Quick return if possible
*
      IF( ( descx( m_ ).LE.0 ).OR.( descx( n_ ).LE.0 ) )
     $   RETURN
*
*     Start the operations
*
      ictxt = descx( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      mpall   = pb_numroc( descx( m_ ), 1, descx( imb_ ), descx( mb_ ),
     $                     myrow, descx( rsrc_ ), nprow )
      nqall   = pb_numroc( descx( n_ ), 1, descx( inb_ ), descx( nb_ ),
     $                     mycol, descx( csrc_ ), npcol )
*
      mbx     = descx( mb_ )
      nbx     = descx( nb_ )
      ldx     = descx( m_ )
      ldpx    = descx( lld_ )
      icurrow = descx( rsrc_ )
      icurcol = descx( csrc_ )
      rowrep  = ( icurrow.EQ.-1 )
      colrep  = ( icurcol.EQ.-1 )
      IF( myrow.EQ.icurrow .OR. rowrep ) THEN
         imbx = descx( imb_ )
      ELSE
         imbx = mbx
      END IF
      IF( mycol.EQ.icurcol .OR. colrep ) THEN
         inbx = descx( inb_ )
      ELSE
         inbx = nbx
      END IF
      IF( rowrep ) THEN
         myrowdist = 0
      ELSE
         myrowdist = mod( myrow - icurrow + nprow, nprow )
      END IF
      IF( colrep ) THEN
         mycoldist = 0
      ELSE
         mycoldist = mod( mycol - icurcol + npcol, npcol )
      END IF
      ii = 1
      jj = 1
*
      IF( incx.EQ.descx( m_ ) ) THEN
*
*        sub( X ) is a row vector
*
         IF( myrow.EQ.icurrow .OR. rowrep ) THEN
*
            i = 1
            IF( mycoldist.EQ.0 ) THEN
               j = 1
            ELSE
               j = descx( inb_ ) + ( mycoldist - 1 ) * nbx + 1
            END IF
            jb = min( max( 0, descx( n_ ) - j + 1 ), inbx )
            ib = min( descx( m_ ), descx( imb_ ) )
*
            DO 20 kk = 0, jb-1
               DO 10 ll = 0, ib-1
                  IF( i+ll.NE.ix .OR. j+kk.LT.jx .OR. j+kk.GT.jx+n-1 )
     $               CALL pderrset( err, errmax,
     $                              x( i+ll+(j+kk-1)*ldx ),
     $                              px( ii+ll+(jj+kk-1)*ldpx ) )
   10          CONTINUE
   20       CONTINUE
            IF( colrep ) THEN
               j = j + inbx
            ELSE
               j = j + inbx + ( npcol - 1 ) * nbx
            END IF
*
            DO 50 jj = inbx+1, nqall, nbx
               jb = min( nqall-jj+1, nbx )
*
               DO 40 kk = 0, jb-1
                  DO 30 ll = 0, ib-1
                     IF( i+ll.NE.ix .OR. j+kk.LT.jx .OR.
     $                   j+kk.GT.jx+n-1 )
     $                  CALL pderrset( err, errmax,
     $                                 x( i+ll+(j+kk-1)*ldx ),
     $                                 px( ii+ll+(jj+kk-1)*ldpx ) )
   30             CONTINUE
   40          CONTINUE
*
               IF( colrep ) THEN
                  j = j + nbx
               ELSE
                  j = j + npcol * nbx
               END IF
*
   50       CONTINUE
*
            ii = ii + ib
*
         END IF
*
         icurrow = mod( icurrow + 1, nprow )
*
         DO 110 i = descx( imb_ ) + 1, descx( m_ ), mbx
            ib = min( descx( m_ ) - i + 1, mbx )
*
            IF( myrow.EQ.icurrow .OR. rowrep ) THEN
*
               IF( mycoldist.EQ.0 ) THEN
                  j = 1
               ELSE
                  j = descx( inb_ ) + ( mycoldist - 1 ) * nbx + 1
               END IF
*
               jj = 1
               jb = min( max( 0, descx( n_ ) - j + 1 ), inbx )
               DO 70 kk = 0, jb-1
                  DO 60 ll = 0, ib-1
                     IF( i+ll.NE.ix .OR. j+kk.LT.jx .OR.
     $                   j+kk.GT.jx+n-1 )
     $                  CALL pderrset( err, errmax,
     $                                 x( i+ll+(j+kk-1)*ldx ),
     $                                 px( ii+ll+(jj+kk-1)*ldpx ) )
   60             CONTINUE
   70          CONTINUE
               IF( colrep ) THEN
                  j = j + inbx
               ELSE
                  j = j + inbx + ( npcol - 1 ) * nbx
               END IF
*
               DO 100 jj = inbx+1, nqall, nbx
                  jb = min( nqall-jj+1, nbx )
*
                  DO 90 kk = 0, jb-1
                     DO 80 ll = 0, ib-1
                        IF( i+ll.NE.ix .OR. j+kk.LT.jx .OR.
     $                      j+kk.GT.jx+n-1 )
     $                     CALL pderrset( err, errmax,
     $                                    x( i+ll+(j+kk-1)*ldx ),
     $                                    px( ii+ll+(jj+kk-1)*ldpx ) )
   80                CONTINUE
   90             CONTINUE
*
                  IF( colrep ) THEN
                     j = j + nbx
                  ELSE
                     j = j + npcol * nbx
                  END IF
*
  100          CONTINUE
*
               ii = ii + ib
*
            END IF
*
            icurrow = mod( icurrow + 1, nprow )
*
  110    CONTINUE
*
      ELSE
*
*        sub( X ) is a column vector
*
         IF( mycol.EQ.icurcol .OR. colrep ) THEN
*
            j = 1
            IF( myrowdist.EQ.0 ) THEN
               i = 1
            ELSE
               i = descx( imb_ ) + ( myrowdist - 1 ) * mbx + 1
            END IF
            ib = min( max( 0, descx( m_ ) - i + 1 ), imbx )
            jb = min( descx( n_ ), descx( inb_ ) )
*
            DO 130 kk = 0, jb-1
               DO 120 ll = 0, ib-1
                  IF( j+kk.NE.jx .OR. i+ll.LT.ix .OR. i+ll.GT.ix+n-1 )
     $               CALL pderrset( err, errmax,
     $                              x( i+ll+(j+kk-1)*ldx ),
     $                              px( ii+ll+(jj+kk-1)*ldpx ) )
  120          CONTINUE
  130       CONTINUE
            IF( rowrep ) THEN
               i = i + imbx
            ELSE
               i = i + imbx + ( nprow - 1 ) * mbx
            END IF
*
            DO 160 ii = imbx+1, mpall, mbx
               ib = min( mpall-ii+1, mbx )
*
               DO 150 kk = 0, jb-1
                  DO 140 ll = 0, ib-1
                     IF( j+kk.NE.jx .OR. i+ll.LT.ix .OR.
     $                   i+ll.GT.ix+n-1 )
     $                  CALL pderrset( err, errmax,
     $                                 x( i+ll+(j+kk-1)*ldx ),
     $                                 px( ii+ll+(jj+kk-1)*ldpx ) )
  140             CONTINUE
  150          CONTINUE
*
               IF( rowrep ) THEN
                  i = i + mbx
               ELSE
                  i = i + nprow * mbx
               END IF
*
  160       CONTINUE
*
            jj = jj + jb
*
         END IF
*
         icurcol = mod( icurcol + 1, npcol )
*
         DO 220 j = descx( inb_ ) + 1, descx( n_ ), nbx
            jb = min( descx( n_ ) - j + 1, nbx )
*
            IF( mycol.EQ.icurcol .OR. colrep ) THEN
*
               IF( myrowdist.EQ.0 ) THEN
                  i = 1
               ELSE
                  i = descx( imb_ ) + ( myrowdist - 1 ) * mbx + 1
               END IF
*
               ii = 1
               ib = min( max( 0, descx( m_ ) - i + 1 ), imbx )
               DO 180 kk = 0, jb-1
                  DO 170 ll = 0, ib-1
                     IF( j+kk.NE.jx .OR. i+ll.LT.ix .OR.
     $                   i+ll.GT.ix+n-1 )
     $                  CALL pderrset( err, errmax,
     $                                 x( i+ll+(j+kk-1)*ldx ),
     $                                 px( ii+ll+(jj+kk-1)*ldpx ) )
  170             CONTINUE
  180          CONTINUE
               IF( rowrep ) THEN
                  i = i + imbx
               ELSE
                  i = i + imbx + ( nprow - 1 ) * mbx
               END IF
*
               DO 210 ii = imbx+1, mpall, mbx
                  ib = min( mpall-ii+1, mbx )
*
                  DO 200 kk = 0, jb-1
                     DO 190 ll = 0, ib-1
                        IF( j+kk.NE.jx .OR. i+ll.LT.ix .OR.
     $                      i+ll.GT.ix+n-1 )
     $                     CALL pderrset( err, errmax,
     $                                    x( i+ll+(j+kk-1)*ldx ),
     $                                    px( ii+ll+(jj+kk-1)*ldpx ) )
  190                CONTINUE
  200             CONTINUE
*
                  IF( rowrep ) THEN
                     i = i + mbx
                  ELSE
                     i = i + nprow * mbx
                  END IF
*
  210          CONTINUE
*
               jj = jj + jb
*
            END IF
*
            icurcol = mod( icurcol + 1, npcol )
*
  220    CONTINUE
*
      END IF
*
      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, errmax, 1, kk, ll, -1,
     $              -1, -1 )
*
      IF( errmax.GT.zero .AND. errmax.LE.eps ) THEN
         info = 1
      ELSE IF( errmax.GT.eps ) THEN
         info = -1
      END IF
*
      RETURN
*
*     End of PDCHKVOUT
*
      END
      SUBROUTINE pdchkmin( ERRMAX, M, N, A, PA, IA, JA, DESCA, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            IA, INFO, JA, M, N
      DOUBLE PRECISION   ERRMAX
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      DOUBLE PRECISION   PA( * ), A( * )
*     ..
*
*  Purpose
*  =======
*
*  PDCHKMIN  checks that the submatrix sub( PA ) remained unchanged. The
*  local  array  entries are compared element by element, and their dif-
*  ference  is tested against 0.0 as well as the epsilon machine. Notice
*  that  this difference should be numerically exactly the zero machine,
*  but  because of the possible fluctuation of some of the data we flag-
*  ged differently a difference less than twice the epsilon machine. The
*  largest error is also returned.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ERRMAX  (global output) DOUBLE PRECISION
*          On exit,  ERRMAX  specifies the largest absolute element-wise
*          difference between sub( A ) and sub( PA ).
*
*  M       (global input) INTEGER
*          On entry,  M  specifies  the  number of rows of the submatrix
*          operand sub( A ). M must be at least zero.
*
*  N       (global input) INTEGER
*          On entry, N  specifies the number of columns of the submatrix
*          operand sub( A ). N must be at least zero.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  PA      (local input) DOUBLE PRECISION array
*          On entry, PA is an array of dimension (DESCA( LLD_ ),*). This
*          array contains the local entries of the matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO = 0, no error has been found,
*          If INFO > 0, the maximum abolute error found is in (0,eps],
*          If INFO < 0, the maximum abolute error found is in (eps,+oo).
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, ROWREP
      INTEGER            H, I, IACOL, IAROW, IB, ICTXT, ICURCOL,
     $                   ICURROW, II, IIA, IN, J, JB, JJ, JJA, JN, K,
     $                   KK, LDA, LDPA, LL, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION   ERR, EPS
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, pb_infog2l, pderrset
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           pdlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod
*     ..
*     .. Executable Statements ..
*
      info   = 0
      errmax = zero
*
*     Quick return if posssible
*
      IF( ( m.EQ.0 ).OR.( n.EQ.0 ) )
     $   RETURN
*
*     Start the operations
*
      ictxt = desca( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      CALL pb_infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia,
     $                 jja, iarow, iacol )
*
      ii      = iia
      jj      = jja
      lda     = desca( m_ )
      ldpa    = desca( lld_ )
      icurrow = iarow
      icurcol = iacol
      rowrep  = ( iarow.EQ.-1 )
      colrep  = ( iacol.EQ.-1 )
*
*     Handle the first block of column separately
*
      jb = desca( inb_ ) - ja  + 1
      IF( jb.LE.0 )
     $   jb = ( ( -jb ) / desca( nb_ ) + 1 ) * desca( nb_ ) + jb
      jb = min( jb, n )
      jn = ja + jb - 1
*
      IF( mycol.EQ.icurcol .OR. colrep ) THEN
*
         DO 40 h = 0, jb-1
            ib = desca( imb_ ) - ia  + 1
            IF( ib.LE.0 )
     $         ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
            ib = min( ib, m )
            in = ia + ib - 1
            IF( myrow.EQ.icurrow .OR. rowrep ) THEN
               DO 10 k = 0, ib-1
                  CALL pderrset( err, errmax, a( ia+k+(ja+h-1)*lda ),
     $                           pa( ii+k+(jj+h-1)*ldpa ) )
   10          CONTINUE
               ii = ii + ib
            END IF
            icurrow = mod( icurrow+1, nprow )
*
*           Loop over remaining block of rows
*
            DO 30 i = in+1, ia+m-1, desca( mb_ )
               ib = min( desca( mb_ ), ia+m-i )
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  DO 20 k = 0, ib-1
                     CALL pderrset( err, errmax, a( i+k+(ja+h-1)*lda ),
     $                              pa( ii+k+(jj+h-1)*ldpa ) )
   20             CONTINUE
                  ii = ii + ib
               END IF
               icurrow = mod( icurrow+1, nprow )
   30       CONTINUE
*
            ii = iia
            icurrow = iarow
   40    CONTINUE
*
         jj = jj + jb
*
      END IF
*
      icurcol = mod( icurcol+1, npcol )
*
*     Loop over remaining column blocks
*
      DO 90 j = jn+1, ja+n-1, desca( nb_ )
         jb = min(  desca( nb_ ), ja+n-j )
         IF( mycol.EQ.icurcol .OR. colrep ) THEN
            DO 80 h = 0, jb-1
               ib = desca( imb_ ) - ia  + 1
               IF( ib.LE.0 )
     $            ib = ( ( -ib ) / desca( mb_ ) + 1 )*desca( mb_ ) + ib
               ib = min( ib, m )
               in = ia + ib - 1
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  DO 50 k = 0, ib-1
                     CALL pderrset( err, errmax, a( ia+k+(j+h-1)*lda ),
     $                              pa( ii+k+(jj+h-1)*ldpa ) )
   50             CONTINUE
                  ii = ii + ib
               END IF
               icurrow = mod( icurrow+1, nprow )
*
*              Loop over remaining block of rows
*
               DO 70 i = in+1, ia+m-1, desca( mb_ )
                  ib = min( desca( mb_ ), ia+m-i )
                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                     DO 60 k = 0, ib-1
                        CALL pderrset( err, errmax,
     $                                 a( i+k+(j+h-1)*lda ),
     $                                 pa( ii+k+(jj+h-1)*ldpa ) )
   60                CONTINUE
                     ii = ii + ib
                  END IF
                  icurrow = mod( icurrow+1, nprow )
   70          CONTINUE
*
               ii = iia
               icurrow = iarow
   80       CONTINUE
*
            jj = jj + jb
         END IF
*
         icurcol = mod( icurcol+1, npcol )
*
   90 CONTINUE
*
      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, errmax, 1, kk, ll, -1,
     $              -1, -1 )
*
      IF( errmax.GT.zero .AND. errmax.LE.eps ) THEN
         info = 1
      ELSE IF( errmax.GT.eps ) THEN
         info = -1
      END IF
*
      RETURN
*
*     End of PDCHKMIN
*
      END
      SUBROUTINE pdchkmout( M, N, A, PA, IA, JA, DESCA, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            IA, INFO, JA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      DOUBLE PRECISION   A( * ), PA( * )
*     ..
*
*  Purpose
*  =======
*
*  PDCHKMOUT  checks  that the matrix PA \ sub( PA ) remained unchanged.
*  The  local array  entries  are compared element by element, and their
*  difference  is tested against 0.0 as well as the epsilon machine. No-
*  tice that this  difference should be numerically exactly the zero ma-
*  chine, but because  of  the  possible movement of some of the data we
*  flagged differently a difference less than twice the epsilon machine.
*  The largest error is reported.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  M       (global input) INTEGER
*          On entry,  M  specifies  the  number of rows of the submatrix
*          sub( PA ). M must be at least zero.
*
*  N       (global input) INTEGER
*          On entry, N specifies the  number of columns of the submatrix
*          sub( PA ). N must be at least zero.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  PA      (local input) DOUBLE PRECISION array
*          On entry, PA is an array of dimension (DESCA( LLD_ ),*). This
*          array contains the local entries of the matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO = 0, no error has been found,
*          If INFO > 0, the maximum abolute error found is in (0,eps],
*          If INFO < 0, the maximum abolute error found is in (eps,+oo).
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, ROWREP
      INTEGER            I, IB, ICTXT, ICURCOL, II, IMBA, J, JB, JJ, KK,
     $                   LDA, LDPA, LL, MPALL, MYCOL, MYROW, MYROWDIST,
     $                   NPCOL, NPROW
      DOUBLE PRECISION   EPS, ERR, ERRMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, pderrset
*     ..
*     .. External Functions ..
      INTEGER            PB_NUMROC
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           PDLAMCH, PB_NUMROC
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min, mod
*     ..
*     .. Executable Statements ..
*
      info = 0
      errmax = zero
*
*     Quick return if possible
*
      IF( ( desca( m_ ).LE.0 ).OR.( desca( n_ ).LE.0 ) )
     $   RETURN
*
*     Start the operations
*
      ictxt = desca( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      mpall = pb_numroc( desca( m_ ), 1, desca( imb_ ), desca( mb_ ),
     $                   myrow, desca( rsrc_ ), nprow )
*
      lda    = desca( m_ )
      ldpa   = desca( lld_ )
*
      ii = 1
      jj = 1
      rowrep  = ( desca( rsrc_ ).EQ.-1 )
      colrep  = ( desca( csrc_ ).EQ.-1 )
      icurcol = desca( csrc_ )
      IF( myrow.EQ.desca( rsrc_ ) .OR. rowrep ) THEN
         imba = desca( imb_ )
      ELSE
         imba = desca( mb_ )
      END IF
      IF( rowrep ) THEN
         myrowdist = 0
      ELSE
         myrowdist = mod( myrow - desca( rsrc_ ) + nprow, nprow )
      END IF
*
      IF( mycol.EQ.icurcol .OR. colrep ) THEN
*
         j = 1
         IF( myrowdist.EQ.0 ) THEN
            i = 1
         ELSE
            i = desca( imb_ ) + ( myrowdist - 1 ) * desca( mb_ ) + 1
         END IF
         ib = min( max( 0, desca( m_ ) - i + 1 ), imba )
         jb = min( desca( n_ ), desca( inb_ ) )
*
         DO 20 kk = 0, jb-1
            DO 10 ll = 0, ib-1
               IF( i+ll.LT.ia .OR. i+ll.GT.ia+m-1 .OR.
     $             j+kk.LT.ja .OR. j+kk.GT.ja+n-1 )
     $            CALL pderrset( err, errmax, a( i+ll+(j+kk-1)*lda ),
     $                           pa( ii+ll+(jj+kk-1)*ldpa ) )
   10       CONTINUE
   20    CONTINUE
         IF( rowrep ) THEN
            i = i + imba
         ELSE
            i = i + imba + ( nprow - 1 ) * desca( mb_ )
         END IF
*
         DO 50 ii = imba + 1, mpall, desca( mb_ )
            ib = min( mpall-ii+1, desca( mb_ ) )
*
            DO 40 kk = 0, jb-1
               DO 30 ll = 0, ib-1
                  IF( i+ll.LT.ia .OR. i+ll.GT.ia+m-1 .OR.
     $                j+kk.LT.ja .OR. j+kk.GT.ja+n-1 )
     $               CALL pderrset( err, errmax,
     $                              a( i+ll+(j+kk-1)*lda ),
     $                              pa( ii+ll+(jj+kk-1)*ldpa ) )
   30          CONTINUE
   40       CONTINUE
*
            IF( rowrep ) THEN
               i = i + desca( mb_ )
            ELSE
               i = i + nprow * desca( mb_ )
            END IF
*
   50    CONTINUE
*
         jj = jj + jb
*
      END IF
*
      icurcol = mod( icurcol + 1, npcol )
*
      DO 110 j = desca( inb_ ) + 1, desca( n_ ), desca( nb_ )
         jb = min( desca( n_ ) - j + 1, desca( nb_ ) )
*
         IF( mycol.EQ.icurcol .OR. colrep ) THEN
*
            IF( myrowdist.EQ.0 ) THEN
               i = 1
            ELSE
               i = desca( imb_ ) + ( myrowdist - 1 ) * desca( mb_ ) + 1
            END IF
*
            ii = 1
            ib = min( max( 0, desca( m_ ) - i + 1 ), imba )
            DO 70 kk = 0, jb-1
               DO 60 ll = 0, ib-1
                  IF( i+ll.LT.ia .OR. i+ll.GT.ia+m-1 .OR.
     $                j+kk.LT.ja .OR. j+kk.GT.ja+n-1 )
     $               CALL pderrset( err, errmax,
     $                              a( i+ll+(j+kk-1)*lda ),
     $                              pa( ii+ll+(jj+kk-1)*ldpa ) )
   60          CONTINUE
   70       CONTINUE
            IF( rowrep ) THEN
               i = i + imba
            ELSE
               i = i + imba + ( nprow - 1 ) * desca( mb_ )
            END IF
*
            DO 100 ii = imba+1, mpall, desca( mb_ )
               ib = min( mpall-ii+1, desca( mb_ ) )
*
               DO 90 kk = 0, jb-1
                  DO 80 ll = 0, ib-1
                     IF( i+ll.LT.ia .OR. i+ll.GT.ia+m-1 .OR.
     $                   j+kk.LT.ja .OR. j+kk.GT.ja+n-1 )
     $                  CALL pderrset( err, errmax,
     $                                 a( i+ll+(j+kk-1)*lda ),
     $                                 pa( ii+ll+(jj+kk-1)*ldpa ) )
   80             CONTINUE
   90          CONTINUE
*
               IF( rowrep ) THEN
                  i = i + desca( mb_ )
               ELSE
                  i = i + nprow * desca( mb_ )
               END IF
*
  100       CONTINUE
*
            jj = jj + jb
*
         END IF
*
         icurcol = mod( icurcol + 1, npcol )
*                                                           INSERT MODE
  110 CONTINUE
*
      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, errmax, 1, kk, ll, -1,
     $              -1, -1 )
*
      IF( errmax.GT.zero .AND. errmax.LE.eps ) THEN
         info = 1
      ELSE IF( errmax.GT.eps ) THEN
         info = -1
      END IF
*
      RETURN
*
*     End of PDCHKMOUT
*
      END
      SUBROUTINE pdmprnt( ICTXT, NOUT, M, N, A, LDA, IRPRNT, ICPRNT,
     $                    CMATNM )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            ICPRNT, ICTXT, IRPRNT, LDA, M, N, NOUT
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)      CMATNM
      DOUBLE PRECISION   A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  PDMPRNT prints to the standard output an array A of size m by n. Only
*  the process of coordinates ( IRPRNT, ICPRNT ) is printing.
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  M       (global input) INTEGER
*          On entry, M  specifies the number of rows of the matrix A.  M
*          must be at least zero.
*
*  N       (global input) INTEGER
*          On entry, N  specifies the number of columns of the matrix A.
*          N must be at least zero.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry,  A  is an array of dimension (LDA,N). The leading m
*          by n part of this array is printed.
*
*  LDA     (local input) INTEGER
*          On entry, LDA  specifies the leading dimension of  the  local
*          array A to be printed. LDA must be at least MAX( 1, M ).
*
*  IRPRNT  (global input) INTEGER
*          On entry, IRPRNT  specifies the process row coordinate of the
*          printing process.
*
*  ICPRNT  (global input) INTEGER
*          On entry,  ICPRNT  specifies the process column coordinate of
*          the printing process.
*
*  CMATNM  (global input) CHARACTER*(*)
*          On entry, CMATNM specifies the identifier of the matrix to be
*          printed.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, J, MYCOL, MYROW, NPCOL, NPROW
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( ( m.LE.0 ).OR.( n.LE.0 ) )
     $   RETURN
*
*     Get grid parameters
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      IF( myrow.EQ.irprnt .AND. mycol.EQ.icprnt ) THEN
*
         WRITE( nout, fmt = * )
         DO 20 j = 1, n
*
            DO 10 i = 1, m
*
               WRITE( nout, fmt = 9999 ) cmatnm, i, j, a( i, j )
*
   10       CONTINUE
*
   20    CONTINUE
*
      END IF
*
 9999 FORMAT( 1x, a, '(', i6, ',', i6, ')=', d30.18 )
*
      RETURN
*
*     End of PDMPRNT
*
      END
      SUBROUTINE pdvprnt( ICTXT, NOUT, N, X, INCX, IRPRNT, ICPRNT,
     $                    CVECNM )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      INTEGER            ICPRNT, ICTXT, INCX, IRPRNT, N, NOUT
*     ..
*     .. Array Arguments ..
      CHARACTER*(*)      CVECNM
      DOUBLE PRECISION   X( * )
*     ..
*
*  Purpose
*  =======
*
*  PDVPRNT  prints  to the standard output an vector x of length n. Only
*  the process of coordinates ( IRPRNT, ICPRNT ) is printing.
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  NOUT    (global input) INTEGER
*          On entry, NOUT specifies the unit number for the output file.
*          When NOUT is 6, output to screen,  when  NOUT is 0, output to
*          stderr. NOUT is only defined for process 0.
*
*  N       (global input) INTEGER
*          On entry, N  specifies the length of the vector X.  N must be
*          at least zero.
*
*  X       (global input) DOUBLE PRECISION array
*          On   entry,   X   is   an   array   of   dimension  at  least
*          ( 1 + ( n - 1 )*abs( INCX ) ).  Before  entry,  the incremen-
*          ted array X must contain the vector x.
*
*  INCX    (global input) INTEGER.
*          On entry, INCX specifies the increment for the elements of X.
*          INCX must not be zero.
*
*  IRPRNT  (global input) INTEGER
*          On entry, IRPRNT  specifies the process row coordinate of the
*          printing process.
*
*  ICPRNT  (global input) INTEGER
*          On entry,  ICPRNT  specifies the process column coordinate of
*          the printing process.
*
*  CVECNM  (global input) CHARACTER*(*)
*          On entry, CVECNM specifies the identifier of the vector to be
*          printed.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, MYCOL, MYROW, NPCOL, NPROW
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( n.LE.0 )
     $   RETURN
*
*     Get grid parameters
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      IF( myrow.EQ.irprnt .AND. mycol.EQ.icprnt ) THEN
*
         WRITE( nout, fmt = * )
         DO 10 i = 1, 1 + ( n-1 )*incx, incx
*
            WRITE( nout, fmt = 9999 ) cvecnm, i, x( i )
*
   10    CONTINUE
*
      END IF
*
 9999 FORMAT( 1x, a, '(', i6, ')=', d30.18 )
*
      RETURN
*
*     End of PDVPRNT
*
      END
      SUBROUTINE pdmvch( ICTXT, TRANS, M, N, ALPHA, A, IA, JA, DESCA,
     $                   X, IX, JX, DESCX, INCX, BETA, Y, PY, IY, JY,
     $                   DESCY, INCY, G, ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        TRANS
      INTEGER            IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
     $                   JY, M, N
      DOUBLE PRECISION   ALPHA, BETA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCX( * ), DESCY( * )
      DOUBLE PRECISION   A( * ), G( * ), PY( * ), X( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  PDMVCH checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  TRANS   (global input) CHARACTER*1
*          On entry,  TRANS  specifies which matrix-vector product is to
*          be computed as follows:
*             If TRANS = 'N',
*                sub( Y ) = BETA * sub( Y ) + sub( A )  * sub( X ),
*             otherwise
*                sub( Y ) = BETA * sub( Y ) + sub( A )' * sub( X ).
*
*  M       (global input) INTEGER
*          On entry,  M  specifies  the  number of rows of the submatrix
*          operand matrix A. M must be at least zero.
*
*  N       (global input) INTEGER
*          On entry,  N  specifies  the  number of columns of the subma-
*          trix operand matrix A. N must be at least zero.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  X       (local input) DOUBLE PRECISION array
*          On entry, X is an array of  dimension  (DESCX( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PX.
*
*  IX      (global input) INTEGER
*          On entry, IX  specifies X's global row index, which points to
*          the beginning of the submatrix sub( X ).
*
*  JX      (global input) INTEGER
*          On entry, JX  specifies X's global column index, which points
*          to the beginning of the submatrix sub( X ).
*
*  DESCX   (global and local input) INTEGER array
*          On entry, DESCX  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix X.
*
*  INCX    (global input) INTEGER
*          On entry,  INCX   specifies  the  global  increment  for  the
*          elements of  X.  Only two values of  INCX   are  supported in
*          this version, namely 1 and M_X. INCX  must not be zero.
*
*  BETA    (global input) DOUBLE PRECISION
*          On entry, BETA specifies the scalar beta.
*
*  Y       (local input/local output) DOUBLE PRECISION array
*          On entry, Y is an array of  dimension  (DESCY( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PY.
*
*  PY      (local input) DOUBLE PRECISION array
*          On entry, PY is an array of dimension (DESCY( LLD_ ),*). This
*          array contains the local entries of the matrix PY.
*
*  IY      (global input) INTEGER
*          On entry, IY  specifies Y's global row index, which points to
*          the beginning of the submatrix sub( Y ).
*
*  JY      (global input) INTEGER
*          On entry, JY  specifies Y's global column index, which points
*          to the beginning of the submatrix sub( Y ).
*
*  DESCY   (global and local input) INTEGER array
*          On entry, DESCY  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix Y.
*
*  INCY    (global input) INTEGER
*          On entry,  INCY   specifies  the  global  increment  for  the
*          elements of  Y.  Only two values of  INCY   are  supported in
*          this version, namely 1 and M_Y. INCY  must not be zero.
*
*  G       (workspace) DOUBLE PRECISION array
*          On entry, G is an array of dimension at least MAX( M, N ).  G
*          is used to compute the gauges.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, ROWREP, TRAN
      INTEGER            I, IB, ICURCOL, ICURROW, IIY, IN, IOFFA, IOFFX,
     $                   IOFFY, IYCOL, IYROW, J, JB, JJY, JN, KK, LDA,
     $                   LDPY, LDX, LDY, ML, MYCOL, MYROW, NL, NPCOL,
     $                   nprow
      DOUBLE PRECISION   EPS, ERRI, GTMP, TBETA, YTMP
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, pb_infog2l
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           lsame, pdlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      IF( m.EQ.0 .OR. n.EQ.0 ) THEN
         tbeta = one
      ELSE
         tbeta = beta
      END IF
*
      tran = lsame( trans, 'T' ).OR.lsame( trans, 'C' )
      IF( tran ) THEN
         ml = n
         nl = m
      ELSE
         ml = m
         nl = n
      END IF
*
      lda = max( 1, desca( m_ ) )
      ldx = max( 1, descx( m_ ) )
      ldy = max( 1, descy( m_ ) )
*
*     Compute expected result in Y using data in A, X and Y.
*     Compute gauges in G. This part of the computation is performed
*     by every process in the grid.
*
      ioffy = iy + ( jy - 1 ) * ldy
      DO 30 i = 1, ml
         ytmp = zero
         gtmp = zero
         ioffx = ix + ( jx - 1 ) * ldx
         IF( tran )THEN
            ioffa = ia + ( ja + i - 2 ) * lda
            DO 10 j = 1, nl
               ytmp = ytmp + a( ioffa ) * x( ioffx )
               gtmp = gtmp + abs( a( ioffa ) * x( ioffx ) )
               ioffa = ioffa + 1
               ioffx = ioffx + incx
   10       CONTINUE
         ELSE
            ioffa = ia + i - 1 + ( ja - 1 ) * lda
            DO 20 j = 1, nl
               ytmp = ytmp + a( ioffa ) * x( ioffx )
               gtmp = gtmp + abs( a( ioffa ) * x( ioffx ) )
               ioffa = ioffa + lda
               ioffx = ioffx + incx
   20       CONTINUE
         END IF
         g( i ) = abs( alpha ) * gtmp + abs( tbeta * y( ioffy ) )
         y( ioffy ) = alpha * ytmp + tbeta * y( ioffy )
         ioffy = ioffy + incy
   30 CONTINUE
*
*     Compute the error ratio for this result.
*
      err  = zero
      info = 0
      ldpy = descy( lld_ )
      ioffy = iy + ( jy - 1 ) * ldy
      CALL pb_infog2l( iy, jy, descy, nprow, npcol, myrow, mycol, iiy,
     $                 jjy, iyrow, iycol )
      icurrow = iyrow
      icurcol = iycol
      rowrep  = ( iyrow.EQ.-1 )
      colrep  = ( iycol.EQ.-1 )
*
      IF( incy.EQ.descy( m_ ) ) THEN
*
*        sub( Y ) is a row vector
*
         jb = descy( inb_ ) - jy + 1
         IF( jb.LE.0 )
     $      jb = ( ( -jb ) / descy( nb_ ) + 1 ) * descy( nb_ ) + jb
         jb = min( jb, ml )
         jn = jy + jb - 1
*
         DO 50 j = jy, jn
*
            IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
     $          ( mycol.EQ.icurcol .OR. colrep ) ) THEN
               erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) ) / eps
               IF( g( j-jy+1 ).NE.zero )
     $            erri = erri / g( j-jy+1 )
               err = max( err, erri )
               IF( err*sqrt( eps ).GE.one )
     $            info = 1
               jjy = jjy + 1
            END IF
*
            ioffy = ioffy + incy
*
   50    CONTINUE
*
         icurcol = mod( icurcol+1, npcol )
*
         DO 70 j = jn+1, jy+ml-1, descy( nb_ )
            jb = min( jy+ml-j, descy( nb_ ) )
*
            DO 60 kk = 0, jb-1
*
               IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
     $             ( mycol.EQ.icurcol .OR. colrep ) ) THEN
                  erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) )/eps
                  IF( g( j+kk-jy+1 ).NE.zero )
     $               erri = erri / g( j+kk-jy+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  jjy = jjy + 1
               END IF
*
               ioffy = ioffy + incy
*
   60       CONTINUE
*
            icurcol = mod( icurcol+1, npcol )
*
   70    CONTINUE
*
      ELSE
*
*        sub( Y ) is a column vector
*
         ib = descy( imb_ ) - iy + 1
         IF( ib.LE.0 )
     $      ib = ( ( -ib ) / descy( mb_ ) + 1 ) * descy( mb_ ) + ib
         ib = min( ib, ml )
         in = iy + ib - 1
*
         DO 80 i = iy, in
*
            IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
     $          ( mycol.EQ.icurcol .OR. colrep ) ) THEN
               erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) ) / eps
               IF( g( i-iy+1 ).NE.zero )
     $            erri = erri / g( i-iy+1 )
               err = max( err, erri )
               IF( err*sqrt( eps ).GE.one )
     $            info = 1
               iiy = iiy + 1
            END IF
*
            ioffy = ioffy + incy
*
   80    CONTINUE
*
         icurrow = mod( icurrow+1, nprow )
*
         DO 100 i = in+1, iy+ml-1, descy( mb_ )
            ib = min( iy+ml-i, descy( mb_ ) )
*
            DO 90 kk = 0, ib-1
*
               IF( ( myrow.EQ.icurrow .OR. rowrep ) .AND.
     $             ( mycol.EQ.icurcol .OR. colrep ) ) THEN
                  erri = abs( py( iiy+(jjy-1)*ldpy ) - y( ioffy ) )/eps
                  IF( g( i+kk-iy+1 ).NE.zero )
     $               erri = erri / g( i+kk-iy+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  iiy = iiy + 1
               END IF
*
               ioffy = ioffy + incy
*
   90       CONTINUE
*
            icurrow = mod( icurrow+1, nprow )
*
  100    CONTINUE
*
      END IF
*
*     If INFO = 0, all results are at least half accurate.
*
      CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $              mycol )
*
      RETURN
*
*     End of PDMVCH
*
      END
      SUBROUTINE pdvmch( ICTXT, UPLO, M, N, ALPHA, X, IX, JX, DESCX,
     $                   INCX, Y, IY, JY, DESCY, INCY, A, PA, IA, JA,
     $                   DESCA, G, ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        UPLO
      INTEGER            IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
     $                   JY, M, N
      DOUBLE PRECISION   ALPHA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCX( * ), DESCY( * )
      DOUBLE PRECISION   A( * ), G( * ), PA( * ), X( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  PDVMCH checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  UPLO    (global input) CHARACTER*1
*          On entry, UPLO specifies which part of the submatrix sub( A )
*          is to be referenced as follows:
*             If UPLO = 'L', only the lower triangular part,
*             If UPLO = 'U', only the upper triangular part,
*             else the entire matrix is to be referenced.
*
*  M       (global input) INTEGER
*          On entry,  M  specifies  the  number of rows of the submatrix
*          operand matrix A. M must be at least zero.
*
*  N       (global input) INTEGER
*          On entry,  N  specifies  the  number of columns of the subma-
*          trix operand matrix A. N must be at least zero.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  X       (local input) DOUBLE PRECISION array
*          On entry, X is an array of  dimension  (DESCX( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PX.
*
*  IX      (global input) INTEGER
*          On entry, IX  specifies X's global row index, which points to
*          the beginning of the submatrix sub( X ).
*
*  JX      (global input) INTEGER
*          On entry, JX  specifies X's global column index, which points
*          to the beginning of the submatrix sub( X ).
*
*  DESCX   (global and local input) INTEGER array
*          On entry, DESCX  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix X.
*
*  INCX    (global input) INTEGER
*          On entry,  INCX   specifies  the  global  increment  for  the
*          elements of  X.  Only two values of  INCX   are  supported in
*          this version, namely 1 and M_X. INCX  must not be zero.
*
*  Y       (local input) DOUBLE PRECISION array
*          On entry, Y is an array of  dimension  (DESCY( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PY.
*
*  IY      (global input) INTEGER
*          On entry, IY  specifies Y's global row index, which points to
*          the beginning of the submatrix sub( Y ).
*
*  JY      (global input) INTEGER
*          On entry, JY  specifies Y's global column index, which points
*          to the beginning of the submatrix sub( Y ).
*
*  DESCY   (global and local input) INTEGER array
*          On entry, DESCY  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix Y.
*
*  INCY    (global input) INTEGER
*          On entry,  INCY   specifies  the  global  increment  for  the
*          elements of  Y.  Only two values of  INCY   are  supported in
*          this version, namely 1 and M_Y. INCY  must not be zero.
*
*  A       (local input/local output) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  PA      (local input) DOUBLE PRECISION array
*          On entry, PA is an array of dimension (DESCA( LLD_ ),*). This
*          array contains the local entries of the matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  G       (workspace) DOUBLE PRECISION array
*          On entry, G is an array of dimension at least MAX( M, N ).  G
*          is used to compute the gauges.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, LOWER, ROWREP, UPPER
      INTEGER            I, IACOL, IAROW, IB, IBEG, ICURROW, IEND, IIA,
     $                   in, ioffa, ioffx, ioffy, j, jja, kk, lda, ldpa,
     $                   ldx, ldy, mycol, myrow, npcol, nprow
      DOUBLE PRECISION   ATMP, EPS, ERRI, GTMP
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, pb_infog2l
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           LSAME, PDLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      upper = lsame( uplo, 'U' )
      lower = lsame( uplo, 'L' )
*
      lda = max( 1, desca( m_ ) )
      ldx = max( 1, descx( m_ ) )
      ldy = max( 1, descy( m_ ) )
*
*     Compute expected result in A using data in A, X and Y.
*     Compute gauges in G. This part of the computation is performed
*     by every process in the grid.
*
      DO 70 j = 1, n
*
         ioffy = iy + ( jy - 1 ) * ldy + ( j - 1 ) * incy
*
         IF( lower ) THEN
            ibeg = j
            iend = m
            DO 10 i = 1, j-1
               g( i ) = zero
   10       CONTINUE
         ELSE IF( upper ) THEN
            ibeg = 1
            iend = j
            DO 20 i = j+1, m
               g( i ) = zero
   20       CONTINUE
         ELSE
            ibeg = 1
            iend = m
         END IF
*
         DO 30 i = ibeg, iend
*
            ioffx = ix + ( jx - 1 ) * ldx + ( i - 1 ) * incx
            ioffa = ia + i - 1 + ( ja + j - 2 ) * lda
            atmp = x( ioffx ) * y( ioffy )
            gtmp = abs( x( ioffx ) * y( ioffy ) )
            g( i ) = abs( alpha ) * gtmp + abs( a( ioffa ) )
            a( ioffa ) = alpha * atmp + a( ioffa )
*
   30    CONTINUE
*
*        Compute the error ratio for this result.
*
         info = 0
         err  = zero
         ldpa = desca( lld_ )
         ioffa = ia + ( ja + j - 2 ) * lda
         CALL pb_infog2l( ia, ja+j-1, desca, nprow, npcol, myrow, mycol,
     $                    iia, jja, iarow, iacol )
         rowrep = ( iarow.EQ.-1 )
         colrep = ( iacol.EQ.-1 )
*
         IF( mycol.EQ.iacol .OR. colrep ) THEN
*
            icurrow = iarow
            ib = desca( imb_ ) - ia + 1
            IF( ib.LE.0 )
     $         ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
            ib = min( ib, m )
            in = ia + ib - 1
*
            DO 40 i = ia, in
*
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  erri = abs( pa( iia+(jja-1)*ldpa ) - a( ioffa ) )/eps
                  IF( g( i-ia+1 ).NE.zero )
     $               erri = erri / g( i-ia+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  iia = iia + 1
               END IF
*
               ioffa = ioffa + 1
*
   40       CONTINUE
*
            icurrow = mod( icurrow+1, nprow )
*
            DO 60 i = in+1, ia+m-1, desca( mb_ )
               ib = min( ia+m-i, desca( mb_ ) )
*
               DO 50 kk = 0, ib-1
*
                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                     erri = abs( pa( iia+(jja-1)*ldpa )-a( ioffa ) )/eps
                     IF( g( i+kk-ia+1 ).NE.zero )
     $                  erri = erri / g( i+kk-ia+1 )
                     err = max( err, erri )
                     IF( err*sqrt( eps ).GE.one )
     $                  info = 1
                     iia = iia + 1
                  END IF
*
                  ioffa = ioffa + 1
*
   50          CONTINUE
*
               icurrow = mod( icurrow+1, nprow )
*
   60       CONTINUE
*
         END IF
*
*        If INFO = 0, all results are at least half accurate.
*
         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
         CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $                 mycol )
         IF( info.NE.0 )
     $      GO TO 80
*
   70 CONTINUE
*
   80 CONTINUE
*
      RETURN
*
*     End of PDVMCH
*
      END
      SUBROUTINE pdvmch2( ICTXT, UPLO, M, N, ALPHA, X, IX, JX, DESCX,
     $                    INCX, Y, IY, JY, DESCY, INCY, A, PA, IA,
     $                    JA, DESCA, G, ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        UPLO
      INTEGER            IA, ICTXT, INCX, INCY, INFO, IX, IY, JA, JX,
     $                   jy, m, n
      DOUBLE PRECISION   ALPHA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCX( * ), DESCY( * )
      DOUBLE PRECISION   A( * ), G( * ), PA( * ), X( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  PDVMCH2 checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  UPLO    (global input) CHARACTER*1
*          On entry, UPLO specifies which part of the submatrix sub( A )
*          is to be referenced as follows:
*             If UPLO = 'L', only the lower triangular part,
*             If UPLO = 'U', only the upper triangular part,
*             else the entire matrix is to be referenced.
*
*  M       (global input) INTEGER
*          On entry,  M  specifies  the  number of rows of the submatrix
*          operand matrix A. M must be at least zero.
*
*  N       (global input) INTEGER
*          On entry,  N  specifies  the  number of columns of the subma-
*          trix operand matrix A. N must be at least zero.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  X       (local input) DOUBLE PRECISION array
*          On entry, X is an array of  dimension  (DESCX( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PX.
*
*  IX      (global input) INTEGER
*          On entry, IX  specifies X's global row index, which points to
*          the beginning of the submatrix sub( X ).
*
*  JX      (global input) INTEGER
*          On entry, JX  specifies X's global column index, which points
*          to the beginning of the submatrix sub( X ).
*
*  DESCX   (global and local input) INTEGER array
*          On entry, DESCX  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix X.
*
*  INCX    (global input) INTEGER
*          On entry,  INCX   specifies  the  global  increment  for  the
*          elements of  X.  Only two values of  INCX   are  supported in
*          this version, namely 1 and M_X. INCX  must not be zero.
*
*  Y       (local input) DOUBLE PRECISION array
*          On entry, Y is an array of  dimension  (DESCY( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PY.
*
*  IY      (global input) INTEGER
*          On entry, IY  specifies Y's global row index, which points to
*          the beginning of the submatrix sub( Y ).
*
*  JY      (global input) INTEGER
*          On entry, JY  specifies Y's global column index, which points
*          to the beginning of the submatrix sub( Y ).
*
*  DESCY   (global and local input) INTEGER array
*          On entry, DESCY  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix Y.
*
*  INCY    (global input) INTEGER
*          On entry,  INCY   specifies  the  global  increment  for  the
*          elements of  Y.  Only two values of  INCY   are  supported in
*          this version, namely 1 and M_Y. INCY  must not be zero.
*
*  A       (local input/local output) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  PA      (local input) DOUBLE PRECISION array
*          On entry, PA is an array of dimension (DESCA( LLD_ ),*). This
*          array contains the local entries of the matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  G       (workspace) DOUBLE PRECISION array
*          On entry, G is an array of dimension at least MAX( M, N ).  G
*          is used to compute the gauges.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, LOWER, ROWREP, UPPER
      INTEGER            I, IACOL, IAROW, IB, IBEG, ICURROW, IEND, IIA,
     $                   IN, IOFFA, IOFFXI, IOFFXJ, IOFFYI, IOFFYJ, J,
     $                   JJA, KK, LDA, LDPA, LDX, LDY, MYCOL, MYROW,
     $                   npcol, nprow
      DOUBLE PRECISION   EPS, ERRI, GTMP, ATMP
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, pb_infog2l
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           lsame, pdlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      upper = lsame( uplo, 'U' )
      lower = lsame( uplo, 'L' )
*
      lda = max( 1, desca( m_ ) )
      ldx = max( 1, descx( m_ ) )
      ldy = max( 1, descy( m_ ) )
*
*     Compute expected result in A using data in A, X and Y.
*     Compute gauges in G. This part of the computation is performed
*     by every process in the grid.
*
      DO 70 j = 1, n
*
         ioffxj = ix + ( jx - 1 ) * ldx + ( j - 1 ) * incx
         ioffyj = iy + ( jy - 1 ) * ldy + ( j - 1 ) * incy
*
         IF( lower ) THEN
            ibeg = j
            iend = m
            DO 10 i = 1, j-1
               g( i ) = zero
   10       CONTINUE
         ELSE IF( upper ) THEN
            ibeg = 1
            iend = j
            DO 20 i = j+1, m
               g( i ) = zero
   20       CONTINUE
         ELSE
            ibeg = 1
            iend = m
         END IF
*
         DO 30 i = ibeg, iend
            ioffa = ia + i - 1 + ( ja + j - 2 ) * lda
            ioffxi = ix + ( jx - 1 ) * ldx + ( i - 1 ) * incx
            ioffyi = iy + ( jy - 1 ) * ldy + ( i - 1 ) * incy
            atmp = x( ioffxi ) * y( ioffyj )
            atmp = atmp + y( ioffyi ) * x( ioffxj )
            gtmp = abs( x( ioffxi ) * y( ioffyj ) )
            gtmp = gtmp + abs( y( ioffyi ) * x( ioffxj ) )
            g( i ) = abs( alpha ) * gtmp + abs( a( ioffa ) )
            a( ioffa ) = alpha*atmp + a( ioffa )
*
   30    CONTINUE
*
*        Compute the error ratio for this result.
*
         info = 0
         err  = zero
         ldpa = desca( lld_ )
         ioffa = ia + ( ja + j - 2 ) * lda
         CALL pb_infog2l( ia, ja+j-1, desca, nprow, npcol, myrow, mycol,
     $                    iia, jja, iarow, iacol )
         rowrep = ( iarow.EQ.-1 )
         colrep = ( iacol.EQ.-1 )
*
         IF( mycol.EQ.iacol .OR. colrep ) THEN
*
            icurrow = iarow
            ib = desca( imb_ ) - ia + 1
            IF( ib.LE.0 )
     $         ib = ( ( -ib ) / desca( mb_ ) + 1 ) * desca( mb_ ) + ib
            ib = min( ib, m )
            in = ia + ib - 1
*
            DO 40 i = ia, in
*
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  erri = abs( pa( iia+(jja-1)*ldpa ) - a( ioffa ) )/eps
                  IF( g( i-ia+1 ).NE.zero )
     $               erri = erri / g( i-ia+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  iia = iia + 1
               END IF
*
               ioffa = ioffa + 1
*
   40       CONTINUE
*
            icurrow = mod( icurrow+1, nprow )
*
            DO 60 i = in+1, ia+m-1, desca( mb_ )
               ib = min( ia+m-i, desca( mb_ ) )
*
               DO 50 kk = 0, ib-1
*
                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                     erri = abs( pa( iia+(jja-1)*ldpa )-a( ioffa ) )/eps
                     IF( g( i+kk-ia+1 ).NE.zero )
     $                  erri = erri / g( i+kk-ia+1 )
                     err = max( err, erri )
                     IF( err*sqrt( eps ).GE.one )
     $                  info = 1
                     iia = iia + 1
                  END IF
*
                  ioffa = ioffa + 1
*
   50          CONTINUE
*
               icurrow = mod( icurrow+1, nprow )
*
   60       CONTINUE
*
         END IF
*
*        If INFO = 0, all results are at least half accurate.
*
         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
         CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $                 mycol )
         IF( info.NE.0 )
     $      GO TO 80
*
   70 CONTINUE
*
   80 CONTINUE
*
      RETURN
*
*     End of PDVMCH2
*
      END
      SUBROUTINE pdmmch( ICTXT, TRANSA, TRANSB, M, N, K, ALPHA, A, IA,
     $                   JA, DESCA, B, IB, JB, DESCB, BETA, C, PC, IC,
     $                   JC, DESCC, CT, G, ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        TRANSA, TRANSB
      INTEGER            IA, IB, IC, ICTXT, INFO, JA, JB, JC, K, M, N
      DOUBLE PRECISION   ALPHA, BETA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCB( * ), DESCC( * )
      DOUBLE PRECISION   A( * ), B( * ), C( * ), CT( * ), G( * ),
     $                   PC( * )
*     ..
*
*  Purpose
*  =======
*
*  PDMMCH checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  TRANSA  (global input) CHARACTER*1
*          On entry, TRANSA specifies if the matrix  operand  A is to be
*          transposed.
*
*  TRANSB  (global input) CHARACTER*1
*          On entry, TRANSB specifies if the matrix  operand  B is to be
*          transposed.
*
*  M       (global input) INTEGER
*          On entry, M specifies the number of rows of C.
*
*  N       (global input) INTEGER
*          On entry, N specifies the number of columns of C.
*
*  K       (global input) INTEGER
*          On entry, K specifies the number of columns (resp. rows) of A
*          when  TRANSA = 'N'  (resp. TRANSA <> 'N')  in PxGEMM, PxSYRK,
*          PxSYR2K, PxHERK and PxHER2K.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  B       (local input) DOUBLE PRECISION array
*          On entry, B is an array of  dimension  (DESCB( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PB.
*
*  IB      (global input) INTEGER
*          On entry, IB  specifies B's global row index, which points to
*          the beginning of the submatrix sub( B ).
*
*  JB      (global input) INTEGER
*          On entry, JB  specifies B's global column index, which points
*          to the beginning of the submatrix sub( B ).
*
*  DESCB   (global and local input) INTEGER array
*          On entry, DESCB  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix B.
*
*  BETA    (global input) DOUBLE PRECISION
*          On entry, BETA specifies the scalar beta.
*
*  C       (local input/local output) DOUBLE PRECISION array
*          On entry, C is an array of  dimension  (DESCC( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PC.
*
*  PC      (local input) DOUBLE PRECISION array
*          On entry, PC is an array of dimension (DESCC( LLD_ ),*). This
*          array contains the local pieces of the matrix PC.
*
*  IC      (global input) INTEGER
*          On entry, IC  specifies C's global row index, which points to
*          the beginning of the submatrix sub( C ).
*
*  JC      (global input) INTEGER
*          On entry, JC  specifies C's global column index, which points
*          to the beginning of the submatrix sub( C ).
*
*  DESCC   (global and local input) INTEGER array
*          On entry, DESCC  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix C.
*
*  CT      (workspace) DOUBLE PRECISION array
*          On entry, CT is an array of dimension at least MAX(M,N,K). CT
*          holds a copy of the current column of C.
*
*  G       (workspace) DOUBLE PRECISION array
*          On entry, G  is  an array of dimension at least MAX(M,N,K). G
*          is used to compute the gauges.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, ROWREP, TRANA, TRANB
      INTEGER            I, IBB, ICCOL, ICROW, ICURROW, IIC, IN, IOFFA,
     $                   IOFFB, IOFFC, J, JJC, KK, LDA, LDB, LDC, LDPC,
     $                   mycol, myrow, npcol, nprow
      DOUBLE PRECISION   EPS, ERRI
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, pb_infog2l
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           LSAME, PDLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      trana = lsame( transa, 'T' ).OR.lsame( transa, 'C' )
      tranb = lsame( transb, 'T' ).OR.lsame( transb, 'C' )
*
      lda = max( 1, desca( m_ ) )
      ldb = max( 1, descb( m_ ) )
      ldc = max( 1, descc( m_ ) )
*
*     Compute expected result in C using data in A, B and C.
*     Compute gauges in G. This part of the computation is performed
*     by every process in the grid.
*
      DO 240 j = 1, n
*
         ioffc = ic + ( jc + j - 2 ) * ldc
         DO 10 i = 1, m
            ct( i ) = zero
            g( i )  = zero
   10    CONTINUE
*
         IF( .NOT.trana .AND. .NOT.tranb ) THEN
            DO 30 kk = 1, k
               ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
               DO 20 i = 1, m
                  ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
                  ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
                  g( i ) = g( i ) + abs( a( ioffa ) ) *
     $                     abs( b( ioffb ) )
   20          CONTINUE
   30       CONTINUE
         ELSE IF( trana .AND. .NOT.tranb ) THEN
            DO 50 kk = 1, k
               ioffb = ib + kk - 1 + ( jb + j - 2 ) * ldb
               DO 40 i = 1, m
                  ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
                  ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
                  g( i ) = g( i ) + abs( a( ioffa ) ) *
     $                     abs( b( ioffb ) )
   40          CONTINUE
   50       CONTINUE
         ELSE IF( .NOT.trana .AND. tranb ) THEN
            DO 70 kk = 1, k
               ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
               DO 60 i = 1, m
                  ioffa = ia + i - 1 + ( ja + kk - 2 ) * lda
                  ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
                  g( i ) = g( i ) + abs( a( ioffa ) ) *
     $                     abs( b( ioffb ) )
   60          CONTINUE
   70       CONTINUE
         ELSE IF( trana .AND. tranb ) THEN
            DO 90 kk = 1, k
               ioffb = ib + j - 1 + ( jb + kk - 2 ) * ldb
               DO 80 i = 1, m
                  ioffa = ia + kk - 1 + ( ja + i - 2 ) * lda
                  ct( i ) = ct( i ) + a( ioffa ) * b( ioffb )
                  g( i ) = g( i ) + abs( a( ioffa ) ) *
     $                     abs( b( ioffb ) )
   80          CONTINUE
   90       CONTINUE
         END IF
*
         DO 200 i = 1, m
            ct( i ) = alpha*ct( i ) + beta * c( ioffc )
            g( i ) = abs( alpha )*g( i ) + abs( beta )*abs( c( ioffc ) )
            c( ioffc ) = ct( i )
            ioffc      = ioffc + 1
  200    CONTINUE
*
*        Compute the error ratio for this result.
*
         err  = zero
         info = 0
         ldpc = descc( lld_ )
         ioffc = ic + ( jc + j - 2 ) * ldc
         CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
     $                    iic, jjc, icrow, iccol )
         icurrow = icrow
         rowrep  = ( icrow.EQ.-1 )
         colrep  = ( iccol.EQ.-1 )
*
         IF( mycol.EQ.iccol .OR. colrep ) THEN
*
            ibb = descc( imb_ ) - ic + 1
            IF( ibb.LE.0 )
     $         ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
            ibb = min( ibb, m )
            in = ic + ibb - 1
*
            DO 210 i = ic, in
*
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  erri = abs( pc( iic+(jjc-1)*ldpc ) -
     $                        c( ioffc ) ) / eps
                  IF( g( i-ic+1 ).NE.zero )
     $               erri = erri / g( i-ic+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  iic = iic + 1
               END IF
*
               ioffc = ioffc + 1
*
  210       CONTINUE
*
            icurrow = mod( icurrow+1, nprow )
*
            DO 230 i = in+1, ic+m-1, descc( mb_ )
               ibb = min( ic+m-i, descc( mb_ ) )
*
               DO 220 kk = 0, ibb-1
*
                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                     erri = abs( pc( iic+(jjc-1)*ldpc ) -
     $                           c( ioffc ) )/eps
                     IF( g( i+kk-ic+1 ).NE.zero )
     $                  erri = erri / g( i+kk-ic+1 )
                     err = max( err, erri )
                     IF( err*sqrt( eps ).GE.one )
     $                  info = 1
                     iic = iic + 1
                  END IF
*
                  ioffc = ioffc + 1
*
  220          CONTINUE
*
               icurrow = mod( icurrow+1, nprow )
*
  230       CONTINUE
*
         END IF
*
*        If INFO = 0, all results are at least half accurate.
*
         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
         CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $                 mycol )
         IF( info.NE.0 )
     $      GO TO 250
*
  240 CONTINUE
*
  250 CONTINUE
*
      RETURN
*
*     End of PDMMCH
*
      END
      SUBROUTINE pdmmch1( ICTXT, UPLO, TRANS, N, K, ALPHA, A, IA, JA,
     $                    DESCA, BETA, C, PC, IC, JC, DESCC, CT, G,
     $                    ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        TRANS, UPLO
      INTEGER            IA, IC, ICTXT, INFO, JA, JC, K, N
      DOUBLE PRECISION   ALPHA, BETA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCC( * )
      DOUBLE PRECISION   A( * ), C( * ), CT( * ), G( * ), PC( * )
*     ..
*
*  Purpose
*  =======
*
*  PDMMCH1 checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  UPLO    (global input) CHARACTER*1
*          On entry,  UPLO  specifies which part of C should contain the
*          result.
*
*  TRANS   (global input) CHARACTER*1
*          On entry,  TRANS  specifies  whether  the  matrix A has to be
*          transposed or not before computing the matrix-matrix product.
*
*  N       (global input) INTEGER
*          On entry, N  specifies  the order  the submatrix operand C. N
*          must be at least zero.
*
*  K       (global input) INTEGER
*          On entry, K specifies the number of columns (resp. rows) of A
*          when  TRANS = 'N'  (resp. TRANS <> 'N').  K  must be at least
*          zero.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  BETA    (global input) DOUBLE PRECISION
*          On entry, BETA specifies the scalar beta.
*
*  C       (local input/local output) DOUBLE PRECISION array
*          On entry, C is an array of  dimension  (DESCC( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PC.
*
*  PC      (local input) DOUBLE PRECISION array
*          On entry, PC is an array of dimension (DESCC( LLD_ ),*). This
*          array contains the local pieces of the matrix PC.
*
*  IC      (global input) INTEGER
*          On entry, IC  specifies C's global row index, which points to
*          the beginning of the submatrix sub( C ).
*
*  JC      (global input) INTEGER
*          On entry, JC  specifies C's global column index, which points
*          to the beginning of the submatrix sub( C ).
*
*  DESCC   (global and local input) INTEGER array
*          On entry, DESCC  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix C.
*
*  CT      (workspace) DOUBLE PRECISION array
*          On entry, CT is an array of dimension at least MAX(M,N,K). CT
*          holds a copy of the current column of C.
*
*  G       (workspace) DOUBLE PRECISION array
*          On entry, G  is  an array of dimension at least MAX(M,N,K). G
*          is used to compute the gauges.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, NOTRAN, ROWREP, TRAN, UPPER
      INTEGER            I, IBB, IBEG, ICCOL, ICROW, ICURROW, IEND, IIC,
     $                   IN, IOFFAK, IOFFAN, IOFFC, J, JJC, KK, LDA,
     $                   LDC, LDPC, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION   EPS, ERRI
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, pb_infog2l
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           lsame, pdlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      upper  = lsame( uplo,  'U' )
      notran = lsame( trans, 'N' )
      tran   = lsame( trans, 'T' )
*
      lda = max( 1, desca( m_ ) )
      ldc = max( 1, descc( m_ ) )
*
*     Compute expected result in C using data in A, B and C.
*     Compute gauges in G. This part of the computation is performed
*     by every process in the grid.
*
      DO 140 j = 1, n
*
         IF( upper ) THEN
            ibeg = 1
            iend = j
         ELSE
            ibeg = j
            iend = n
         END IF
*
         DO 10 i = 1, n
            ct( i ) = zero
            g( i )  = zero
   10    CONTINUE
*
         IF( notran ) THEN
            DO 30 kk = 1, k
               ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
               DO 20 i = ibeg, iend
                  ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
                  ct( i ) = ct( i ) + a( ioffak ) * a( ioffan )
                  g( i ) = g( i ) + abs( a( ioffak ) ) *
     $                     abs( a( ioffan ) )
   20          CONTINUE
   30       CONTINUE
         ELSE IF( tran ) THEN
            DO 50 kk = 1, k
               ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
               DO 40 i = ibeg, iend
                  ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
                  ct( i ) = ct( i ) + a( ioffak ) * a( ioffan )
                  g( i ) = g( i ) + abs( a( ioffak ) ) *
     $                     abs( a( ioffan ) )
   40          CONTINUE
   50       CONTINUE
         END IF
*
         ioffc = ic + ibeg - 1 + ( jc + j - 2 ) * ldc
*
         DO 100 i = ibeg, iend
            ct( i ) = alpha*ct( i ) + beta * c( ioffc )
            g( i ) = abs( alpha )*g( i ) + abs( beta )*abs( c( ioffc ) )
            c( ioffc ) = ct( i )
            ioffc = ioffc + 1
  100    CONTINUE
*
*        Compute the error ratio for this result.
*
         err  = zero
         info = 0
         ldpc = descc( lld_ )
         ioffc = ic + ( jc + j - 2 ) * ldc
         CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
     $                    iic, jjc, icrow, iccol )
         icurrow = icrow
         rowrep  = ( icrow.EQ.-1 )
         colrep  = ( iccol.EQ.-1 )
*
         IF( mycol.EQ.iccol .OR. colrep ) THEN
*
            ibb = descc( imb_ ) - ic + 1
            IF( ibb.LE.0 )
     $         ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
            ibb = min( ibb, n )
            in = ic + ibb - 1
*
            DO 110 i = ic, in
*
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  erri = abs( pc( iic+(jjc-1)*ldpc ) -
     $                        c( ioffc ) ) / eps
                  IF( g( i-ic+1 ).NE.zero )
     $               erri = erri / g( i-ic+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  iic = iic + 1
               END IF
*
               ioffc = ioffc + 1
*
  110       CONTINUE
*
            icurrow = mod( icurrow+1, nprow )
*
            DO 130 i = in+1, ic+n-1, descc( mb_ )
               ibb = min( ic+n-i, descc( mb_ ) )
*
               DO 120 kk = 0, ibb-1
*
                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                     erri = abs( pc( iic+(jjc-1)*ldpc ) -
     $                           c( ioffc ) )/eps
                     IF( g( i+kk-ic+1 ).NE.zero )
     $                  erri = erri / g( i+kk-ic+1 )
                     err = max( err, erri )
                     IF( err*sqrt( eps ).GE.one )
     $                  info = 1
                     iic = iic + 1
                  END IF
*
                  ioffc = ioffc + 1
*
  120          CONTINUE
*
               icurrow = mod( icurrow+1, nprow )
*
  130       CONTINUE
*
         END IF
*
*        If INFO = 0, all results are at least half accurate.
*
         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
         CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $                 mycol )
         IF( info.NE.0 )
     $      GO TO 150
*
  140 CONTINUE
*
  150 CONTINUE
*
      RETURN
*
*     End of PDMMCH1
*
      END
      SUBROUTINE pdmmch2( ICTXT, UPLO, TRANS, N, K, ALPHA, A, IA, JA,
     $                    DESCA, B, IB, JB, DESCB, BETA, C, PC, IC,
     $                    JC, DESCC, CT, G, ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        TRANS, UPLO
      INTEGER            IA, IB, IC, ICTXT, INFO, JA, JB, JC, K, N
      DOUBLE PRECISION   ALPHA, BETA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCB( * ), DESCC( * )
      DOUBLE PRECISION   A( * ), B( * ), C( * ), CT( * ), G( * ),
     $                   pc( * )
*     ..
*
*  Purpose
*  =======
*
*  PDMMCH2 checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  ICTXT   (local input) INTEGER
*          On entry,  ICTXT  specifies the BLACS context handle, indica-
*          ting the global  context of the operation. The context itself
*          is global, but the value of ICTXT is local.
*
*  UPLO    (global input) CHARACTER*1
*          On entry,  UPLO  specifies which part of C should contain the
*          result.
*
*  TRANS   (global input) CHARACTER*1
*          On entry,  TRANS  specifies whether the matrices A and B have
*          to  be  transposed  or not before computing the matrix-matrix
*          product.
*
*  N       (global input) INTEGER
*          On entry, N  specifies  the order  the submatrix operand C. N
*          must be at least zero.
*
*  K       (global input) INTEGER
*          On entry, K specifies the number of columns (resp. rows) of A
*          and B when  TRANS = 'N' (resp. TRANS <> 'N').  K  must  be at
*          least zero.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  B       (local input) DOUBLE PRECISION array
*          On entry, B is an array of  dimension  (DESCB( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PB.
*
*  IB      (global input) INTEGER
*          On entry, IB  specifies B's global row index, which points to
*          the beginning of the submatrix sub( B ).
*
*  JB      (global input) INTEGER
*          On entry, JB  specifies B's global column index, which points
*          to the beginning of the submatrix sub( B ).
*
*  DESCB   (global and local input) INTEGER array
*          On entry, DESCB  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix B.
*
*  BETA    (global input) DOUBLE PRECISION
*          On entry, BETA specifies the scalar beta.
*
*  C       (local input/local output) DOUBLE PRECISION array
*          On entry, C is an array of  dimension  (DESCC( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PC.
*
*  PC      (local input) DOUBLE PRECISION array
*          On entry, PC is an array of dimension (DESCC( LLD_ ),*). This
*          array contains the local pieces of the matrix PC.
*
*  IC      (global input) INTEGER
*          On entry, IC  specifies C's global row index, which points to
*          the beginning of the submatrix sub( C ).
*
*  JC      (global input) INTEGER
*          On entry, JC  specifies C's global column index, which points
*          to the beginning of the submatrix sub( C ).
*
*  DESCC   (global and local input) INTEGER array
*          On entry, DESCC  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix C.
*
*  CT      (workspace) DOUBLE PRECISION array
*          On entry, CT is an array of dimension at least MAX(M,N,K). CT
*          holds a copy of the current column of C.
*
*  G       (workspace) DOUBLE PRECISION array
*          On entry, G  is  an array of dimension at least MAX(M,N,K). G
*          is used to compute the gauges.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, NOTRAN, ROWREP, TRAN, UPPER
      INTEGER            I, IBB, IBEG, ICCOL, ICROW, ICURROW, IEND, IIC,
     $                   IN, IOFFAK, IOFFAN, IOFFBK, IOFFBN, IOFFC, J,
     $                   JJC, KK, LDA, LDB, LDC, LDPC, MYCOL, MYROW,
     $                   NPCOL, NPROW
      DOUBLE PRECISION   EPS, ERRI
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgamx2d, igsum2d, pb_infog2l
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           LSAME, PDLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      eps = pdlamch( ictxt, 'eps' )
*
      upper = lsame( uplo, 'U' )
      notran = lsame( trans, 'N' )
      tran = lsame( trans, 'T' )
*
      lda = max( 1, desca( m_ ) )
      ldb = max( 1, descb( m_ ) )
      ldc = max( 1, descc( m_ ) )
*
*     Compute expected result in C using data in A, B and C.
*     Compute gauges in G. This part of the computation is performed
*     by every process in the grid.
*
      DO 140 j = 1, n
*
         IF( upper ) THEN
            ibeg = 1
            iend = j
         ELSE
            ibeg = j
            iend = n
         END IF
*
         DO 10 i = 1, n
            ct( i ) = zero
            g( i )  = zero
   10    CONTINUE
*
         IF( notran ) THEN
            DO 30 kk = 1, k
               ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
               ioffbk = ib + j - 1 + ( jb + kk - 2 ) * ldb
               DO 20 i = ibeg, iend
                  ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
                  ioffbn = ib + i - 1 + ( jb + kk - 2 ) * ldb
                  ct( i ) = ct( i ) + alpha * (
     $                      a( ioffan ) * b( ioffbk ) +
     $                      b( ioffbn ) * a( ioffak ) )
                  g( i ) = g( i ) + abs( alpha ) * (
     $                     abs( a( ioffan ) ) * abs( b( ioffbk ) ) +
     $                     abs( b( ioffbn ) ) * abs( a( ioffak ) ) )
   20          CONTINUE
   30       CONTINUE
         ELSE IF( tran ) THEN
            DO 50 kk = 1, k
               ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
               ioffbk = ib + kk - 1 + ( jb + j - 2 ) * ldb
               DO 40 i = ibeg, iend
                  ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
                  ioffbn = ib + kk - 1 + ( jb + i - 2 ) * ldb
                  ct( i ) = ct( i ) + alpha * (
     $                      a( ioffan ) * b( ioffbk ) +
     $                      b( ioffbn ) * a( ioffak ) )
                  g( i ) = g( i ) + abs( alpha ) * (
     $                     abs( a( ioffan ) ) * abs( b( ioffbk ) ) +
     $                     abs( b( ioffbn ) ) * abs( a( ioffak ) ) )
   40          CONTINUE
   50       CONTINUE
         END IF
*
         ioffc = ic + ibeg - 1 + ( jc + j - 2 ) * ldc
*
         DO 100 i = ibeg, iend
            ct( i ) = ct( i ) + beta * c( ioffc )
            g( i ) = g( i ) + abs( beta )*abs( c( ioffc ) )
            c( ioffc ) = ct( i )
            ioffc = ioffc + 1
  100    CONTINUE
*
*        Compute the error ratio for this result.
*
         err  = zero
         info = 0
         ldpc = descc( lld_ )
         ioffc = ic + ( jc + j - 2 ) * ldc
         CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
     $                    iic, jjc, icrow, iccol )
         icurrow = icrow
         rowrep  = ( icrow.EQ.-1 )
         colrep  = ( iccol.EQ.-1 )
*
         IF( mycol.EQ.iccol .OR. colrep ) THEN
*
            ibb = descc( imb_ ) - ic + 1
            IF( ibb.LE.0 )
     $         ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
            ibb = min( ibb, n )
            in = ic + ibb - 1
*
            DO 110 i = ic, in
*
               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                  erri = abs( pc( iic+(jjc-1)*ldpc ) -
     $                        c( ioffc ) ) / eps
                  IF( g( i-ic+1 ).NE.zero )
     $               erri = erri / g( i-ic+1 )
                  err = max( err, erri )
                  IF( err*sqrt( eps ).GE.one )
     $               info = 1
                  iic = iic + 1
               END IF
*
               ioffc = ioffc + 1
*
  110       CONTINUE
*
            icurrow = mod( icurrow+1, nprow )
*
            DO 130 i = in+1, ic+n-1, descc( mb_ )
               ibb = min( ic+n-i, descc( mb_ ) )
*
               DO 120 kk = 0, ibb-1
*
                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
                     erri = abs( pc( iic+(jjc-1)*ldpc ) -
     $                           c( ioffc ) )/eps
                     IF( g( i+kk-ic+1 ).NE.zero )
     $                  erri = erri / g( i+kk-ic+1 )
                     err = max( err, erri )
                     IF( err*sqrt( eps ).GE.one )
     $                  info = 1
                     iic = iic + 1
                  END IF
*
                  ioffc = ioffc + 1
*
  120          CONTINUE
*
               icurrow = mod( icurrow+1, nprow )
*
  130       CONTINUE
*
         END IF
*
*        If INFO = 0, all results are at least half accurate.
*
         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
         CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $                 mycol )
         IF( info.NE.0 )
     $      GO TO 150
*
  140 CONTINUE
*
  150 CONTINUE
*
      RETURN
*
*     End of PDMMCH2
*
      END
      SUBROUTINE pdmmch3( UPLO, TRANS, M, N, ALPHA, A, IA, JA, DESCA,
     $                    BETA, C, PC, IC, JC, DESCC, ERR, INFO )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        TRANS, UPLO
      INTEGER            IA, IC, INFO, JA, JC, M, N
      DOUBLE PRECISION   ALPHA, BETA, ERR
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCC( * )
      DOUBLE PRECISION   A( * ), C( * ), PC( * )
*     ..
*
*  Purpose
*  =======
*
*  PDMMCH3 checks the results of the computational tests.
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESCA:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
*  CTXT_A  (global) DESCA( CTXT_  ) The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is distributed over.  The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA( M_     ) The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA( N_     ) The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA( IMB_   ) The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA( INB_   ) The  number  of columns of the upper
*                                   left   block   of   the   matrix  A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA( MB_    ) The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A rows of  A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA( NB_    ) The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA( RSRC_  ) The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA( CSRC_  ) The  process  column  over which the
*                                   first  column of  A  is distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA( LLD_   ) The  leading  dimension of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_NUMROC:
*  Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  UPLO    (global input) CHARACTER*1
*          On entry,  UPLO  specifies which part of C should contain the
*          result.
*
*  TRANS   (global input) CHARACTER*1
*          On entry,  TRANS  specifies  whether  the  matrix A has to be
*          transposed  or not  before computing the  matrix-matrix addi-
*          tion.
*
*  M       (global input) INTEGER
*          On entry, M specifies the number of rows of C.
*
*  N       (global input) INTEGER
*          On entry, N specifies the number of columns of C.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  A       (local input) DOUBLE PRECISION array
*          On entry, A is an array of  dimension  (DESCA( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PA.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  BETA    (global input) DOUBLE PRECISION
*          On entry, BETA specifies the scalar beta.
*
*  C       (local input/local output) DOUBLE PRECISION array
*          On entry, C is an array of  dimension  (DESCC( M_ ),*).  This
*          array contains a local copy of the initial entire matrix PC.
*
*  PC      (local input) DOUBLE PRECISION array
*          On entry, PC is an array of dimension (DESCC( LLD_ ),*). This
*          array contains the local pieces of the matrix PC.
*
*  IC      (global input) INTEGER
*          On entry, IC  specifies C's global row index, which points to
*          the beginning of the submatrix sub( C ).
*
*  JC      (global input) INTEGER
*          On entry, JC  specifies C's global column index, which points
*          to the beginning of the submatrix sub( C ).
*
*  DESCC   (global and local input) INTEGER array
*          On entry, DESCC  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix C.
*
*  ERR     (global output) DOUBLE PRECISION
*          On exit, ERR specifies the largest error in absolute value.
*
*  INFO    (global output) INTEGER
*          On exit, if INFO <> 0, the result is less than half accurate.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
     $                   RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D_INB = 2, dlen_ = 11,
     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLREP, LOWER, NOTRAN, ROWREP, UPPER
      INTEGER            I, ICCOL, ICROW, ICTXT, IIC, IOFFA, IOFFC, J,
     $                   JJC, LDA, LDC, LDPC, MYCOL, MYROW, NPCOL,
     $                   NPROW
      DOUBLE PRECISION   ERR0, ERRI, PREC
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, DGAMX2D, IGSUM2D, PB_INFOG2L,
     $                   pderraxpby
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           LSAME, PDLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
      ictxt = descc( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      prec   = pdlamch( ictxt, 'eps' )
*
      upper  = lsame( uplo,  'U' )
      lower  = lsame( uplo,  'L' )
      notran = lsame( trans, 'N' )
*
*     Compute expected result in C using data in A and C. This part of
*     the computation is performed by every process in the grid.
*
      info   = 0
      err    = zero
*
      lda    = max( 1, desca( m_   ) )
      ldc    = max( 1, descc( m_   ) )
      ldpc   = max( 1, descc( lld_ ) )
      rowrep = ( descc( rsrc_ ).EQ.-1 )
      colrep = ( descc( csrc_ ).EQ.-1 )
*
      IF( notran ) THEN
*
         DO 20 j = jc, jc + n - 1
*
            ioffc = ic + ( j  - 1          ) * ldc
            ioffa = ia + ( ja - 1 + j - jc ) * lda
*
            DO 10 i = ic, ic + m - 1
*
               IF( upper ) THEN
                  IF( ( j - jc ).GE.( i - ic ) ) THEN
                     CALL pderraxpby( erri, alpha, a( ioffa ), beta,
     $                                c( ioffc ), prec )
                  ELSE
                     erri = zero
                  END IF
               ELSE IF( lower ) THEN
                  IF( ( j - jc ).LE.( i - ic ) ) THEN
                     CALL pderraxpby( erri, alpha, a( ioffa ), beta,
     $                                c( ioffc ), prec )
                  ELSE
                     erri = zero
                  END IF
               ELSE
                  CALL pderraxpby( erri, alpha, a( ioffa ), beta,
     $                             c( ioffc ), prec )
               END IF
*
               CALL pb_infog2l( i, j, descc, nprow, npcol, myrow, mycol,
     $                          iic, jjc, icrow, iccol )
               IF( ( myrow.EQ.icrow .OR. rowrep ) .AND.
     $             ( mycol.EQ.iccol .OR. colrep ) ) THEN
                  err0 = abs( pc( iic+(jjc-1)*ldpc )-c( ioffc ) )
                  IF( err0.GT.erri )
     $               info = 1
                  err = max( err, err0 )
               END IF
*
               ioffa = ioffa + 1
               ioffc = ioffc + 1
*
   10       CONTINUE
*
   20    CONTINUE
*
      ELSE
*
         DO 40 j = jc, jc + n - 1
*
            ioffc = ic +              ( j  - 1 ) * ldc
            ioffa = ia + ( j - jc ) + ( ja - 1 ) * lda
*
            DO 30 i = ic, ic + m - 1
*
               IF( upper ) THEN
                  IF( ( j - jc ).GE.( i - ic ) ) THEN
                     CALL pderraxpby( erri, alpha, a( ioffa ), beta,
     $                               c( ioffc ), prec )
                  ELSE
                     erri = zero
                  END IF
               ELSE IF( lower ) THEN
                  IF( ( j - jc ).LE.( i - ic ) ) THEN
                     CALL pderraxpby( erri, alpha, a( ioffa ), beta,
     $                               c( ioffc ), prec )
                  ELSE
                     erri = zero
                  END IF
               ELSE
                  CALL pderraxpby( erri, alpha, a( ioffa ), beta,
     $                            c( ioffc ), prec )
               END IF
*
               CALL pb_infog2l( i, j, descc, nprow, npcol, myrow, mycol,
     $                          iic, jjc, icrow, iccol )
               IF( ( myrow.EQ.icrow .OR. rowrep ) .AND.
     $             ( mycol.EQ.iccol .OR. colrep ) ) THEN
                  err0 = abs( pc( iic+(jjc-1)*ldpc )-c( ioffc ) )
                  IF( err0.GT.erri )
     $               info = 1
                  err = max( err, err0 )
               END IF
*
               ioffc = ioffc + 1
               ioffa = ioffa + lda
*
   30       CONTINUE
*
   40    CONTINUE
*
      END IF
*
*     If INFO = 0, all results are at least half accurate.
*
      CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
      CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
     $              mycol )
*
      RETURN
*
*     End of PDMMCH3
*
      END
      SUBROUTINE pderraxpby( ERRBND, ALPHA, X, BETA, Y, PREC )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   ALPHA, BETA, ERRBND, PREC, X, Y
*     ..
*
*  Purpose
*  =======
*
*  PDERRAXPBY  serially  computes  y := beta*y + alpha * x and returns a
*  scaled relative acceptable error bound on the result.
*
*  Arguments
*  =========
*
*  ERRBND  (global output) DOUBLE PRECISION
*          On exit, ERRBND  specifies the scaled relative acceptable er-
*          ror bound.
*
*  ALPHA   (global input) DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.
*
*  X       (global input) DOUBLE PRECISION
*          On entry, X  specifies the scalar x to be scaled.
*
*  BETA    (global input) DOUBLE PRECISION
*          On entry, BETA specifies the scalar beta.
*
*  Y       (global input/global output) DOUBLE PRECISION
*          On entry,  Y  specifies  the scalar y to be added. On exit, Y
*          contains the resulting scalar y.
*
*  PREC    (global input) DOUBLE PRECISION
*          On entry, PREC specifies the machine precision.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University  of  Tennessee, Knoxville 37996, USA.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, TWO, ZERO
      PARAMETER          ( ONE = 1.0d+0, two = 2.0d+0,
     $                   zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   ADDBND, FACT, SUMPOS, SUMNEG, TMP
*     ..
*     .. Intrinsic Functions ..
*     ..
*     .. Executable Statements ..
*
      SUMPOS = zero
      sumneg = zero
      fact = one + two * prec
      addbnd = two * two * two * prec
*
      tmp = alpha * x
      IF( tmp.GE.zero ) THEN
         sumpos = sumpos + tmp * fact
      ELSE
         sumneg = sumneg - tmp * fact
      END IF
*
      tmp = beta * y
      IF( tmp.GE.zero ) THEN
         sumpos = sumpos + tmp * fact
      ELSE
         sumneg = sumneg - tmp * fact
      END IF
*
      y = ( beta * y ) + ( alpha * x )
*
      errbnd = addbnd * max( sumpos, sumneg )
*
      RETURN
*
*     End of PDERRAXPBY
*
      END
      DOUBLE PRECISION   FUNCTION pdlamch( ICTXT, CMACH )
*
*  -- PBLAS test routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*     .. Scalar Arguments ..
      CHARACTER*1        cmach
      INTEGER            ictxt
*     ..
*
*  Purpose
*  =======
*
*  PDLAMCH determines double precision machine parameters.
*
&#