pcsepsubtst.f

Go to the documentation of this file.

*
*
      SUBROUTINE pcsepsubtst( WKNOWN, JOBZ, RANGE, UPLO, N, VL, VU, IL,
     $                        IU, THRESH, ABSTOL, A, COPYA, Z, IA, JA,
     $                        DESCA, WIN, WNEW, IFAIL, ICLUSTR, GAP,
     $                        IPREPAD, IPOSTPAD, WORK, LWORK, RWORK,
     $                        LRWORK, LWORK1, IWORK, LIWORK, RESULT,
     $                        TSTNRM, QTQNRM, NOUT )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*     .. Scalar Arguments ..
      LOGICAL            WKNOWN
      CHARACTER          JOBZ, RANGE, UPLO
      INTEGER            IA, IL, IPOSTPAD, IPREPAD, IU, JA, LIWORK,
     $                   LRWORK, LWORK, LWORK1, N, NOUT, RESULT
      REAL               ABSTOL, QTQNRM, THRESH, TSTNRM, VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), ICLUSTR( * ), IFAIL( * ),
     $                   IWORK( * )
      REAL               GAP( * ), RWORK( * ), WIN( * ), WNEW( * )
      COMPLEX            A( * ), COPYA( * ), WORK( * ), Z( * )
*     ..
*
*  Purpose
*  =======
*
*  PCSEPSUBTST calls PCHEEVX and then tests the output of
*  PCHEEVX
*  If JOBZ = 'V' then the following two tests are performed:
*     |AQ -QL| / (abstol + eps * norm(A) ) < N*THRESH
*     |QT * Q - I| / eps * norm(A) < N*THRESH
*  If WKNOWN then
*     we check to make sure that the eigenvalues match expectations
*     i.e. |WIN - WNEW(1+IPREPAD)| / (eps * |WIN|) < THRESH
*     where WIN is the array of eigenvalues as computed by
*     PCHEEVX when eigenvectors are requested
*
*  Arguments
*  =========
*
*     NP = the number of rows local to a given process.
*     NQ = the number of columns local to a given process.
*
*  WKNOWN  (global input) INTEGER
*          .FALSE.:  WIN does not contain the eigenvalues
*          .TRUE.:   WIN does contain the eigenvalues
*
*  JOBZ    (global input) CHARACTER*1
*          Specifies whether or not to compute the eigenvectors:
*          = 'N':  Compute eigenvalues only.
*          = 'V':  Compute eigenvalues and eigenvectors.
*          Must be 'V' on first call to PCSEPSUBTST
*
*  RANGE   (global input) CHARACTER*1
*          = 'A': all eigenvalues will be found.
*          = 'V': all eigenvalues in the interval [VL,VU]
*                 will be found.
*          = 'I': the IL-th through IU-th eigenvalues will be found.
*          Must be 'A' on first call to PCSEPSUBTST
*
*  UPLO    (global input) CHARACTER*1
*          Specifies whether the upper or lower triangular part of the
*          Hermitian matrix A is stored:
*          = 'U':  Upper triangular
*          = 'L':  Lower triangular
*
*  N       (global input) INTEGER
*          Size of the matrix to be tested.  (global size)
*
*  VL      (global input) REAL
*          If RANGE='V', the lower bound of the interval to be searched
*          for eigenvalues.  Not referenced if RANGE = 'A' or 'I'.
*
*  VU      (global input) REAL
*          If RANGE='V', the upper bound of the interval to be searched
*          for eigenvalues.  Not referenced if RANGE = 'A' or 'I'.
*
*  IL      (global input) INTEGER
*          If RANGE='I', the index (from smallest to largest) of the
*          smallest eigenvalue to be returned.  IL >= 1.
*          Not referenced if RANGE = 'A' or 'V'.
*
*  IU      (global input) INTEGER
*          If RANGE='I', the index (from smallest to largest) of the
*          largest eigenvalue to be returned.  min(IL,N) <= IU <= N.
*          Not referenced if RANGE = 'A' or 'V'.
*
*  THRESH  (global input) REAL
*          A test will count as "failed" if the "error", computed as
*          described below, exceeds THRESH.  Note that the error
*          is scaled to be O(1), so THRESH should be a reasonably
*          small multiple of 1, e.g., 10 or 100.  In particular,
*          it should not depend on the precision (single vs. double)
*          or the size of the matrix.  It must be at least zero.
*
*  ABSTOL  (global input) REAL
*          The absolute tolerance for the eigenvalues. An
*          eigenvalue is considered to be located if it has
*          been determined to lie in an interval whose width
*          is "abstol" or less. If "abstol" is less than or equal
*          to zero, then ulp*|T| will be used, where |T| is
*          the 1-norm of the matrix. If eigenvectors are
*          desired later by inverse iteration ("PCSTEIN"),
*          "abstol" MUST NOT be bigger than ulp*|T|.
*
*  A       (local workspace) COMPLEX array
*          global dimension (N, N), local dimension (DESCA(DLEN_), NQ)
*            A is distributed in a block cyclic manner over both rows
*            and columns.
*            See PCHEEVX for a description of block cyclic layout.
*          The test matrix, which is then modified by PCHEEVX
*          A has already been padded front and back, use A(1+IPREPAD)
*
*  COPYA   (local input) COMPLEX array, dimension(N*N)
*          COPYA holds a copy of the original matrix A
*          identical in both form and content to A
*
*  Z       (local workspace) COMPLEX array, dim (N*N)
*          Z is distributed in the same manner as A
*          Z contains the eigenvector matrix
*          Z is used as workspace by the test routines
*          PCSEPCHK and PCSEPQTQ.
*          Z has already been padded front and back, use Z(1+IPREPAD)
*
*  IA      (global input) INTEGER
*          On entry, IA specifies the global row index of the submatrix
*          of the global matrix A, COPYA and Z to operate on.
*
*  JA      (global input) INTEGER
*          On entry, IA specifies the global column index of the submat
*          of the global matrix A, COPYA and Z to operate on.
*
*  DESCA   (global/local input) INTEGER array of dimension 8
*          The array descriptor for the matrix A, COPYA and Z.
*
*  WIN     (global input) REAL array, dimension (N)
*          If .not. WKNOWN, WIN is ignored on input
*          Otherwise, WIN() is taken as the standard by which the
*          eigenvalues are to be compared against.
*
*  WNEW    (global workspace)  REAL array, dimension (N)
*          The eigenvalues as copmuted by this call to PCHEEVX
*          If JOBZ <> 'V' or RANGE <> 'A' these eigenvalues are
*          compared against those in WIN().
*          WNEW has already been padded front and back,
*          use WNEW(1+IPREPAD)
*
*  IFAIL   (global output) INTEGER array, dimension (N)
*          If JOBZ = 'V', then on normal exit, the first M elements of
*          IFAIL are zero.  If INFO > 0 on exit, then IFAIL contains the
*          indices of the eigenvectors that failed to converge.
*          If JOBZ = 'N', then IFAIL is not referenced.
*          IFAIL has already been padded front and back,
*          use IFAIL(1+IPREPAD)
*
*  ICLUSTR (global workspace) integer array, dimension (2*NPROW*NPCOL)
*
*  GAP     (global workspace) REAL array,
*          dimension (NPROW*NPCOL)
*
*  WORK    (local workspace) COMPLEX array, dimension (LWORK)
*          WORK has already been padded front and back,
*          use WORK(1+IPREPAD)
*
*  LWORK   (local input) INTEGER
*          The actual length of the array WORK after padding.
*
*  RWORK   (local workspace) REAL array, dimension (LRWORK)
*          RWORK has already been padded front and back,
*          use RWORK(1+IPREPAD)
*
*  LRWORK   (local input) INTEGER
*          The actual length of the array RWORK after padding.
*
*  LWORK1  (local input) INTEGER
*          The amount of real workspace to pass to PCHEEVX
*
*  IWORK   (local workspace) INTEGER array, dimension (LIWORK)
*          IWORK has already been padded front and back,
*          use IWORK(1+IPREPAD)
*
*  LIWORK  (local input) INTEGER
*          The length of the array IWORK after padding.
*
*  RESULT  (global output) INTEGER
*          The result of this call to PCHEEVX
*          RESULT = -3   =>  This process did not participate
*          RESULT = 0    =>  All tests passed
*          RESULT = 1    =>  ONe or more tests failed
*
*  TSTNRM  (global output) REAL
*          |AQ- QL| / (ABSTOL+EPS*|A|)*N
*
*  QTQNRM  (global output) REAL
*          |QTQ -I| / N*EPS
*
*     .. Parameters ..
*
      INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
     $                   MB_, NB_, RSRC_, CSRC_, LLD_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
     $                   ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                   rsrc_ = 7, csrc_ = 8, lld_ = 9 )
      REAL               PADVAL, FIVE, NEGONE
      PARAMETER          ( PADVAL = 13.5285e+0, five = 5.0e+0,
     $                   negone = -1.0e+0 )
      COMPLEX            CPADVAL
      PARAMETER          ( CPADVAL = ( 13.989e+0, 1.93e+0 ) )
      INTEGER            IPADVAL
      parameter( ipadval = 927 )
*     ..
*     .. Local Scalars ..
      LOGICAL            MISSLARGEST, MISSSMALLEST
      INTEGER            I, IAM, INDIWRK, INFO, ISIZEHEEVX, ISIZESUBTST,
     $                   isizetst, j, m, maxeigs, maxil, maxiu, maxsize,
     $                   minil, mq, mycol, myil, myrow, nclusters, np,
     $                   npcol, nprow, nq, nz, oldil, oldiu, oldnz, res,
     $                   rsizechk, rsizeheevx, rsizeqtq, rsizesubtst,
     $                   rsizetst, sizeheevx, sizemqrleft, sizemqrright,
     $                   sizeqrf, sizesubtst, sizetms, sizetst, valsize,
     $                   vecsize, sizeheevd, rsizeheevd, isizeheevd
      REAL               EPS, EPSNORMA, ERROR, MAXERROR, MAXVU,
     $                   MINERROR, MINVL, NORMWIN, OLDVL, OLDVU, ORFAC,
     $                   SAFMIN
*     ..
*     .. Local Arrays ..
      INTEGER            DESCZ( DLEN_ ), DSEED( 4 ), ITMP( 2 )
*     ..
*     .. External Functions ..
*
      LOGICAL            LSAME
      INTEGER            NUMROC
      REAL               PCLANHE, PSLAMCH
      EXTERNAL           LSAME, NUMROC, PCLANHE, PSLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, clacpy, descinit, igamn2d,
     $                   igamx2d, pcchekpad, pcelset, pcfillpad,
     $                   pcheevx, pclasizeheevx, pclasizesep, pcsepchk,
     $                   pcsepqtq, pichekpad, pifillpad, pschekpad,
     $                   psfillpad, sgamn2d, sgamx2d, slboot, sltimer
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, mod
*     ..
*     .. Executable Statements ..
*       This is just to keep ftnchek happy
      IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
     $    rsrc_.LT.0 )RETURN
      CALL pclasizesep( desca, iprepad, ipostpad, sizemqrleft,
     $                  sizemqrright, sizeqrf, sizetms, rsizeqtq,
     $                  rsizechk, sizeheevx, rsizeheevx, isizeheevx,
     $                  sizeheevd, rsizeheevd, isizeheevd,
     $                  sizesubtst, rsizesubtst, isizesubtst, sizetst,
     $                  rsizetst, isizetst )
*
      tstnrm = negone
      qtqnrm = negone
      eps = pslamch( desca( ctxt_ ), 'Eps' )
      safmin = pslamch( desca( ctxt_ ), 'Safe min' )
*
      normwin = safmin / eps
      IF( n.GE.1 )
     $   normwin = max( abs( win( 1 ) ), abs( win( n ) ), normwin )
*
*     Make sure that we aren't using information from previous calls
*
      nz = -13
      oldnz = nz
      oldil = il
      oldiu = iu
      oldvl = vl
      oldvu = vu
*
      DO 10 i = 1, lwork1, 1
         rwork( i+iprepad ) = 14.3e+0
   10 CONTINUE
      DO 20 i = 1, liwork, 1
         iwork( i+iprepad ) = 14
   20 CONTINUE
      DO 30 i = 1, lwork, 1
         work( i+iprepad ) = ( 15.63e+0, 1.1e+0 )
   30 CONTINUE
*
      DO 40 i = 1, n
         wnew( i+iprepad ) = 3.14159e+0
   40 CONTINUE
*
      iclustr( 1+iprepad ) = 139
*
      IF( lsame( jobz, 'N' ) ) THEN
         maxeigs = 0
      ELSE
         IF( lsame( range, 'A' ) ) THEN
            maxeigs = n
         ELSE IF( lsame( range, 'I' ) ) THEN
            maxeigs = iu - il + 1
         ELSE
            minvl = vl - normwin*five*eps - abstol
            maxvu = vu + normwin*five*eps + abstol
            minil = 1
            maxiu = 0
            DO 50 i = 1, n
               IF( win( i ).LT.minvl )
     $            minil = minil + 1
               IF( win( i ).LE.maxvu )
     $            maxiu = maxiu + 1
   50       CONTINUE
*
            maxeigs = maxiu - minil + 1
         END IF
      END IF
*
*
      CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
     $               desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
     $               desca( ctxt_ ), desca( lld_ ), info )
*
      CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
      indiwrk = 1 + iprepad + nprow*npcol + 1
*
      iam = 1
      IF( myrow.EQ.0 .AND. mycol.EQ.0 )
     $   iam = 0
*
*     If this process is not involved in this test, bail out now
*
      result = -3
      IF( myrow.GE.nprow .OR. myrow.LT.0 )
     $   GO TO 160
      result = 0
*
*
*     DSEED is not used in this call to PCLASIZEHEEVX, the
*     following line just makes ftnchek happy.
*
      dseed( 1 ) = 1
*
      CALL pclasizeheevx( wknown, range, n, desca, vl, vu, il, iu,
     $                    dseed, win, maxsize, vecsize, valsize )
*
      np = numroc( n, desca( mb_ ), myrow, 0, nprow )
      nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
      mq = numroc( maxeigs, desca( nb_ ), mycol, 0, npcol )
*
      CALL clacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
     $             desca( lld_ ) )
*
      CALL pcfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
     $                ipostpad, cpadval )
*
      CALL pcfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
     $                ipostpad, cpadval+1.0e+0 )
*
      CALL psfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
     $                padval+2.0e+0 )
*
      CALL psfillpad( desca( ctxt_ ), nprow*npcol, 1, gap, nprow*npcol,
     $                iprepad, ipostpad, padval+3.0e+0 )
*
      CALL psfillpad( desca( ctxt_ ), lwork1, 1, rwork, lwork1, iprepad,
     $                ipostpad, padval+4.0e+0 )
*
      CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
     $                ipostpad, ipadval )
*
      CALL pifillpad( desca( ctxt_ ), n, 1, ifail, n, iprepad, ipostpad,
     $                ipadval )
*
      CALL pifillpad( desca( ctxt_ ), 2*nprow*npcol, 1, iclustr,
     $                2*nprow*npcol, iprepad, ipostpad, ipadval )
*
      CALL pcfillpad( desca( ctxt_ ), lwork, 1, work, lwork, iprepad,
     $                ipostpad, cpadval+4.1e+0 )
*
*     Make sure that PCHEEVX does not cheat (i.e. use answers
*     already computed.)
*
      DO 70 i = 1, n, 1
         DO 60 j = 1, maxeigs, 1
            CALL pcelset( z( 1+iprepad ), i, j, desca,
     $                    ( 13.0e+0, 1.34e+0 ) )
   60    CONTINUE
   70 CONTINUE
*
      orfac = -1.0e+0
*
      CALL slboot
      CALL sltimer( 1 )
      CALL sltimer( 6 )
      CALL pcheevx( jobz, range, uplo, n, a( 1+iprepad ), ia, ja, desca,
     $              vl, vu, il, iu, abstol, m, nz, wnew( 1+iprepad ),
     $              orfac, z( 1+iprepad ), ia, ja, desca,
     $              work( 1+iprepad ), sizeheevx, rwork( 1+iprepad ),
     $              lwork1, iwork( 1+iprepad ), liwork,
     $              ifail( 1+iprepad ), iclustr( 1+iprepad ),
     $              gap( 1+iprepad ), info )
      CALL sltimer( 6 )
      CALL sltimer( 1 )
*
      IF( thresh.LE.0 ) THEN
         result = 0
      ELSE
         CALL pcchekpad( desca( ctxt_ ), 'PCHEEVX-A', np, nq, a,
     $                   desca( lld_ ), iprepad, ipostpad, cpadval )
*
         CALL pcchekpad( descz( ctxt_ ), 'PCHEEVX-Z', np, mq, z,
     $                   descz( lld_ ), iprepad, ipostpad,
     $                   cpadval+1.0e+0 )
*
         CALL pschekpad( desca( ctxt_ ), 'PCHEEVX-WNEW', n, 1, wnew, n,
     $                   iprepad, ipostpad, padval+2.0e+0 )
*
         CALL pschekpad( desca( ctxt_ ), 'PCHEEVX-GAP', nprow*npcol, 1,
     $                   gap, nprow*npcol, iprepad, ipostpad,
     $                   padval+3.0e+0 )
*
         CALL pschekpad( desca( ctxt_ ), 'PCHEEVX-rWORK', lwork1, 1,
     $                   rwork, lwork1, iprepad, ipostpad,
     $                   padval+4.0e+0 )
*
         CALL pcchekpad( desca( ctxt_ ), 'PCHEEVX-WORK', lwork, 1, work,
     $                   lwork, iprepad, ipostpad, cpadval+4.1e+0 )
*
         CALL pichekpad( desca( ctxt_ ), 'PCHEEVX-IWORK', liwork, 1,
     $                   iwork, liwork, iprepad, ipostpad, ipadval )
*
         CALL pichekpad( desca( ctxt_ ), 'PCHEEVX-IFAIL', n, 1, ifail,
     $                   n, iprepad, ipostpad, ipadval )
*
         CALL pichekpad( desca( ctxt_ ), 'PCHEEVX-ICLUSTR',
     $                   2*nprow*npcol, 1, iclustr, 2*nprow*npcol,
     $                   iprepad, ipostpad, ipadval )
*
*
*     Since we now know the spectrum, we can potentially reduce MAXSIZE.
*
         IF( lsame( range, 'A' ) ) THEN
            CALL pclasizeheevx( .true., range, n, desca, vl, vu, il, iu,
     $                          dseed, wnew( 1+iprepad ), maxsize,
     $                          vecsize, valsize )
         END IF
*
*
*     Check INFO
*
*
*     Make sure that all processes return the same value of INFO
*
         itmp( 1 ) = info
         itmp( 2 ) = info
*
         CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
     $                 -1, -1, 0 )
         CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
     $                 1, -1, -1, 0 )
*
*
         IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = * )
     $         'Different processes return different INFO'
            result = 1
         ELSE IF( mod( info, 2 ).EQ.1 .OR. info.GT.7 .OR. info.LT.0 )
     $             THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9999 )info
            result = 1
         ELSE IF( mod( info / 2, 2 ).EQ.1 .AND. lwork1.GE.maxsize ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9996 )info
            result = 1
         ELSE IF( mod( info / 4, 2 ).EQ.1 .AND. lwork1.GE.vecsize ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9996 )info
            result = 1
         END IF
*
*
         IF( lsame( jobz, 'V' ) .AND. ( iclustr( 1+iprepad ).NE.
     $       0 ) .AND. ( mod( info / 2, 2 ).NE.1 ) ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9995 )
            result = 1
         END IF
*
*     Check M
*
         IF( ( m.LT.0 ) .OR. ( m.GT.n ) ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9994 )
            result = 1
         ELSE IF( lsame( range, 'A' ) .AND. ( m.NE.n ) ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9993 )
            result = 1
         ELSE IF( lsame( range, 'I' ) .AND. ( m.NE.iu-il+1 ) ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9992 )
            result = 1
         ELSE IF( lsame( jobz, 'V' ) .AND.
     $            ( .NOT.( lsame( range, 'V' ) ) ) .AND. ( m.NE.nz ) )
     $             THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9991 )
            result = 1
         END IF
*
*        Check NZ
*
         IF( lsame( jobz, 'V' ) ) THEN
            IF( lsame( range, 'V' ) ) THEN
               IF( nz.GT.m ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9990 )
                  result = 1
               END IF
               IF( nz.LT.m .AND. mod( info / 4, 2 ).NE.1 ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9989 )
                  result = 1
               END IF
            ELSE
               IF( nz.NE.m ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9988 )
                  result = 1
               END IF
            END IF
         END IF
         IF( result.EQ.0 ) THEN
*
*     Make sure that all processes return the same # of eigenvalues
*
            itmp( 1 ) = m
            itmp( 2 ) = m
*
            CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
     $                    -1, -1, 0 )
            CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
     $                    1, 1, -1, -1, 0 )
*
            IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
               IF( iam.EQ.0 )
     $            WRITE( nout, fmt = 9987 )
               result = 1
            ELSE
*
*     Make sure that different processes return the same eigenvalues
*
               DO 80 i = 1, m
                  rwork( i ) = wnew( i+iprepad )
                  rwork( i+m ) = wnew( i+iprepad )
   80          CONTINUE
*
               CALL sgamn2d( desca( ctxt_ ), 'a', ' ', m, 1, rwork, m,
     $                       1, 1, -1, -1, 0 )
               CALL sgamx2d( desca( ctxt_ ), 'a', ' ', m, 1,
     $                       rwork( 1+m ), m, 1, 1, -1, -1, 0 )
*
               DO 90 i = 1, m
*
                  IF( result.EQ.0 .AND. ( abs( rwork( i )-rwork( m+
     $                i ) ).GT.five*eps*abs( rwork( i ) ) ) ) THEN
                     IF( iam.EQ.0 )
     $                  WRITE( nout, fmt = 9986 )
                     result = 1
                  END IF
   90          CONTINUE
            END IF
         END IF
*
*     Make sure that all processes return the same # of clusters
*
         IF( lsame( jobz, 'V' ) ) THEN
            nclusters = 0
            DO 100 i = 0, nprow*npcol - 1
               IF( iclustr( 1+iprepad+2*i ).EQ.0 )
     $            GO TO 110
               nclusters = nclusters + 1
  100       CONTINUE
  110       CONTINUE
            itmp( 1 ) = nclusters
            itmp( 2 ) = nclusters
*
            CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
     $                    -1, -1, 0 )
            CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
     $                    1, 1, -1, -1, 0 )
*
            IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
               IF( iam.EQ.0 )
     $            WRITE( nout, fmt = 9985 )
               result = 1
            ELSE
*
*     Make sure that different processes return the same clusters
*
               DO 120 i = 1, nclusters
                  iwork( indiwrk+i ) = iclustr( i+iprepad )
                  iwork( indiwrk+i+nclusters ) = iclustr( i+iprepad )
  120          CONTINUE
               CALL igamn2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
     $                       iwork( indiwrk+1 ), nclusters*2+1, 1, 1,
     $                       -1, -1, 0 )
               CALL igamx2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
     $                       iwork( indiwrk+1+nclusters ),
     $                       nclusters*2+1, 1, 1, -1, -1, 0 )
*
*
               DO 130 i = 1, nclusters
                  IF( result.EQ.0 .AND. iwork( indiwrk+i ).NE.
     $                iwork( indiwrk+nclusters+i ) ) THEN
                     IF( iam.EQ.0 )
     $                  WRITE( nout, fmt = 9984 )
                     result = 1
                  END IF
  130          CONTINUE
*
               IF( iclustr( 1+iprepad+nclusters*2 ).NE.0 ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9983 )
                  result = 1
               END IF
            END IF
         END IF
*
*
         CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1,
     $                 -1, -1, 0 )
         IF( result.NE.0 )
     $      GO TO 160
*
*     Compute eps * norm(A)
*
         IF( n.EQ.0 ) THEN
            epsnorma = eps
         ELSE
            epsnorma = pclanhe( 'I', uplo, n, copya, ia, ja, desca,
     $                 rwork )*eps
         END IF
*
*     Note that a couple key variables get redefined in PCSEPCHK
*     as described by this table:
*
*     PCSEPTST name         PCSEPCHK name
*     -------------         -------------
*     COPYA                 A
*     Z                     Q
*     A                     C
*
*
         IF( lsame( jobz, 'V' ) ) THEN
*
*     Perform the |AQ - QE| test
*
            CALL psfillpad( desca( ctxt_ ), rsizechk, 1, rwork,
     $                      rsizechk, iprepad, ipostpad, 4.3e+0 )
*
            CALL pcsepchk( n, nz, copya, ia, ja, desca,
     $                     max( abstol+epsnorma, safmin ), thresh,
     $                     z( 1+iprepad ), ia, ja, descz,
     $                     a( 1+iprepad ), ia, ja, desca,
     $                     wnew( 1+iprepad ), rwork( 1+iprepad ),
     $                     rsizechk, tstnrm, res )
*
            CALL pschekpad( desca( ctxt_ ), 'PCSEPCHK-rWORK', rsizechk,
     $                      1, rwork, rsizechk, iprepad, ipostpad,
     $                      4.3e+0 )
*
            IF( res.NE.0 )
     $         result = 1
*
*     Perform the |QTQ - I| test
*
            CALL psfillpad( desca( ctxt_ ), rsizeqtq, 1, rwork,
     $                      rsizeqtq, iprepad, ipostpad, 4.3e+0 )
*
*
            CALL pcsepqtq( n, nz, thresh, z( 1+iprepad ), ia, ja, descz,
     $                     a( 1+iprepad ), ia, ja, desca,
     $                     iwork( 1+iprepad+1 ), iclustr( 1+iprepad ),
     $                     gap( 1+iprepad ), rwork( iprepad+1 ),
     $                     rsizeqtq, qtqnrm, info, res )
*
            CALL pschekpad( desca( ctxt_ ), 'PCSEPQTQ-rWORK', rsizeqtq,
     $                      1, rwork, rsizeqtq, iprepad, ipostpad,
     $                      4.3e+0 )
*
            IF( res.NE.0 )
     $         result = 1
*
            IF( info.NE.0 ) THEN
               IF( iam.EQ.0 )
     $            WRITE( nout, fmt = 9998 )info
               result = 1
            END IF
         END IF
*
*     Check to make sure that we have the right eigenvalues
*
         IF( wknown ) THEN
*
*        Set up MYIL if necessary
*
            myil = il
*
            IF( lsame( range, 'V' ) ) THEN
               myil = 1
               minil = 1
               maxil = n - m + 1
            ELSE
               IF( lsame( range, 'A' ) ) THEN
                  myil = 1
               END IF
               minil = myil
               maxil = myil
            END IF
*
*     Find the largest difference between the computed
*     and expected eigenvalues
*
            minerror = normwin
*
            DO 150 myil = minil, maxil
               maxerror = 0
*
*     Make sure that we aren't skipping any important eigenvalues
*
               misssmallest = .true.
               IF( .NOT.lsame( range, 'V' ) .OR. ( myil.EQ.1 ) )
     $            misssmallest = .false.
               IF( misssmallest .AND. ( win( myil-1 ).LT.vl+normwin*
     $             five*thresh*eps ) )misssmallest = .false.
               misslargest = .true.
               IF( .NOT.lsame( range, 'V' ) .OR. ( myil.EQ.maxil ) )
     $            misslargest = .false.
               IF( misslargest .AND. ( win( myil+m ).GT.vu-normwin*five*
     $             thresh*eps ) )misslargest = .false.
               IF( .NOT.misssmallest ) THEN
                  IF( .NOT.misslargest ) THEN
*
*     Make sure that the eigenvalues that we report are OK
*
                     DO 140 i = 1, m
                        error = abs( win( i+myil-1 )-wnew( i+iprepad ) )
                        maxerror = max( maxerror, error )
  140                CONTINUE
*
                     minerror = min( maxerror, minerror )
                  END IF
               END IF
  150       CONTINUE
*
*
*           If JOBZ = 'V' and RANGE='A', we might be comparing
*           against our estimate of what the eigenvalues ought to
*           be, rather than comparing against what PxHEEVX computed
*           last time around, so we have to be more generous.
*
            IF( lsame( jobz, 'V' ) .AND. lsame( range, 'A' ) ) THEN
               IF( minerror.GT.normwin*five*five*thresh*eps ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9997 )minerror, normwin
                  result = 1
               END IF
            ELSE
               IF( minerror.GT.normwin*five*thresh*eps ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9997 )minerror, normwin
                  result = 1
               END IF
            END IF
         END IF
*
*
*     Make sure that the IL, IU, VL and VU were not altered
*
         IF( il.NE.oldil .OR. iu.NE.oldiu .OR. vl.NE.oldvl .OR. vu.NE.
     $       oldvu ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9982 )
            result = 1
         END IF
*
         IF( lsame( jobz, 'N' ) .AND. ( nz.NE.oldnz ) ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9981 )
            result = 1
         END IF
*
      END IF
*
*     All processes should report the same result
*
      CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
     $              -1, 0 )
*
  160 CONTINUE
*
*
      RETURN
*
 9999 FORMAT( 'PCHEEVX returned INFO=', i7 )
 9998 FORMAT( 'PCSEPQTQ returned INFO=', i7 )
 9997 FORMAT( 'PCSEPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
 9996 FORMAT( 'PCHEEVX returned INFO=', i7,
     $      ' despite adequate workspace' )
 9995 FORMAT( .NE..NE.'ICLUSTR(1)0 but mod(INFO/2,2)1' )
 9994 FORMAT( 'M not in the range 0 to N' )
 9993 FORMAT( 'M not equal to N' )
 9992 FORMAT( 'M not equal to IU-IL+1' )
 9991 FORMAT( 'M not equal to NZ' )
 9990 FORMAT( 'NZ > M' )
 9989 FORMAT( 'NZ < M' )
 9988 FORMAT( 'NZ not equal to M' )
 9987 FORMAT( 'Different processes return different values for M' )
 9986 FORMAT( 'Different processes return different eigenvalues' )
 9985 FORMAT( 'Different processes return ',
     $      'different numbers of clusters' )
 9984 FORMAT( 'Different processes return different clusters' )
 9983 FORMAT( 'ICLUSTR not zero terminated' )
 9982 FORMAT( 'IL, IU, VL or VU altered by PCHEEVX' )
 9981 FORMAT( 'NZ altered by PCHEEVX with JOBZ=N' )
*
*     End of PCSEPSUBTST
*
      END