OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pcstein.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine pcstein (n, d, e, m, w, iblock, isplit, orfac, z, iz, jz, descz, work, lwork, iwork, liwork, ifail, iclustr, gap, info)

Function/Subroutine Documentation

◆ pcstein()

subroutine pcstein ( integer n,
real, dimension( * ) d,
real, dimension( * ) e,
integer m,
real, dimension( * ) w,
integer, dimension( * ) iblock,
integer, dimension( * ) isplit,
real orfac,
complex, dimension( * ) z,
integer iz,
integer jz,
integer, dimension( * ) descz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer, dimension( * ) ifail,
integer, dimension( * ) iclustr,
real, dimension( * ) gap,
integer info )

Definition at line 1 of file pcstein.f.

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 REAL ORFAC
13* ..
14* .. Array Arguments ..
15 INTEGER DESCZ( * ), IBLOCK( * ), ICLUSTR( * ),
16 $ IFAIL( * ), ISPLIT( * ), IWORK( * )
17 REAL D( * ), E( * ), GAP( * ), W( * ), WORK( * )
18 COMPLEX Z( * )
19* ..
20*
21* Purpose
22* =======
23*
24* PCSTEIN computes the eigenvectors of a symmetric tridiagonal matrix
25* in parallel, using inverse iteration. The eigenvectors found
26* correspond to user specified eigenvalues. PCSTEIN 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. PCSTEIN decides on the allocation of work among the
31* processes and then calls SSTEIN2 (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) REAL array, dimension (N)
103* The n diagonal elements of the tridiagonal matrix T.
104*
105* E (global input) REAL 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) REAL 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 PSSTEBZ 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* SLAMCH('U') --- ABSTOL is an input parameter
125* to PSSTEBZ )
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 PSSTEBZ 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 PSSTEBZ is expected here.)
141*
142* ORFAC (global input) REAL
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) COMPLEX 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* SSTEIN2 ).
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) REAL 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 CSTEIN2).
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 CSTEIN), 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) REAL 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, max, min, mod, real
265* ..
266* .. External Functions ..
267 INTEGER ICEIL, NUMROC
268 EXTERNAL iceil, numroc
269* ..
270* .. External Subroutines ..
271 EXTERNAL blacs_gridinfo, chk1mat, igamn2d, igebr2d,
272 $ igebs2d, pchk1mat, pclaevswp, pxerbla, sgebr2d,
273 $ sgebs2d, slasrt2, sstein2
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 REAL ZERO, NEGONE, ODM1, FIVE, ODM3, ODM18
282 parameter( zero = 0.0e+0, negone = -1.0e+0,
283 $ odm1 = 1.0e-1, five = 5.0e+0, odm3 = 1.0e-3,
284 $ odm18 = 1.0e-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 REAL 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 sgebs2d( descz( ctxt_ ), 'ALL', ' ', 1, 1, TMPFAC,
337 $ 1 )
338 ELSE
339 CALL SGEBR2D( DESCZ( CTXT_ ), 'all', ' ', 1, 1, TMPFAC,
340 $ 1, 0, 0 )
341 END IF
342*
343.EQ..OR..EQ. LQUERY = ( LWORK-1 LIWORK-1 )
344.LT..OR..GT. IF( M0 MN ) THEN
345 INFO = -4
346.LT..AND..NOT. ELSE IF( MAXVECLOAD LQUERY ) THEN
347 INFO = -14
348.LT..AND..NOT. ELSE IF( LIWORK3*N+P+1 LQUERY ) THEN
349 INFO = -16
350 ELSE
351 DO 10 I = 2, M
352.LT. IF( IBLOCK( I )IBLOCK( I-1 ) ) THEN
353 INFO = -6
354 GO TO 20
355 END IF
356.EQ..AND..LT. IF( IBLOCK( I )IBLOCK( I-1 ) W( I )
357 $ W( I-1 ) ) THEN
358 INFO = -5
359 GO TO 20
360 END IF
361 10 CONTINUE
362 20 CONTINUE
363.EQ. IF( INFO0 ) THEN
364.GT. IF( ABS( TMPFAC-ORFAC )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 ) = REAL( MAX( 5*N, NP00*MQ00 )+ICEIL( M, P )*N )
375 IWORK( 1 ) = 3*N + P + 1
376 END IF
377.NE. IF( INFO0 ) THEN
378 CALL PXERBLA( DESCZ( CTXT_ ), 'pcstein', -INFO )
379 RETURN
380.EQ..OR..EQ. ELSE IF( LWORK-1 LIWORK-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.EQ..OR..EQ. IF( N0 M0 )
400 $ RETURN
401*
402.GE. IF( ORFACZERO ) 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.EQ. IF( MOD( M, LOAD )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.GT. IF( JNVS ) THEN
423 NBLK = IBLOCK( NEXT )
424.EQ..AND..NE. IF( NBLKIBLOCK( NEXT-1 ) NBLKOLNBLK ) THEN
425*
426* Compute orthogonalization criterion
427*
428.EQ. IF( NBLK1 ) 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.GT. IF( TMPFACODM18 ) THEN
447 ORTOL = TMPFAC*ONENRM
448 DO 80 J = NEXT - 1, MIN( TILL, M-1 )
449.NE..OR. IF( IBLOCK( J+1 )IBLOCK( J ) W( J+1 )-
450.GE. $ W( J )ORTOL ) THEN
451 GO TO 90
452 END IF
453 80 CONTINUE
454.EQ..AND..GE. IF( JM TILLM )
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.EQ. IF( SELFI )
463 $ IM = MAX( 0, J-NVS )
464*
465 IWORK( I+1 ) = NVS
466 NVS = MAX( J, NVS )
467 100 CONTINUE
468.EQ. IF( SELFILAST )
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.NE. IF( IBLOCK( I )NBLK ) THEN
483 NBLK = IBLOCK( I )
484.EQ. IF( NBLK1 ) 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.GT. IF( I1 ) THEN
500 DIFF = W( I ) - W( I-1 )
501.NE..OR..EQ..OR..GT. IF( IBLOCK( I )IBLOCK( I-1 ) IM DIFF
502 $ ORGFAC*ONENRM ) THEN
503 IFIRST = ILAST
504.EQ. IF( IM ) THEN
505.NE..OR..GT. IF( IBLOCK( M )IBLOCK( M-1 ) DIFFORGFAC*
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.GT. IF( CLSIZ1 ) THEN
516.LT. IF( LGCLSIZCLSIZ )
517 $ LGCLSIZ = CLSIZ
518 MINGAP = ONENRM
519 130 CONTINUE
520.GT. IF( BNDRYP+1 )
521 $ GO TO 150
522.GT..AND..LT. IF( IWORK( BNDRY )IFIRST IWORK( BNDRY )
523 $ ILAST ) THEN
524 MINGAP = MIN( W( IWORK( BNDRY )+1 )-
525 $ W( IWORK( BNDRY ) ), MINGAP )
526.GE. ELSE IF( IWORK( BNDRY )ILAST ) THEN
527.LT. IF( MINGAPONENRM ) 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 SSTEIN2 to find the eigenvectors
545*
546 CALL SSTEIN2( 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 SLASRT2( '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.EQ. IF( SELFI-1 ) THEN
584 CALL IGEBS2D( DESCZ( CTXT_ ), 'all', ' ', 1, 1, LOCINFO, 1 )
585.NE. IF( LOCINFO0 ) 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.NE. IF( NERR0 ) 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 PCLAEVSWP( 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.GT. IF( IWORK( I )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.NOT. IF( 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 PCSTEIN
642*
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sgebs2d
subroutine sgebs2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1072
chk1mat
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition mpi.f:1577
pxerbla
subroutine pxerbla(contxt, srname, info)
Definition mpi.f:1600
sgebr2d
subroutine sgebr2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1113
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition mpi.f:786
pchk1mat
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
Definition pchkxmat.f:3
pclaevswp
subroutine pclaevswp(n, zin, ldzi, z, iz, jz, descz, nvs, key, rwork, lrwork)
Definition pclaevswp.f:5
pcstein
subroutine pcstein(n, d, e, m, w, iblock, isplit, orfac, z, iz, jz, descz, work, lwork, iwork, liwork, ifail, iclustr, gap, info)
Definition pcstein.f:4
slasrt2
subroutine slasrt2(id, n, d, key, info)
Definition slasrt2.f:4
sstein2
subroutine sstein2(n, d, e, m, w, iblock, isplit, orfac, z, ldz, work, iwork, ifail, info)
Definition sstein2.f:5