psstein.f

Go to the documentation of this file.

      SUBROUTINE psstein( 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
      REAL               ORFAC
*     ..
*     .. Array Arguments ..
      INTEGER            DESCZ( * ), IBLOCK( * ), ICLUSTR( * ),
     $                   IFAIL( * ), ISPLIT( * ), IWORK( * )
      REAL               D( * ), E( * ), GAP( * ), W( * ), WORK( * ),
     $                   Z( * )
*     ..
*
*  Purpose
*  =======
*
*  PSSTEIN computes the eigenvectors of a symmetric tridiagonal matrix
*  in parallel, using inverse iteration. The eigenvectors found
*  correspond to user specified eigenvalues. PSSTEIN 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. PSSTEIN decides on the allocation of work among the
*  processes and then calls SSTEIN2 (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) REAL array, dimension (N)
*          The n diagonal elements of the tridiagonal matrix T.
*
*  E       (global input) REAL 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) REAL 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 PSSTEBZ 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 =
*          SLAMCH('U') --- ABSTOL is an input parameter
*          to PSSTEBZ )
*
*  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 PSSTEBZ 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 PSSTEBZ is expected here.)
*
*  ORFAC   (global input) REAL
*          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) REAL 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
*          SSTEIN2 ).
*          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) REAL 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 SSTEIN2).
*          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 SSTEIN), 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) REAL 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, max, min, mod, real
*     ..
*     .. External Functions ..
      INTEGER            ICEIL, NUMROC
      EXTERNAL           ICEIL, NUMROC
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, chk1mat, igamn2d, igebr2d,
     $                   igebs2d, pchk1mat, pslaevswp, pxerbla, sgebr2d,
     $                   sgebs2d, slasrt2, sstein2
*     ..
*     .. 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               ZERO, NEGONE, ODM1, FIVE, ODM3, ODM18
      PARAMETER          ( ZERO = 0.0e+0, negone = -1.0e+0,
     $                   odm1 = 1.0e-1, five = 5.0e+0, odm3 = 1.0e-3,
     $                   odm18 = 1.0e-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
      REAL               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 sgebs2d( descz( ctxt_ ), 'ALL', ' ', 1, 1, tmpfac,
     $                       1 )
            ELSE
               CALL sgebr2d( 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 ) = real( 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_ ), 'PSSTEIN', -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 SSTEIN2 to find the eigenvectors
*
      CALL sstein2( 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 slasrt2( '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 pslaevswp( 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 PSSTEIN
*
      END