pdstein.f

Go to the documentation of this file.

      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 )
.NE.            IF( N0 )
     $         MAXVEC = ( LLWORK-INDRW ) / N
            LOAD = ICEIL( M, P )
.EQ..AND..EQ.            IF( MYROW0  MYCOL0 ) 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
*
.EQ..OR..EQ.            LQUERY = ( LWORK-1  LIWORK-1 )
.LT..OR..GT.            IF( M0  MN ) THEN
               INFO = -4
.LT..AND..NOT.            ELSE IF( MAXVECLOAD  LQUERY ) THEN
               INFO = -14
.LT..AND..NOT.            ELSE IF( LIWORK3*N+P+1  LQUERY ) THEN
               INFO = -16
            ELSE
               DO 10 I = 2, M
.LT.                  IF( IBLOCK( I )IBLOCK( I-1 ) ) THEN
                     INFO = -6
                     GO TO 20
                  END IF
.EQ..AND..LT.                  IF( IBLOCK( I )IBLOCK( I-1 )  W( I )
     $                W( I-1 ) ) THEN
                     INFO = -5
                     GO TO 20
                  END IF
   10          CONTINUE
   20          CONTINUE
.EQ.               IF( INFO0 ) THEN
.GT.                  IF( ABS( TMPFAC-ORFAC )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
.NE.      IF( INFO0 ) THEN
         CALL PXERBLA( DESCZ( CTXT_ ), 'pdstein', -INFO )
         RETURN
.EQ..OR..EQ.      ELSE IF( LWORK-1  LIWORK-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
*
.EQ..OR..EQ.      IF( N0  M0 )
     $   RETURN
*
.GE.      IF( ORFACZERO ) THEN
         TMPFAC = ORFAC
      ELSE
         TMPFAC = ODM3
      END IF
      ORGFAC = TMPFAC
*
*     Allocate the work among the processes
*
      ILAST = M / LOAD
.EQ.      IF( MOD( M, LOAD )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
.GT.         IF( JNVS ) THEN
            NBLK = IBLOCK( NEXT )
.EQ..AND..NE.            IF( NBLKIBLOCK( NEXT-1 )  NBLKOLNBLK ) THEN
*
*              Compute orthogonalization criterion
*
.EQ.               IF( NBLK1 ) 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
.GT.            IF( TMPFACODM18 ) THEN
               ORTOL = TMPFAC*ONENRM
               DO 80 J = NEXT - 1, MIN( TILL, M-1 )
.NE..OR.                  IF( IBLOCK( J+1 )IBLOCK( J )  W( J+1 )-
.GE.     $                W( J )ORTOL ) THEN
                     GO TO 90
                  END IF
   80          CONTINUE
.EQ..AND..GE.               IF( JM  TILLM )
     $            GO TO 90
               TMPFAC = TMPFAC*ODM1
               GO TO 70
            END IF
   90       CONTINUE
            J = MIN( J, TILL )
         END IF
.EQ.         IF( SELFI )
     $      IM = MAX( 0, J-NVS )
*
         IWORK( I+1 ) = NVS
         NVS = MAX( J, NVS )
  100 CONTINUE
.EQ.      IF( SELFILAST )
     $   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
.NE.         IF( IBLOCK( I )NBLK ) THEN
            NBLK = IBLOCK( I )
.EQ.            IF( NBLK1 ) 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
.GT.         IF( I1 ) THEN
            DIFF = W( I ) - W( I-1 )
.NE..OR..EQ..OR..GT.            IF( IBLOCK( I )IBLOCK( I-1 )  IM  DIFF
     $          ORGFAC*ONENRM ) THEN
               IFIRST = ILAST
.EQ.               IF( IM ) THEN
.NE..OR..GT.                  IF( IBLOCK( M )IBLOCK( M-1 )  DIFFORGFAC*
     $                ONENRM ) THEN
                     ILAST = M - 1
                  ELSE
                     ILAST = M
                  END IF
               ELSE
                  ILAST = I - 1
               END IF
               CLSIZ = ILAST - IFIRST
.GT.               IF( CLSIZ1 ) THEN
.LT.                  IF( LGCLSIZCLSIZ )
     $               LGCLSIZ = CLSIZ
                  MINGAP = ONENRM
  130             CONTINUE
.GT.                  IF( BNDRYP+1 )
     $               GO TO 150
.GT..AND..LT.                  IF( IWORK( BNDRY )IFIRST  IWORK( BNDRY )
     $                ILAST ) THEN
                     MINGAP = MIN( W( IWORK( BNDRY )+1 )-
     $                        W( IWORK( BNDRY ) ), MINGAP )
.GE.                  ELSE IF( IWORK( BNDRY )ILAST ) THEN
.LT.                     IF( MINGAPONENRM ) 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
.EQ.         IF( SELFI-1 ) THEN
            CALL IGEBS2D( DESCZ( CTXT_ ), 'all', ' ', 1, 1, LOCINFO, 1 )
.NE.            IF( LOCINFO0 ) 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 )
.NE.            IF( NERR0 ) 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
.GT.         IF( IWORK( I )IWORK( I+1 ) ) THEN
            ITMP = IWORK( I+1 )
            IWORK( I+1 ) = IWORK( I )
            IWORK( I ) = ITMP
            SORTED = .FALSE.
         END IF
  240 CONTINUE
.NOT.      IF( 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