SUBROUTINE PDSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, ORFAC, Z, IZ,
     $                    JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, IFAIL,
     $                    ICLUSTR, GAP, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*     .. Scalar Arguments ..
      INTEGER            INFO, IZ, JZ, LIWORK, LWORK, M, N
      DOUBLE PRECISION   ORFAC
*     ..
*     .. Array Arguments ..
      INTEGER            DESCZ( * ), IBLOCK( * ), ICLUSTR( * ),
     $                   IFAIL( * ), ISPLIT( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), GAP( * ), W( * ), WORK( * ),
     $                   Z( * )
*     ..
*
*  Purpose
*  =======
*
*  PDSTEIN computes the eigenvectors of a symmetric tridiagonal matrix
*  in parallel, using inverse iteration. The eigenvectors found
*  correspond to user specified eigenvalues. PDSTEIN does not
*  orthogonalize vectors that are on different processes. The extent
*  of orthogonalization is controlled by the input parameter LWORK.
*  Eigenvectors that are to be orthogonalized are computed by the same
*  process. PDSTEIN decides on the allocation of work among the
*  processes and then calls DSTEIN2 (modified LAPACK routine) on each
*  individual process. If insufficient workspace is allocated, the
*  expected orthogonalization may not be done.
*
*  Note : If the eigenvectors obtained are not orthogonal, increase
*         LWORK and run the code again.
*
*  Notes
*  =====
*
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension r x c.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the r processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the c processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*
*  Arguments
*  =========
*
*          P = NPROW * NPCOL is the total number of processes
*
*  N       (global input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*  D       (global input) DOUBLE PRECISION array, dimension (N)
*          The n diagonal elements of the tridiagonal matrix T.
*
*  E       (global input) DOUBLE PRECISION array, dimension (N-1)
*          The (n-1) off-diagonal elements of the tridiagonal matrix T.
*
*  M       (global input) INTEGER
*          The total number of eigenvectors to be found. 0 <= M <= N.
*
*  W       (global input/global output) DOUBLE PRECISION array, dim (M)
*          On input, the first M elements of W contain all the
*          eigenvalues for which eigenvectors are to be computed. The
*          eigenvalues should be grouped by split-off block and ordered
*          from smallest to largest within the block (The output array
*          W from PDSTEBZ with ORDER='b' is expected here). This
*          array should be replicated on all processes.
*          On output, the first M elements contain the input
*          eigenvalues in ascending order.
*
*          Note : To obtain orthogonal vectors, it is best if
*          eigenvalues are computed to highest accuracy ( this can be
*          done by setting ABSTOL to the underflow threshold =
*          DLAMCH('U') --- ABSTOL is an input parameter
*          to PDSTEBZ )
*
*  IBLOCK  (global input) INTEGER array, dimension (N)
*          The submatrix indices associated with the corresponding
*          eigenvalues in W -- 1 for eigenvalues belonging to the
*          first submatrix from the top, 2 for those belonging to
*          the second submatrix, etc. (The output array IBLOCK
*          from PDSTEBZ is expected here).
*
*  ISPLIT  (global input) INTEGER array, dimension (N)
*          The splitting points, at which T breaks up into submatrices.
*          The first submatrix consists of rows/columns 1 to ISPLIT(1),
*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
*          etc., and the NSPLIT-th consists of rows/columns
*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N (The output array
*          ISPLIT from PDSTEBZ is expected here.)
*
*  ORFAC   (global input) DOUBLE PRECISION
*          ORFAC specifies which eigenvectors should be orthogonalized.
*          Eigenvectors that correspond to eigenvalues which are within
*          ORFAC*||T|| of each other are to be orthogonalized.
*          However, if the workspace is insufficient (see LWORK), this
*          tolerance may be decreased until all eigenvectors to be
*          orthogonalized can be stored in one process.
*          No orthogonalization will be done if ORFAC equals zero.
*          A default value of 10^-3 is used if ORFAC is negative.
*          ORFAC should be identical on all processes.
*
*  Z       (local output) DOUBLE PRECISION array,
*          dimension (DESCZ(DLEN_), N/npcol + NB)
*          Z contains the computed eigenvectors associated with the
*          specified eigenvalues. Any vector which fails to converge is
*          set to its current iterate after MAXITS iterations ( See
*          DSTEIN2 ).
*          On output, Z is distributed across the P processes in block
*          cyclic format.
*
*  IZ      (global input) INTEGER
*          Z's global row index, which points to the beginning of the
*          submatrix which is to be operated on.
*
*  JZ      (global input) INTEGER
*          Z's global column index, which points to the beginning of
*          the submatrix which is to be operated on.
*
*  DESCZ   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix Z.
*
*  WORK    (local workspace/global output) DOUBLE PRECISION array,
*          dimension ( LWORK )
*          On output, WORK(1) gives a lower bound on the
*          workspace ( LWORK ) that guarantees the user desired
*          orthogonalization (see ORFAC).
*          Note that this may overestimate the minimum workspace needed.
*
*  LWORK   (local input) integer
*          LWORK  controls the extent of orthogonalization which can be
*          done. The number of eigenvectors for which storage is
*          allocated on each process is
*                NVEC = floor(( LWORK- max(5*N,NP00*MQ00) )/N).
*          Eigenvectors corresponding to eigenvalue clusters of size
*          NVEC - ceil(M/P) + 1 are guaranteed to be orthogonal ( the
*          orthogonality is similar to that obtained from DSTEIN2).
*          Note : LWORK  must be no smaller than:
*               max(5*N,NP00*MQ00) + ceil(M/P)*N,
*          and should have the same input value on all processes.
*          It is the minimum value of LWORK input on different processes
*          that is significant.
*
*          If LWORK = -1, then LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  IWORK   (local workspace/global output) INTEGER array,
*          dimension ( 3*N+P+1 )
*          On return, IWORK(1) contains the amount of integer workspace
*          required.
*          On return, the IWORK(2) through IWORK(P+2) indicate
*          the eigenvectors computed by each process. Process I computes
*          eigenvectors indexed IWORK(I+2)+1 thru' IWORK(I+3).
*
*  LIWORK  (local input) INTEGER
*          Size of array IWORK.  Must be >= 3*N + P + 1
*
*          If LIWORK = -1, then LIWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  IFAIL   (global output) integer array, dimension (M)
*          On normal exit, all elements of IFAIL are zero.
*          If one or more eigenvectors fail to converge after MAXITS
*          iterations (as in DSTEIN), then INFO > 0 is returned.
*          If mod(INFO,M+1)>0, then
*          for I=1 to mod(INFO,M+1), the eigenvector
*          corresponding to the eigenvalue W(IFAIL(I)) failed to
*          converge ( W refers to the array of eigenvalues on output ).
*
*  ICLUSTR (global output) integer array, dimension (2*P)
*          This output array contains indices of eigenvectors
*          corresponding to a cluster of eigenvalues that could not be
*          orthogonalized due to insufficient workspace (see LWORK,
*          ORFAC and INFO). Eigenvectors corresponding to clusters of
*          eigenvalues indexed ICLUSTR(2*I-1) to ICLUSTR(2*I), I = 1 to
*          INFO/(M+1), could not be orthogonalized due to lack of
*          workspace. Hence the eigenvectors corresponding to these
*          clusters may not be orthogonal. ICLUSTR is a zero terminated
*          array --- ( ICLUSTR(2*K).NE.0 .AND. ICLUSTR(2*K+1).EQ.0 )
*          if and only if K is the number of clusters.
*
*  GAP     (global output) DOUBLE PRECISION array, dimension (P)
*          This output array contains the gap between eigenvalues whose
*          eigenvectors could not be orthogonalized. The INFO/M output
*          values in this array correspond to the INFO/(M+1) clusters
*          indicated by the array ICLUSTR. As a result, the dot product
*          between eigenvectors corresponding to the I^th cluster may be
*          as high as ( O(n)*macheps ) / GAP(I).
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*          < 0 :  if INFO = -I, the I-th argument had an illegal value
*          > 0 :  if mod(INFO,M+1) = I, then I eigenvectors failed to
*                 converge in MAXITS iterations. Their indices are
*                 stored in the array IFAIL.
*                 if INFO/(M+1) = I, then eigenvectors corresponding to
*                 I clusters of eigenvalues could not be orthogonalized
*                 due to insufficient workspace. The indices of the
*                 clusters are stored in the array ICLUSTR.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, MAX, MIN, MOD
*     ..
*     .. External Functions ..
      INTEGER            ICEIL, NUMROC
      EXTERNAL           ICEIL, NUMROC
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DGEBR2D, DGEBS2D,
     $                   DLASRT2, DSTEIN2, IGAMN2D, IGEBR2D, IGEBS2D,
     $                   PCHK1MAT, PDLAEVSWP, PXERBLA
*     ..
*     .. 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 )
      DOUBLE PRECISION   ZERO, NEGONE, ODM1, FIVE, ODM3, ODM18
      PARAMETER          ( ZERO = 0.0D+0, NEGONE = -1.0D+0,
     $                   ODM1 = 1.0D-1, FIVE = 5.0D+0, ODM3 = 1.0D-3,
     $                   ODM18 = 1.0D-18 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, SORTED
      INTEGER            B1, BN, BNDRY, CLSIZ, COL, I, IFIRST, IINFO,
     $                   ILAST, IM, INDRW, ITMP, J, K, LGCLSIZ, LLWORK,
     $                   LOAD, LOCINFO, MAXVEC, MQ00, MYCOL, MYROW,
     $                   NBLK, NERR, NEXT, NP00, NPCOL, NPROW, NVS,
     $                   OLNBLK, P, ROW, SELF, TILL, TOTERR
      DOUBLE PRECISION   DIFF, MINGAP, ONENRM, ORGFAC, ORTOL, TMPFAC
*     ..
*     .. Local Arrays ..
      INTEGER            IDUM1( 1 ), IDUM2( 1 )
*     ..
*     .. 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 BLACS_GRIDINFO( DESCZ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
      SELF = MYROW*NPCOL + MYCOL
*
*     Make sure that we belong to this context (before calling PCHK1MAT)
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -( 1200+CTXT_ )
      ELSE
*
*        Make sure that NPROW>0 and NPCOL>0 before calling NUMROC
*
         CALL CHK1MAT( N, 1, N, 1, IZ, JZ, DESCZ, 12, INFO )
         IF( INFO.EQ.0 ) THEN
*
*           Now we know that our context is good enough to
*           perform the rest of the checks
*
            NP00 = NUMROC( N, DESCZ( MB_ ), 0, 0, NPROW )
            MQ00 = NUMROC( M, DESCZ( NB_ ), 0, 0, NPCOL )
            P = NPROW*NPCOL
*
*           Compute the maximum number of vectors per process
*
            LLWORK = LWORK
            CALL IGAMN2D( DESCZ( CTXT_ ), 'A', ' ', 1, 1, LLWORK, 1, 1,
     $                    1, -1, -1, -1 )
            INDRW = MAX( 5*N, NP00*MQ00 )
            IF( N.NE.0 )
     $         MAXVEC = ( LLWORK-INDRW ) / N
            LOAD = ICEIL( M, P )
            IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
               TMPFAC = ORFAC
               CALL DGEBS2D( DESCZ( CTXT_ ), 'ALL', ' ', 1, 1, TMPFAC,
     $                       1 )
            ELSE
               CALL DGEBR2D( DESCZ( CTXT_ ), 'ALL', ' ', 1, 1, TMPFAC,
     $                       1, 0, 0 )
            END IF
*
            LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
            IF( M.LT.0 .OR. M.GT.N ) THEN
               INFO = -4
            ELSE IF( MAXVEC.LT.LOAD .AND. .NOT.LQUERY ) THEN
               INFO = -14
            ELSE IF( LIWORK.LT.3*N+P+1 .AND. .NOT.LQUERY ) THEN
               INFO = -16
            ELSE
               DO 10 I = 2, M
                  IF( IBLOCK( I ).LT.IBLOCK( I-1 ) ) THEN
                     INFO = -6
                     GO TO 20
                  END IF
                  IF( IBLOCK( I ).EQ.IBLOCK( I-1 ) .AND. W( I ).LT.
     $                W( I-1 ) ) THEN
                     INFO = -5
                     GO TO 20
                  END IF
   10          CONTINUE
   20          CONTINUE
               IF( INFO.EQ.0 ) THEN
                  IF( ABS( TMPFAC-ORFAC ).GT.FIVE*ABS( TMPFAC ) )
     $               INFO = -8
               END IF
            END IF
*
         END IF
         IDUM1( 1 ) = M
         IDUM2( 1 ) = 4
         CALL PCHK1MAT( N, 1, N, 1, IZ, JZ, DESCZ, 12, 1, IDUM1, IDUM2,
     $                  INFO )
         WORK( 1 ) = DBLE( MAX( 5*N, NP00*MQ00 )+ICEIL( M, P )*N )
         IWORK( 1 ) = 3*N + P + 1
      END IF
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( DESCZ( CTXT_ ), 'PDSTEIN', -INFO )
         RETURN
      ELSE IF( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) THEN
         RETURN
      END IF
*
      DO 30 I = 1, M
         IFAIL( I ) = 0
   30 CONTINUE
      DO 40 I = 1, P + 1
         IWORK( I ) = 0
   40 CONTINUE
      DO 50 I = 1, P
         GAP( I ) = NEGONE
         ICLUSTR( 2*I-1 ) = 0
         ICLUSTR( 2*I ) = 0
   50 CONTINUE
*
*
*     Quick return if possible
*
      IF( N.EQ.0 .OR. M.EQ.0 )
     $   RETURN
*
      IF( ORFAC.GE.ZERO ) THEN
         TMPFAC = ORFAC
      ELSE
         TMPFAC = ODM3
      END IF
      ORGFAC = TMPFAC
*
*     Allocate the work among the processes
*
      ILAST = M / LOAD
      IF( MOD( M, LOAD ).EQ.0 )
     $   ILAST = ILAST - 1
      OLNBLK = -1
      NVS = 0
      NEXT = 1
      IM = 0
      ONENRM = ZERO
      DO 100 I = 0, ILAST - 1
         NEXT = NEXT + LOAD
         J = NEXT - 1
         IF( J.GT.NVS ) THEN
            NBLK = IBLOCK( NEXT )
            IF( NBLK.EQ.IBLOCK( NEXT-1 ) .AND. NBLK.NE.OLNBLK ) THEN
*
*              Compute orthogonalization criterion
*
               IF( NBLK.EQ.1 ) THEN
                  B1 = 1
               ELSE
                  B1 = ISPLIT( NBLK-1 ) + 1
               END IF
               BN = ISPLIT( NBLK )
*
               ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
               ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
               DO 60 J = B1 + 1, BN - 1
                  ONENRM = MAX( ONENRM, ABS( D( J ) )+ABS( E( J-1 ) )+
     $                     ABS( E( J ) ) )
   60          CONTINUE
               OLNBLK = NBLK
            END IF
            TILL = NVS + MAXVEC
   70       CONTINUE
            J = NEXT - 1
            IF( TMPFAC.GT.ODM18 ) THEN
               ORTOL = TMPFAC*ONENRM
               DO 80 J = NEXT - 1, MIN( TILL, M-1 )
                  IF( IBLOCK( J+1 ).NE.IBLOCK( J ) .OR. W( J+1 )-
     $                W( J ).GE.ORTOL ) THEN
                     GO TO 90
                  END IF
   80          CONTINUE
               IF( J.EQ.M .AND. TILL.GE.M )
     $            GO TO 90
               TMPFAC = TMPFAC*ODM1
               GO TO 70
            END IF
   90       CONTINUE
            J = MIN( J, TILL )
         END IF
         IF( SELF.EQ.I )
     $      IM = MAX( 0, J-NVS )
*
         IWORK( I+1 ) = NVS
         NVS = MAX( J, NVS )
  100 CONTINUE
      IF( SELF.EQ.ILAST )
     $   IM = M - NVS
      IWORK( ILAST+1 ) = NVS
      DO 110 I = ILAST + 2, P + 1
         IWORK( I ) = M
  110 CONTINUE
*
      CLSIZ = 1
      LGCLSIZ = 1
      ILAST = 0
      NBLK = 0
      BNDRY = 2
      K = 1
      DO 140 I = 1, M
         IF( IBLOCK( I ).NE.NBLK ) THEN
            NBLK = IBLOCK( I )
            IF( NBLK.EQ.1 ) THEN
               B1 = 1
            ELSE
               B1 = ISPLIT( NBLK-1 ) + 1
            END IF
            BN = ISPLIT( NBLK )
*
            ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
            ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
            DO 120 J = B1 + 1, BN - 1
               ONENRM = MAX( ONENRM, ABS( D( J ) )+ABS( E( J-1 ) )+
     $                  ABS( E( J ) ) )
  120       CONTINUE
*
         END IF
         IF( I.GT.1 ) THEN
            DIFF = W( I ) - W( I-1 )
            IF( IBLOCK( I ).NE.IBLOCK( I-1 ) .OR. I.EQ.M .OR. DIFF.GT.
     $          ORGFAC*ONENRM ) THEN
               IFIRST = ILAST
               IF( I.EQ.M ) THEN
                  IF( IBLOCK( M ).NE.IBLOCK( M-1 ) .OR. DIFF.GT.ORGFAC*
     $                ONENRM ) THEN
                     ILAST = M - 1
                  ELSE
                     ILAST = M
                  END IF
               ELSE
                  ILAST = I - 1
               END IF
               CLSIZ = ILAST - IFIRST
               IF( CLSIZ.GT.1 ) THEN
                  IF( LGCLSIZ.LT.CLSIZ )
     $               LGCLSIZ = CLSIZ
                  MINGAP = ONENRM
  130             CONTINUE
                  IF( BNDRY.GT.P+1 )
     $               GO TO 150
                  IF( IWORK( BNDRY ).GT.IFIRST .AND. IWORK( BNDRY ).LT.
     $                ILAST ) THEN
                     MINGAP = MIN( W( IWORK( BNDRY )+1 )-
     $                        W( IWORK( BNDRY ) ), MINGAP )
                  ELSE IF( IWORK( BNDRY ).GE.ILAST ) THEN
                     IF( MINGAP.LT.ONENRM ) THEN
                        ICLUSTR( 2*K-1 ) = IFIRST + 1
                        ICLUSTR( 2*K ) = ILAST
                        GAP( K ) = MINGAP / ONENRM
                        K = K + 1
                     END IF
                     GO TO 140
                  END IF
                  BNDRY = BNDRY + 1
                  GO TO 130
               END IF
            END IF
         END IF
  140 CONTINUE
  150 CONTINUE
      INFO = ( K-1 )*( M+1 )
*
*     Call DSTEIN2 to find the eigenvectors
*
      CALL DSTEIN2( N, D, E, IM, W( IWORK( SELF+1 )+1 ),
     $              IBLOCK( IWORK( SELF+1 )+1 ), ISPLIT, ORGFAC,
     $              WORK( INDRW+1 ), N, WORK, IWORK( P+2 ),
     $              IFAIL( IWORK( SELF+1 )+1 ), LOCINFO )
*
*     Redistribute the eigenvector matrix to conform with the block
*     cyclic distribution of the input matrix
*
*
      DO 160 I = 1, M
         IWORK( P+1+I ) = I
  160 CONTINUE
*
      CALL DLASRT2( 'I', M, W, IWORK( P+2 ), IINFO )
*
      DO 170 I = 1, M
         IWORK( M+P+1+IWORK( P+1+I ) ) = I
  170 CONTINUE
*
*
      DO 180 I = 1, LOCINFO
         ITMP = IWORK( SELF+1 ) + I
         IFAIL( ITMP ) = IFAIL( ITMP ) + ITMP - I
         IFAIL( ITMP ) = IWORK( M+P+1+IFAIL( ITMP ) )
  180 CONTINUE
*
      DO 190 I = 1, K - 1
         ICLUSTR( 2*I-1 ) = IWORK( M+P+1+ICLUSTR( 2*I-1 ) )
         ICLUSTR( 2*I ) = IWORK( M+P+1+ICLUSTR( 2*I ) )
  190 CONTINUE
*
*
* Still need to apply the above permutation to IFAIL
*
*
      TOTERR = 0
      DO 210 I = 1, P
         IF( SELF.EQ.I-1 ) THEN
            CALL IGEBS2D( DESCZ( CTXT_ ), 'ALL', ' ', 1, 1, LOCINFO, 1 )
            IF( LOCINFO.NE.0 ) THEN
               CALL IGEBS2D( DESCZ( CTXT_ ), 'ALL', ' ', LOCINFO, 1,
     $                       IFAIL( IWORK( I )+1 ), LOCINFO )
               DO 200 J = 1, LOCINFO
                  IFAIL( TOTERR+J ) = IFAIL( IWORK( I )+J )
  200          CONTINUE
               TOTERR = TOTERR + LOCINFO
            END IF
         ELSE
*
            ROW = ( I-1 ) / NPCOL
            COL = MOD( I-1, NPCOL )
*
            CALL IGEBR2D( DESCZ( CTXT_ ), 'ALL', ' ', 1, 1, NERR, 1,
     $                    ROW, COL )
            IF( NERR.NE.0 ) THEN
               CALL IGEBR2D( DESCZ( CTXT_ ), 'ALL', ' ', NERR, 1,
     $                       IFAIL( TOTERR+1 ), NERR, ROW, COL )
               TOTERR = TOTERR + NERR
            END IF
         END IF
  210 CONTINUE
      INFO = INFO + TOTERR
*
*
      CALL PDLAEVSWP( N, WORK( INDRW+1 ), N, Z, IZ, JZ, DESCZ, IWORK,
     $                IWORK( M+P+2 ), WORK, INDRW )
*
      DO 220 I = 2, P
         IWORK( I ) = IWORK( M+P+1+IWORK( I ) )
  220 CONTINUE
*
*
*     Sort the IWORK array
*
*
  230 CONTINUE
      SORTED = .TRUE.
      DO 240 I = 2, P - 1
         IF( IWORK( I ).GT.IWORK( I+1 ) ) THEN
            ITMP = IWORK( I+1 )
            IWORK( I+1 ) = IWORK( I )
            IWORK( I ) = ITMP
            SORTED = .FALSE.
         END IF
  240 CONTINUE
      IF( .NOT.SORTED )
     $   GO TO 230
*
      DO 250 I = P + 1, 1, -1
         IWORK( I+1 ) = IWORK( I )
  250 CONTINUE
*
      WORK( 1 ) = ( LGCLSIZ+LOAD-1 )*N + INDRW
      IWORK( 1 ) = 3*N + P + 1
*
*     End of PDSTEIN
*
      END