Routine: PSSTEIN()  File: SRC\psstein.f

 
 
# lines: 643
  # code: 643
  # comment: 0
  # blank:0
# Variables:77
# Callers:1
# Callings:3
# Words:374
# Keywords:245
 

 

 
..
     .. Array Arguments ..
     ..
  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 ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PSSTEIN( N , D , E , M , W , IBLOCK , ISPLIT , ORFAC , Z , IZ ,
002       $JZ , DESCZ , WORK , LWORK , IWORK , LIWORK , IFAIL ,
003       $ICLUSTR , GAP , INFO )
004  
005  *     -- ScaLAPACK routine(version 1.7) --
006  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
007  *     and University of California , Berkeley.
008  *     November 15 , 1997
009  
010  *     .. Scalar Arguments ..
011        INTEGER INFO , IZ , JZ , LIWORK , LWORK , M , N
012        REAL ORFAC
013        INTRINSIC ABS , MAX , MIN , MOD , REAL
014  *     ..
015  *     .. External Functions ..
016        INTEGER ICEIL , NUMROC
017        EXTERNAL ICEIL , NUMROC
018  *     ..
019  *     .. External Subroutines ..
020        EXTERNAL BLACS_GRIDINFO , CHK1MAT , IGAMN2D , IGEBR2D ,
021       $IGEBS2D , PCHK1MAT , PSLAEVSWP , PXERBLA , SGEBR2D ,
022       $SGEBS2D , SLASRT2 , SSTEIN2  
023  *     ..
024  *     .. Parameters ..
025        INTEGER BLOCK_CYCLIC_2D , DLEN_ , DTYPE_ , CTXT_ , M_ , N_ ,
026       $MB_ , NB_ , RSRC_ , CSRC_ , LLD_
027        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
028       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
029       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
030        REAL ZERO , NEGONE , ODM1 , FIVE , ODM3 , ODM18
031        PARAMETER( ZERO = 0.0E + 0 , NEGONE = - 1.0E + 0 ,
032       $ODM1 = 1.0E - 1 , FIVE = 5.0E + 0 , ODM3 = 1.0E - 3 ,
033       $ODM18 = 1.0E - 18 )
034  *     ..
035  *     .. Local Scalars ..
036        LOGICAL LQUERY , SORTED
037        INTEGER B1 , BN , BNDRY , CLSIZ , COL , I , IFIRST , IINFO ,
038       $ILAST , IM , INDRW , ITMP , J , K , LGCLSIZ , LLWORK ,
039       $LOAD , LOCINFO , MAXVEC , MQ00 , MYCOL , MYROW ,
040       $NBLK , NERR , NEXT , NP00 , NPCOL , NPROW , NVS ,
041       $OLNBLK , P , ROW , SELF , TILL , TOTERR
042        REAL DIFF , MINGAP , ONENRM , ORGFAC , ORTOL , TMPFAC
043  *     ..
044  *     .. Local Arrays ..
045        INTEGER IDUM1( 1 ) , IDUM2( 1 )
046  *     ..
047  *     .. Executable Statements ..
048  *     This is just to keep ftnchek happy
049        IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
050       $    RSRC_.LT.0 )RETURN
051  
052            CALL BLACS_GRIDINFO( DESCZ( CTXT_ ) , NPROW , NPCOL , MYROW , MYCOL )
053            SELF = MYROW*NPCOL + MYCOL
054  
055  *         Make sure that we belong to this context(before calling PCHK1MAT)
056  
057            INFO = 0
058            IF( NPROW.EQ. - 1 ) THEN
059                INFO = - ( 1200 + CTXT_ )
060            ELSE
061  
062  *             Make sure that NPROW > 0 and NPCOL > 0 before calling NUMROC
063  
064                CALL CHK1MAT( N , 1 , N , 1 , IZ , JZ , DESCZ , 12 , INFO )
065                IF( INFO.EQ.0 ) THEN
066  
067  *                 Now we know that our context is good enough to
068  *                 perform the rest of the checks
069  
070                    NP00 = NUMROC( N , DESCZ( MB_ ) , 0 , 0 , NPROW )
071                    MQ00 = NUMROC( M , DESCZ( NB_ ) , 0 , 0 , NPCOL )
072                    P = NPROW*NPCOL
073  
074  *                 Compute the maximum number of vectors per process
075  
076                    LLWORK = LWORK
077                    CALL IGAMN2D( DESCZ( CTXT_ ) , 'A' , ' ' , 1 , 1 , LLWORK , 1 , 1 ,
078       $            1 , - 1 , - 1 , - 1 )
079                    INDRW = MAX( 5*N , NP00*MQ00 )
080                    IF( N.NE.0 )
081       $                MAXVEC =( LLWORK - INDRW ) / N
082                        LOAD = ICEIL( M , P )
083                        IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
084                            TMPFAC = ORFAC
085                            CALL SGEBS2D( DESCZ( CTXT_ ) , 'ALL' , ' ' , 1 , 1 , TMPFAC ,
086       $                    1 )
087                        ELSE
088                            CALL SGEBR2D( DESCZ( CTXT_ ) , 'ALL' , ' ' , 1 , 1 , TMPFAC ,
089       $                    1 , 0 , 0 )
090                        END IF
091  
092                        LQUERY =( LWORK.EQ. - 1 .OR. LIWORK.EQ. - 1 )
093                        IF( M.LT.0 .OR. M.GT.N ) THEN
094                            INFO = - 4
095                        ELSE IF( MAXVEC.LT.LOAD .AND. .NOT.LQUERY ) THEN
096                            INFO = - 14
097                        ELSE IF( LIWORK.LT.3*N + P + 1 .AND. .NOT.LQUERY ) THEN
098                            INFO = - 16
099                        ELSE
100                            DO 10 I = 2 , M
101                                IF( IBLOCK( I ).LT.IBLOCK( I - 1 ) ) THEN
102                                    INFO = - 6
103                                    GO TO 20
104                                END IF
105                                IF( IBLOCK( I ).EQ.IBLOCK( I - 1 ) .AND. W( I ).LT.
106       $                        W( I - 1 ) ) THEN
107                                INFO = - 5
108                                GO TO 20
109                            END IF
110     10                     CONTINUE
111     20 CONTINUE
112        IF( INFO.EQ.0 ) THEN
113            IF( ABS( TMPFAC - ORFAC ).GT.FIVE*ABS( TMPFAC ) )
114       $        INFO = - 8
115            END IF
116        END IF
117  
118        END IF
119        IDUM1( 1 ) = M
120        IDUM2( 1 ) = 4
121        CALL PCHK1MAT( N , 1 , N , 1 , IZ , JZ , DESCZ , 12 , 1 , IDUM1 , IDUM2 ,
122       $INFO )
123        WORK( 1 ) = REAL( MAX( 5*N , NP00*MQ00 ) + ICEIL( M , P )*N )
124        IWORK( 1 ) = 3*N + P + 1
125        END IF
126        IF( INFO.NE.0 ) THEN
127            CALL PXERBLA( DESCZ( CTXT_ ) , 'PSSTEIN' , - INFO )
128            RETURN
129        ELSE IF( LWORK.EQ. - 1 .OR. LIWORK.EQ. - 1 ) THEN
130            RETURN
131        END IF
132  
133        DO 30 I = 1 , M
134            IFAIL( I ) = 0
135     30 CONTINUE
136        DO 40 I = 1 , P + 1
137            IWORK( I ) = 0
138     40 CONTINUE
139        DO 50 I = 1 , P
140            GAP( I ) = NEGONE
141            ICLUSTR( 2*I - 1 ) = 0
142            ICLUSTR( 2*I ) = 0
143     50 CONTINUE
144  
145  *     Quick return if possible
146  
147        IF( N.EQ.0 .OR. M.EQ.0 )
148       $    RETURN
149  
150            IF( ORFAC.GE.ZERO ) THEN
151                TMPFAC = ORFAC
152            ELSE
153                TMPFAC = ODM3
154            END IF
155            ORGFAC = TMPFAC
156  
157  *         Allocate the work among the processes
158  
159            ILAST = M / LOAD
160            IF( MOD( M , LOAD ).EQ.0 )
161       $        ILAST = ILAST - 1
162                OLNBLK = - 1
163                NVS = 0
164            NEXT = 1
165            IM = 0
166            ONENRM = ZERO
167            DO 100 I = 0 , ILAST - 1
168            NEXT = NEXT + LOAD
169            J = NEXT - 1
170            IF( J.GT.NVS ) THEN
171                NBLK = IBLOCK( NEXT )
172                IF( NBLK.EQ.IBLOCK( NEXT - 1 ) .AND. NBLK.NE.OLNBLK ) THEN
173  
174  *                 Compute orthogonalization criterion
175  
176                    IF( NBLK.EQ.1 ) THEN
177                        B1 = 1
178                    ELSE
179                        B1 = ISPLIT( NBLK - 1 ) + 1
180                    END IF
181                    BN = ISPLIT( NBLK )
182  
183                    ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
184                    ONENRM = MAX( ONENRM , ABS( D( BN ) ) + ABS( E( BN - 1 ) ) )
185                    DO 60 J = B1 + 1 , BN - 1
186                        ONENRM = MAX( ONENRM , ABS( D( J ) ) + ABS( E( J - 1 ) ) +
187       $                ABS( E( J ) ) )
188     60             CONTINUE
189                    OLNBLK = NBLK
190                END IF
191                TILL = NVS + MAXVEC
192     70 CONTINUE
193        J = NEXT - 1
194        IF( TMPFAC.GT.ODM18 ) THEN
195            ORTOL = TMPFAC*ONENRM
196            DO 80 J = NEXT - 1 , MIN( TILL , M - 1 )
197                IF( IBLOCK( J + 1 ).NE.IBLOCK( J ) .OR. W( J + 1 ) -
198       $        W( J ).GE.ORTOL ) THEN
199                GO TO 90
200            END IF
201     80     CONTINUE
202            IF( J.EQ.M .AND. TILL.GE.M )
203       $        GO TO 90
204                TMPFAC = TMPFAC*ODM1
205                GO TO 70
206            END IF
207     90 CONTINUE
208        J = MIN( J , TILL )
209        END IF
210        IF( SELF.EQ.I )
211       $    IM = MAX( 0 , J - NVS )
212  
213            IWORK( I + 1 ) = NVS
214            NVS = MAX( J , NVS )
215    100 CONTINUE
216        IF( SELF.EQ.ILAST )
217       $    IM = M - NVS
218            IWORK( ILAST + 1 ) = NVS
219            DO 110 I = ILAST + 2 , P + 1
220                IWORK( I ) = M
221    110     CONTINUE
222  
223            CLSIZ = 1
224            LGCLSIZ = 1
225            ILAST = 0
226            NBLK = 0
227            BNDRY = 2
228            K = 1
229            DO 140 I = 1 , M
230                IF( IBLOCK( I ).NE.NBLK ) THEN
231                    NBLK = IBLOCK( I )
232                    IF( NBLK.EQ.1 ) THEN
233                        B1 = 1
234                    ELSE
235                        B1 = ISPLIT( NBLK - 1 ) + 1
236                    END IF
237                    BN = ISPLIT( NBLK )
238  
239                    ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
240                    ONENRM = MAX( ONENRM , ABS( D( BN ) ) + ABS( E( BN - 1 ) ) )
241                    DO 120 J = B1 + 1 , BN - 1
242                        ONENRM = MAX( ONENRM , ABS( D( J ) ) + ABS( E( J - 1 ) ) +
243       $                ABS( E( J ) ) )
244    120             CONTINUE
245  
246                END IF
247                IF( I.GT.1 ) THEN
248                    DIFF = W( I ) - W( I - 1 )
249                    IF( IBLOCK( I ).NE.IBLOCK( I - 1 ) .OR. I.EQ.M .OR. DIFF.GT.
250       $            ORGFAC*ONENRM ) THEN
251                    IFIRST = ILAST
252                    IF( I.EQ.M ) THEN
253                        IF( IBLOCK( M ).NE.IBLOCK( M - 1 ) .OR. DIFF.GT.ORGFAC*
254       $                ONENRM ) THEN
255                        ILAST = M - 1
256                    ELSE
257                        ILAST = M
258                    END IF
259                ELSE
260                    ILAST = I - 1
261                END IF
262                CLSIZ = ILAST - IFIRST
263                IF( CLSIZ.GT.1 ) THEN
264                    IF( LGCLSIZ.LT.CLSIZ )
265       $                LGCLSIZ = CLSIZ
266                        MINGAP = ONENRM
267    130 CONTINUE
268        IF( BNDRY.GT.P + 1 )
269       $    GO TO 150
270            IF( IWORK( BNDRY ).GT.IFIRST .AND. IWORK( BNDRY ).LT.
271       $    ILAST ) THEN
272            MINGAP = MIN( W( IWORK( BNDRY ) + 1 ) -
273       $    W( IWORK( BNDRY ) ) , MINGAP )
274        ELSE IF( IWORK( BNDRY ).GE.ILAST ) THEN
275            IF( MINGAP.LT.ONENRM ) THEN
276                ICLUSTR( 2*K - 1 ) = IFIRST + 1
277                ICLUSTR( 2*K ) = ILAST
278                GAP( K ) = MINGAP / ONENRM
279                K = K + 1
280            END IF
281            GO TO 140
282        END IF
283        BNDRY = BNDRY + 1
284        GO TO 130
285        END IF
286        END IF
287        END IF
288    140 CONTINUE
289    150 CONTINUE
290        INFO =( K - 1 )*( M + 1 )
291  
292  *     Call SSTEIN2 to find the eigenvectors
293  
294        CALL SSTEIN2 ( N , D , E , IM , W( IWORK( SELF + 1 ) + 1 ) ,
295       $IBLOCK( IWORK( SELF + 1 ) + 1 ) , ISPLIT , ORGFAC ,
296       $WORK( INDRW + 1 ) , N , WORK , IWORK( P + 2 ) ,
297       $IFAIL( IWORK( SELF + 1 ) + 1 ) , LOCINFO )
298  
299  *     Redistribute the eigenvector matrix to conform with the block
300  *     cyclic distribution of the input matrix
301  
302        DO 160 I = 1 , M
303            IWORK( P + 1 + I ) = I
304    160 CONTINUE
305  
306        CALL SLASRT2 ( 'I' , M , W , IWORK( P + 2 ) , IINFO )
307  
308        DO 170 I = 1 , M
309            IWORK( M + P + 1 + IWORK( P + 1 + I ) ) = I
310    170 CONTINUE
311  
312        DO 180 I = 1 , LOCINFO
313            ITMP = IWORK( SELF + 1 ) + I
314            IFAIL( ITMP ) = IFAIL( ITMP ) + ITMP - I
315            IFAIL( ITMP ) = IWORK( M + P + 1 + IFAIL( ITMP ) )
316    180 CONTINUE
317  
318        DO 190 I = 1 , K - 1
319            ICLUSTR( 2*I - 1 ) = IWORK( M + P + 1 + ICLUSTR( 2*I - 1 ) )
320            ICLUSTR( 2*I ) = IWORK( M + P + 1 + ICLUSTR( 2*I ) )
321    190 CONTINUE
322  
323  *     Still need to apply the above permutation to IFAIL
324  
325        TOTERR = 0
326        DO 210 I = 1 , P
327            IF( SELF.EQ.I - 1 ) THEN
328                CALL IGEBS2D( DESCZ( CTXT_ ) , 'ALL' , ' ' , 1 , 1 , LOCINFO , 1 )
329                IF( LOCINFO.NE.0 ) THEN
330                    CALL IGEBS2D( DESCZ( CTXT_ ) , 'ALL' , ' ' , LOCINFO , 1 ,
331       $            IFAIL( IWORK( I ) + 1 ) , LOCINFO )
332                    DO 200 J = 1 , LOCINFO
333                        IFAIL( TOTERR + J ) = IFAIL( IWORK( I ) + J )
334    200             CONTINUE
335                    TOTERR = TOTERR + LOCINFO
336                END IF
337            ELSE
338  
339                ROW =( I - 1 ) / NPCOL
340                COL = MOD( I - 1 , NPCOL )
341  
342                CALL IGEBR2D( DESCZ( CTXT_ ) , 'ALL' , ' ' , 1 , 1 , NERR , 1 ,
343       $        ROW , COL )
344                IF( NERR.NE.0 ) THEN
345                    CALL IGEBR2D( DESCZ( CTXT_ ) , 'ALL' , ' ' , NERR , 1 ,
346       $            IFAIL( TOTERR + 1 ) , NERR , ROW , COL )
347                    TOTERR = TOTERR + NERR
348                END IF
349            END IF
350    210 CONTINUE
351        INFO = INFO + TOTERR
352  
353        CALL PSLAEVSWP ( N , WORK( INDRW + 1 ) , N , Z , IZ , JZ , DESCZ , IWORK ,
354       $IWORK( M + P + 2 ) , WORK , INDRW )
355  
356        DO 220 I = 2 , P
357            IWORK( I ) = IWORK( M + P + 1 + IWORK( I ) )
358    220 CONTINUE
359  
360  *     Sort the IWORK array
361  
362    230 CONTINUE
363        SORTED = .TRUE.
364        DO 240 I = 2 , P - 1
365            IF( IWORK( I ).GT.IWORK( I + 1 ) ) THEN
366                ITMP = IWORK( I + 1 )
367                IWORK( I + 1 ) = IWORK( I )
368                IWORK( I ) = ITMP
369                SORTED = .FALSE.
370            END IF
371    240 CONTINUE
372        IF( .NOT.SORTED )
373       $    GO TO 230
374  
375            DO 250 I = P + 1 , 1 , - 1
376                IWORK( I + 1 ) = IWORK( I )
377    250     CONTINUE
378  
379            WORK( 1 ) =( LGCLSIZ + LOAD - 1 )*N + INDRW
380            IWORK( 1 ) = 3*N + P + 1
381  
382  *         End of PSSTEIN
383  
384        END