ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdstein.f
Go to the documentation of this file.
1  SUBROUTINE pdstein( N, D, E, M, W, IBLOCK, ISPLIT, ORFAC, Z, IZ,
2  $ JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, IFAIL,
3  $ ICLUSTR, GAP, INFO )
4 *
5 * -- ScaLAPACK routine (version 1.7) --
6 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7 * and University of California, Berkeley.
8 * November 15, 1997
9 *
10 * .. Scalar Arguments ..
11  INTEGER INFO, IZ, JZ, LIWORK, LWORK, M, N
12  DOUBLE PRECISION ORFAC
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCZ( * ), IBLOCK( * ), ICLUSTR( * ),
16  $ IFAIL( * ), ISPLIT( * ), IWORK( * )
17  DOUBLE PRECISION D( * ), E( * ), GAP( * ), W( * ), WORK( * ),
18  $ Z( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * PDSTEIN computes the eigenvectors of a symmetric tridiagonal matrix
25 * in parallel, using inverse iteration. The eigenvectors found
26 * correspond to user specified eigenvalues. PDSTEIN does not
27 * orthogonalize vectors that are on different processes. The extent
28 * of orthogonalization is controlled by the input parameter LWORK.
29 * Eigenvectors that are to be orthogonalized are computed by the same
30 * process. PDSTEIN decides on the allocation of work among the
31 * processes and then calls DSTEIN2 (modified LAPACK routine) on each
32 * individual process. If insufficient workspace is allocated, the
33 * expected orthogonalization may not be done.
34 *
35 * Note : If the eigenvectors obtained are not orthogonal, increase
36 * LWORK and run the code again.
37 *
38 * Notes
39 * =====
40 *
41 *
42 * Each global data object is described by an associated description
43 * vector. This vector stores the information required to establish
44 * the mapping between an object element and its corresponding process
45 * and memory location.
46 *
47 * Let A be a generic term for any 2D block cyclicly distributed array.
48 * Such a global array has an associated description vector DESCA.
49 * In the following comments, the character _ should be read as
50 * "of the global array".
51 *
52 * NOTATION STORED IN EXPLANATION
53 * --------------- -------------- --------------------------------------
54 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
55 * DTYPE_A = 1.
56 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
57 * the BLACS process grid A is distribu-
58 * ted over. The context itself is glo-
59 * bal, but the handle (the integer
60 * value) may vary.
61 * M_A (global) DESCA( M_ ) The number of rows in the global
62 * array A.
63 * N_A (global) DESCA( N_ ) The number of columns in the global
64 * array A.
65 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
66 * the rows of the array.
67 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
68 * the columns of the array.
69 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
70 * row of the array A is distributed.
71 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
72 * first column of the array A is
73 * distributed.
74 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
75 * array. LLD_A >= MAX(1,LOCr(M_A)).
76 *
77 * Let K be the number of rows or columns of a distributed matrix,
78 * and assume that its process grid has dimension r x c.
79 * LOCr( K ) denotes the number of elements of K that a process
80 * would receive if K were distributed over the r processes of its
81 * process column.
82 * Similarly, LOCc( K ) denotes the number of elements of K that a
83 * process would receive if K were distributed over the c processes of
84 * its process row.
85 * The values of LOCr() and LOCc() may be determined via a call to the
86 * ScaLAPACK tool function, NUMROC:
87 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
88 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
89 * An upper bound for these quantities may be computed by:
90 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
91 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
92 *
93 *
94 * Arguments
95 * =========
96 *
97 * P = NPROW * NPCOL is the total number of processes
98 *
99 * N (global input) INTEGER
100 * The order of the tridiagonal matrix T. N >= 0.
101 *
102 * D (global input) DOUBLE PRECISION array, dimension (N)
103 * The n diagonal elements of the tridiagonal matrix T.
104 *
105 * E (global input) DOUBLE PRECISION array, dimension (N-1)
106 * The (n-1) off-diagonal elements of the tridiagonal matrix T.
107 *
108 * M (global input) INTEGER
109 * The total number of eigenvectors to be found. 0 <= M <= N.
110 *
111 * W (global input/global output) DOUBLE PRECISION array, dim (M)
112 * On input, the first M elements of W contain all the
113 * eigenvalues for which eigenvectors are to be computed. The
114 * eigenvalues should be grouped by split-off block and ordered
115 * from smallest to largest within the block (The output array
116 * W from PDSTEBZ with ORDER='b' is expected here). This
117 * array should be replicated on all processes.
118 * On output, the first M elements contain the input
119 * eigenvalues in ascending order.
120 *
121 * Note : To obtain orthogonal vectors, it is best if
122 * eigenvalues are computed to highest accuracy ( this can be
123 * done by setting ABSTOL to the underflow threshold =
124 * DLAMCH('U') --- ABSTOL is an input parameter
125 * to PDSTEBZ )
126 *
127 * IBLOCK (global input) INTEGER array, dimension (N)
128 * The submatrix indices associated with the corresponding
129 * eigenvalues in W -- 1 for eigenvalues belonging to the
130 * first submatrix from the top, 2 for those belonging to
131 * the second submatrix, etc. (The output array IBLOCK
132 * from PDSTEBZ is expected here).
133 *
134 * ISPLIT (global input) INTEGER array, dimension (N)
135 * The splitting points, at which T breaks up into submatrices.
136 * The first submatrix consists of rows/columns 1 to ISPLIT(1),
137 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
138 * etc., and the NSPLIT-th consists of rows/columns
139 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N (The output array
140 * ISPLIT from PDSTEBZ is expected here.)
141 *
142 * ORFAC (global input) DOUBLE PRECISION
143 * ORFAC specifies which eigenvectors should be orthogonalized.
144 * Eigenvectors that correspond to eigenvalues which are within
145 * ORFAC*||T|| of each other are to be orthogonalized.
146 * However, if the workspace is insufficient (see LWORK), this
147 * tolerance may be decreased until all eigenvectors to be
148 * orthogonalized can be stored in one process.
149 * No orthogonalization will be done if ORFAC equals zero.
150 * A default value of 10^-3 is used if ORFAC is negative.
151 * ORFAC should be identical on all processes.
152 *
153 * Z (local output) DOUBLE PRECISION array,
154 * dimension (DESCZ(DLEN_), N/npcol + NB)
155 * Z contains the computed eigenvectors associated with the
156 * specified eigenvalues. Any vector which fails to converge is
157 * set to its current iterate after MAXITS iterations ( See
158 * DSTEIN2 ).
159 * On output, Z is distributed across the P processes in block
160 * cyclic format.
161 *
162 * IZ (global input) INTEGER
163 * Z's global row index, which points to the beginning of the
164 * submatrix which is to be operated on.
165 *
166 * JZ (global input) INTEGER
167 * Z's global column index, which points to the beginning of
168 * the submatrix which is to be operated on.
169 *
170 * DESCZ (global and local input) INTEGER array of dimension DLEN_.
171 * The array descriptor for the distributed matrix Z.
172 *
173 * WORK (local workspace/global output) DOUBLE PRECISION array,
174 * dimension ( LWORK )
175 * On output, WORK(1) gives a lower bound on the
176 * workspace ( LWORK ) that guarantees the user desired
177 * orthogonalization (see ORFAC).
178 * Note that this may overestimate the minimum workspace needed.
179 *
180 * LWORK (local input) integer
181 * LWORK controls the extent of orthogonalization which can be
182 * done. The number of eigenvectors for which storage is
183 * allocated on each process is
184 * NVEC = floor(( LWORK- max(5*N,NP00*MQ00) )/N).
185 * Eigenvectors corresponding to eigenvalue clusters of size
186 * NVEC - ceil(M/P) + 1 are guaranteed to be orthogonal ( the
187 * orthogonality is similar to that obtained from DSTEIN2).
188 * Note : LWORK must be no smaller than:
189 * max(5*N,NP00*MQ00) + ceil(M/P)*N,
190 * and should have the same input value on all processes.
191 * It is the minimum value of LWORK input on different processes
192 * that is significant.
193 *
194 * If LWORK = -1, then LWORK is global input and a workspace
195 * query is assumed; the routine only calculates the minimum
196 * and optimal size for all work arrays. Each of these
197 * values is returned in the first entry of the corresponding
198 * work array, and no error message is issued by PXERBLA.
199 *
200 * IWORK (local workspace/global output) INTEGER array,
201 * dimension ( 3*N+P+1 )
202 * On return, IWORK(1) contains the amount of integer workspace
203 * required.
204 * On return, the IWORK(2) through IWORK(P+2) indicate
205 * the eigenvectors computed by each process. Process I computes
206 * eigenvectors indexed IWORK(I+2)+1 thru' IWORK(I+3).
207 *
208 * LIWORK (local input) INTEGER
209 * Size of array IWORK. Must be >= 3*N + P + 1
210 *
211 * If LIWORK = -1, then LIWORK is global input and a workspace
212 * query is assumed; the routine only calculates the minimum
213 * and optimal size for all work arrays. Each of these
214 * values is returned in the first entry of the corresponding
215 * work array, and no error message is issued by PXERBLA.
216 *
217 * IFAIL (global output) integer array, dimension (M)
218 * On normal exit, all elements of IFAIL are zero.
219 * If one or more eigenvectors fail to converge after MAXITS
220 * iterations (as in DSTEIN), then INFO > 0 is returned.
221 * If mod(INFO,M+1)>0, then
222 * for I=1 to mod(INFO,M+1), the eigenvector
223 * corresponding to the eigenvalue W(IFAIL(I)) failed to
224 * converge ( W refers to the array of eigenvalues on output ).
225 *
226 * ICLUSTR (global output) integer array, dimension (2*P)
227 * This output array contains indices of eigenvectors
228 * corresponding to a cluster of eigenvalues that could not be
229 * orthogonalized due to insufficient workspace (see LWORK,
230 * ORFAC and INFO). Eigenvectors corresponding to clusters of
231 * eigenvalues indexed ICLUSTR(2*I-1) to ICLUSTR(2*I), I = 1 to
232 * INFO/(M+1), could not be orthogonalized due to lack of
233 * workspace. Hence the eigenvectors corresponding to these
234 * clusters may not be orthogonal. ICLUSTR is a zero terminated
235 * array --- ( ICLUSTR(2*K).NE.0 .AND. ICLUSTR(2*K+1).EQ.0 )
236 * if and only if K is the number of clusters.
237 *
238 * GAP (global output) DOUBLE PRECISION array, dimension (P)
239 * This output array contains the gap between eigenvalues whose
240 * eigenvectors could not be orthogonalized. The INFO/M output
241 * values in this array correspond to the INFO/(M+1) clusters
242 * indicated by the array ICLUSTR. As a result, the dot product
243 * between eigenvectors corresponding to the I^th cluster may be
244 * as high as ( O(n)*macheps ) / GAP(I).
245 *
246 * INFO (global output) INTEGER
247 * = 0: successful exit
248 * < 0: If the i-th argument is an array and the j-entry had
249 * an illegal value, then INFO = -(i*100+j), if the i-th
250 * argument is a scalar and had an illegal value, then
251 * INFO = -i.
252 * < 0 : if INFO = -I, the I-th argument had an illegal value
253 * > 0 : if mod(INFO,M+1) = I, then I eigenvectors failed to
254 * converge in MAXITS iterations. Their indices are
255 * stored in the array IFAIL.
256 * if INFO/(M+1) = I, then eigenvectors corresponding to
257 * I clusters of eigenvalues could not be orthogonalized
258 * due to insufficient workspace. The indices of the
259 * clusters are stored in the array ICLUSTR.
260 *
261 * =====================================================================
262 *
263 * .. Intrinsic Functions ..
264  INTRINSIC abs, dble, max, min, mod
265 * ..
266 * .. External Functions ..
267  INTEGER ICEIL, NUMROC
268  EXTERNAL ICEIL, NUMROC
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL blacs_gridinfo, chk1mat, dgebr2d, dgebs2d,
272  $ dlasrt2, dstein2, igamn2d, igebr2d, igebs2d,
273  $ pchk1mat, pdlaevswp, pxerbla
274 * ..
275 * .. Parameters ..
276  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
277  $ MB_, NB_, RSRC_, CSRC_, LLD_
278  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
279  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
280  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
281  DOUBLE PRECISION ZERO, NEGONE, ODM1, FIVE, ODM3, ODM18
282  PARAMETER ( ZERO = 0.0d+0, negone = -1.0d+0,
283  $ odm1 = 1.0d-1, five = 5.0d+0, odm3 = 1.0d-3,
284  $ odm18 = 1.0d-18 )
285 * ..
286 * .. Local Scalars ..
287  LOGICAL LQUERY, SORTED
288  INTEGER B1, BN, BNDRY, CLSIZ, COL, I, IFIRST, IINFO,
289  $ ilast, im, indrw, itmp, j, k, lgclsiz, llwork,
290  $ load, locinfo, maxvec, mq00, mycol, myrow,
291  $ nblk, nerr, next, np00, npcol, nprow, nvs,
292  $ olnblk, p, row, self, till, toterr
293  DOUBLE PRECISION DIFF, MINGAP, ONENRM, ORGFAC, ORTOL, TMPFAC
294 * ..
295 * .. Local Arrays ..
296  INTEGER IDUM1( 1 ), IDUM2( 1 )
297 * ..
298 * .. Executable Statements ..
299 * This is just to keep ftnchek happy
300  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
301  $ rsrc_.LT.0 )RETURN
302 *
303  CALL blacs_gridinfo( descz( ctxt_ ), nprow, npcol, myrow, mycol )
304  self = myrow*npcol + mycol
305 *
306 * Make sure that we belong to this context (before calling PCHK1MAT)
307 *
308  info = 0
309  IF( nprow.EQ.-1 ) THEN
310  info = -( 1200+ctxt_ )
311  ELSE
312 *
313 * Make sure that NPROW>0 and NPCOL>0 before calling NUMROC
314 *
315  CALL chk1mat( n, 1, n, 1, iz, jz, descz, 12, info )
316  IF( info.EQ.0 ) THEN
317 *
318 * Now we know that our context is good enough to
319 * perform the rest of the checks
320 *
321  np00 = numroc( n, descz( mb_ ), 0, 0, nprow )
322  mq00 = numroc( m, descz( nb_ ), 0, 0, npcol )
323  p = nprow*npcol
324 *
325 * Compute the maximum number of vectors per process
326 *
327  llwork = lwork
328  CALL igamn2d( descz( ctxt_ ), 'A', ' ', 1, 1, llwork, 1, 1,
329  $ 1, -1, -1, -1 )
330  indrw = max( 5*n, np00*mq00 )
331  IF( n.NE.0 )
332  $ maxvec = ( llwork-indrw ) / n
333  load = iceil( m, p )
334  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
335  tmpfac = orfac
336  CALL dgebs2d( descz( ctxt_ ), 'ALL', ' ', 1, 1, tmpfac,
337  $ 1 )
338  ELSE
339  CALL dgebr2d( descz( ctxt_ ), 'ALL', ' ', 1, 1, tmpfac,
340  $ 1, 0, 0 )
341  END IF
342 *
343  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
344  IF( m.LT.0 .OR. m.GT.n ) THEN
345  info = -4
346  ELSE IF( maxvec.LT.load .AND. .NOT.lquery ) THEN
347  info = -14
348  ELSE IF( liwork.LT.3*n+p+1 .AND. .NOT.lquery ) THEN
349  info = -16
350  ELSE
351  DO 10 i = 2, m
352  IF( iblock( i ).LT.iblock( i-1 ) ) THEN
353  info = -6
354  GO TO 20
355  END IF
356  IF( iblock( i ).EQ.iblock( i-1 ) .AND. w( i ).LT.
357  $ w( i-1 ) ) THEN
358  info = -5
359  GO TO 20
360  END IF
361  10 CONTINUE
362  20 CONTINUE
363  IF( info.EQ.0 ) THEN
364  IF( abs( tmpfac-orfac ).GT.five*abs( tmpfac ) )
365  $ info = -8
366  END IF
367  END IF
368 *
369  END IF
370  idum1( 1 ) = m
371  idum2( 1 ) = 4
372  CALL pchk1mat( n, 1, n, 1, iz, jz, descz, 12, 1, idum1, idum2,
373  $ info )
374  work( 1 ) = dble( max( 5*n, np00*mq00 )+iceil( m, p )*n )
375  iwork( 1 ) = 3*n + p + 1
376  END IF
377  IF( info.NE.0 ) THEN
378  CALL pxerbla( descz( ctxt_ ), 'PDSTEIN', -info )
379  RETURN
380  ELSE IF( lwork.EQ.-1 .OR. liwork.EQ.-1 ) THEN
381  RETURN
382  END IF
383 *
384  DO 30 i = 1, m
385  ifail( i ) = 0
386  30 CONTINUE
387  DO 40 i = 1, p + 1
388  iwork( i ) = 0
389  40 CONTINUE
390  DO 50 i = 1, p
391  gap( i ) = negone
392  iclustr( 2*i-1 ) = 0
393  iclustr( 2*i ) = 0
394  50 CONTINUE
395 *
396 *
397 * Quick return if possible
398 *
399  IF( n.EQ.0 .OR. m.EQ.0 )
400  $ RETURN
401 *
402  IF( orfac.GE.zero ) THEN
403  tmpfac = orfac
404  ELSE
405  tmpfac = odm3
406  END IF
407  orgfac = tmpfac
408 *
409 * Allocate the work among the processes
410 *
411  ilast = m / load
412  IF( mod( m, load ).EQ.0 )
413  $ ilast = ilast - 1
414  olnblk = -1
415  nvs = 0
416  next = 1
417  im = 0
418  onenrm = zero
419  DO 100 i = 0, ilast - 1
420  next = next + load
421  j = next - 1
422  IF( j.GT.nvs ) THEN
423  nblk = iblock( next )
424  IF( nblk.EQ.iblock( next-1 ) .AND. nblk.NE.olnblk ) THEN
425 *
426 * Compute orthogonalization criterion
427 *
428  IF( nblk.EQ.1 ) THEN
429  b1 = 1
430  ELSE
431  b1 = isplit( nblk-1 ) + 1
432  END IF
433  bn = isplit( nblk )
434 *
435  onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
436  onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
437  DO 60 j = b1 + 1, bn - 1
438  onenrm = max( onenrm, abs( d( j ) )+abs( e( j-1 ) )+
439  $ abs( e( j ) ) )
440  60 CONTINUE
441  olnblk = nblk
442  END IF
443  till = nvs + maxvec
444  70 CONTINUE
445  j = next - 1
446  IF( tmpfac.GT.odm18 ) THEN
447  ortol = tmpfac*onenrm
448  DO 80 j = next - 1, min( till, m-1 )
449  IF( iblock( j+1 ).NE.iblock( j ) .OR. w( j+1 )-
450  $ w( j ).GE.ortol ) THEN
451  GO TO 90
452  END IF
453  80 CONTINUE
454  IF( j.EQ.m .AND. till.GE.m )
455  $ GO TO 90
456  tmpfac = tmpfac*odm1
457  GO TO 70
458  END IF
459  90 CONTINUE
460  j = min( j, till )
461  END IF
462  IF( self.EQ.i )
463  $ im = max( 0, j-nvs )
464 *
465  iwork( i+1 ) = nvs
466  nvs = max( j, nvs )
467  100 CONTINUE
468  IF( self.EQ.ilast )
469  $ im = m - nvs
470  iwork( ilast+1 ) = nvs
471  DO 110 i = ilast + 2, p + 1
472  iwork( i ) = m
473  110 CONTINUE
474 *
475  clsiz = 1
476  lgclsiz = 1
477  ilast = 0
478  nblk = 0
479  bndry = 2
480  k = 1
481  DO 140 i = 1, m
482  IF( iblock( i ).NE.nblk ) THEN
483  nblk = iblock( i )
484  IF( nblk.EQ.1 ) THEN
485  b1 = 1
486  ELSE
487  b1 = isplit( nblk-1 ) + 1
488  END IF
489  bn = isplit( nblk )
490 *
491  onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
492  onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
493  DO 120 j = b1 + 1, bn - 1
494  onenrm = max( onenrm, abs( d( j ) )+abs( e( j-1 ) )+
495  $ abs( e( j ) ) )
496  120 CONTINUE
497 *
498  END IF
499  IF( i.GT.1 ) THEN
500  diff = w( i ) - w( i-1 )
501  IF( iblock( i ).NE.iblock( i-1 ) .OR. i.EQ.m .OR. diff.GT.
502  $ orgfac*onenrm ) THEN
503  ifirst = ilast
504  IF( i.EQ.m ) THEN
505  IF( iblock( m ).NE.iblock( m-1 ) .OR. diff.GT.orgfac*
506  $ onenrm ) THEN
507  ilast = m - 1
508  ELSE
509  ilast = m
510  END IF
511  ELSE
512  ilast = i - 1
513  END IF
514  clsiz = ilast - ifirst
515  IF( clsiz.GT.1 ) THEN
516  IF( lgclsiz.LT.clsiz )
517  $ lgclsiz = clsiz
518  mingap = onenrm
519  130 CONTINUE
520  IF( bndry.GT.p+1 )
521  $ GO TO 150
522  IF( iwork( bndry ).GT.ifirst .AND. iwork( bndry ).LT.
523  $ ilast ) THEN
524  mingap = min( w( iwork( bndry )+1 )-
525  $ w( iwork( bndry ) ), mingap )
526  ELSE IF( iwork( bndry ).GE.ilast ) THEN
527  IF( mingap.LT.onenrm ) THEN
528  iclustr( 2*k-1 ) = ifirst + 1
529  iclustr( 2*k ) = ilast
530  gap( k ) = mingap / onenrm
531  k = k + 1
532  END IF
533  GO TO 140
534  END IF
535  bndry = bndry + 1
536  GO TO 130
537  END IF
538  END IF
539  END IF
540  140 CONTINUE
541  150 CONTINUE
542  info = ( k-1 )*( m+1 )
543 *
544 * Call DSTEIN2 to find the eigenvectors
545 *
546  CALL dstein2( n, d, e, im, w( iwork( self+1 )+1 ),
547  $ iblock( iwork( self+1 )+1 ), isplit, orgfac,
548  $ work( indrw+1 ), n, work, iwork( p+2 ),
549  $ ifail( iwork( self+1 )+1 ), locinfo )
550 *
551 * Redistribute the eigenvector matrix to conform with the block
552 * cyclic distribution of the input matrix
553 *
554 *
555  DO 160 i = 1, m
556  iwork( p+1+i ) = i
557  160 CONTINUE
558 *
559  CALL dlasrt2( 'I', m, w, iwork( p+2 ), iinfo )
560 *
561  DO 170 i = 1, m
562  iwork( m+p+1+iwork( p+1+i ) ) = i
563  170 CONTINUE
564 *
565 *
566  DO 180 i = 1, locinfo
567  itmp = iwork( self+1 ) + i
568  ifail( itmp ) = ifail( itmp ) + itmp - i
569  ifail( itmp ) = iwork( m+p+1+ifail( itmp ) )
570  180 CONTINUE
571 *
572  DO 190 i = 1, k - 1
573  iclustr( 2*i-1 ) = iwork( m+p+1+iclustr( 2*i-1 ) )
574  iclustr( 2*i ) = iwork( m+p+1+iclustr( 2*i ) )
575  190 CONTINUE
576 *
577 *
578 * Still need to apply the above permutation to IFAIL
579 *
580 *
581  toterr = 0
582  DO 210 i = 1, p
583  IF( self.EQ.i-1 ) THEN
584  CALL igebs2d( descz( ctxt_ ), 'ALL', ' ', 1, 1, locinfo, 1 )
585  IF( locinfo.NE.0 ) THEN
586  CALL igebs2d( descz( ctxt_ ), 'ALL', ' ', locinfo, 1,
587  $ ifail( iwork( i )+1 ), locinfo )
588  DO 200 j = 1, locinfo
589  ifail( toterr+j ) = ifail( iwork( i )+j )
590  200 CONTINUE
591  toterr = toterr + locinfo
592  END IF
593  ELSE
594 *
595  row = ( i-1 ) / npcol
596  col = mod( i-1, npcol )
597 *
598  CALL igebr2d( descz( ctxt_ ), 'ALL', ' ', 1, 1, nerr, 1,
599  $ row, col )
600  IF( nerr.NE.0 ) THEN
601  CALL igebr2d( descz( ctxt_ ), 'ALL', ' ', nerr, 1,
602  $ ifail( toterr+1 ), nerr, row, col )
603  toterr = toterr + nerr
604  END IF
605  END IF
606  210 CONTINUE
607  info = info + toterr
608 *
609 *
610  CALL pdlaevswp( n, work( indrw+1 ), n, z, iz, jz, descz, iwork,
611  $ iwork( m+p+2 ), work, indrw )
612 *
613  DO 220 i = 2, p
614  iwork( i ) = iwork( m+p+1+iwork( i ) )
615  220 CONTINUE
616 *
617 *
618 * Sort the IWORK array
619 *
620 *
621  230 CONTINUE
622  sorted = .true.
623  DO 240 i = 2, p - 1
624  IF( iwork( i ).GT.iwork( i+1 ) ) THEN
625  itmp = iwork( i+1 )
626  iwork( i+1 ) = iwork( i )
627  iwork( i ) = itmp
628  sorted = .false.
629  END IF
630  240 CONTINUE
631  IF( .NOT.sorted )
632  $ GO TO 230
633 *
634  DO 250 i = p + 1, 1, -1
635  iwork( i+1 ) = iwork( i )
636  250 CONTINUE
637 *
638  work( 1 ) = ( lgclsiz+load-1 )*n + indrw
639  iwork( 1 ) = 3*n + p + 1
640 *
641 * End of PDSTEIN
642 *
643  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdstein
subroutine pdstein(N, D, E, M, W, IBLOCK, ISPLIT, ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO)
Definition: pdstein.f:4
dlasrt2
subroutine dlasrt2(ID, N, D, KEY, INFO)
Definition: dlasrt2.f:4
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pdlaevswp
subroutine pdlaevswp(N, ZIN, LDZI, Z, IZ, JZ, DESCZ, NVS, KEY, WORK, LWORK)
Definition: pdlaevswp.f:5
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
dstein2
subroutine dstein2(N, D, E, M, W, IBLOCK, ISPLIT, ORFAC, Z, LDZ, WORK, IWORK, IFAIL, INFO)
Definition: dstein2.f:5
min
#define min(A, B)
Definition: pcgemr.c:181