ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzheevr.f
Go to the documentation of this file.
1  SUBROUTINE pzheevr( JOBZ, RANGE, UPLO, N, A, IA, JA,
2  $ DESCA, VL, VU, IL, IU, M, NZ, W, Z, IZ,
3  $ JZ, DESCZ,
4  $ WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK,
5  $ INFO )
6 
7  IMPLICIT NONE
8 *
9 * -- ScaLAPACK routine (version 2.0.2) --
10 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
11 * May 1 2012
12 *
13 * .. Scalar Arguments ..
14  CHARACTER JOBZ, RANGE, UPLO
15 
16  INTEGER IA, IL, INFO, IU, IZ, JA, JZ, LIWORK, LRWORK,
17  $ LWORK, M, N, NZ
18  DOUBLE PRECISION VL, VU
19 * ..
20 * .. Array Arguments ..
21  INTEGER DESCA( * ), DESCZ( * ), IWORK( * )
22  DOUBLE PRECISION W( * ), RWORK( * )
23  COMPLEX*16 A( * ), WORK( * ), Z( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * PZHEEVR computes selected eigenvalues and, optionally, eigenvectors
30 * of a complex Hermitian matrix A distributed in 2D blockcyclic format
31 * by calling the recommended sequence of ScaLAPACK routines.
32 *
33 * First, the matrix A is reduced to real symmetric tridiagonal form.
34 * Then, the eigenproblem is solved using the parallel MRRR algorithm.
35 * Last, if eigenvectors have been computed, a backtransformation is done.
36 *
37 * Upon successful completion, each processor stores a copy of all computed
38 * eigenvalues in W. The eigenvector matrix Z is stored in
39 * 2D blockcyclic format distributed over all processors.
40 *
41 * For constructive feedback and comments, please contact cvoemel@lbl.gov
42 * C. Voemel
43 *
44 *
45 * Arguments
46 * =========
47 *
48 * JOBZ (global input) CHARACTER*1
49 * Specifies whether or not to compute the eigenvectors:
50 * = 'N': Compute eigenvalues only.
51 * = 'V': Compute eigenvalues and eigenvectors.
52 *
53 * RANGE (global input) CHARACTER*1
54 * = 'A': all eigenvalues will be found.
55 * = 'V': all eigenvalues in the interval [VL,VU] will be found.
56 * = 'I': the IL-th through IU-th eigenvalues will be found.
57 *
58 * UPLO (global input) CHARACTER*1
59 * Specifies whether the upper or lower triangular part of the
60 * symmetric matrix A is stored:
61 * = 'U': Upper triangular
62 * = 'L': Lower triangular
63 *
64 * N (global input) INTEGER
65 * The number of rows and columns of the matrix A. N >= 0
66 *
67 * A (local input/workspace) 2D block cyclic COMPLEX*16 array,
68 * global dimension (N, N),
69 * local dimension ( LLD_A, LOCc(JA+N-1) )
70 * (see Notes below for more detailed explanation of 2d arrays)
71 *
72 * On entry, the symmetric matrix A. If UPLO = 'U', only the
73 * upper triangular part of A is used to define the elements of
74 * the symmetric matrix. If UPLO = 'L', only the lower
75 * triangular part of A is used to define the elements of the
76 * symmetric matrix.
77 *
78 * On exit, the lower triangle (if UPLO='L') or the upper
79 * triangle (if UPLO='U') of A, including the diagonal, is
80 * destroyed.
81 *
82 * IA (global input) INTEGER
83 * A's global row index, which points to the beginning of the
84 * submatrix which is to be operated on.
85 * It should be set to 1 when operating on a full matrix.
86 *
87 * JA (global input) INTEGER
88 * A's global column index, which points to the beginning of
89 * the submatrix which is to be operated on.
90 * It should be set to 1 when operating on a full matrix.
91 *
92 * DESCA (global and local input) INTEGER array of dimension DLEN_.
93 * (The ScaLAPACK descriptor length is DLEN_ = 9.)
94 * The array descriptor for the distributed matrix A.
95 * The descriptor stores details about the 2D block-cyclic
96 * storage, see the notes below.
97 * If DESCA is incorrect, PZHEEVR cannot work correctly.
98 * Also note the array alignment requirements specified below
99 *
100 * VL (global input) DOUBLE PRECISION
101 * If RANGE='V', the lower bound of the interval to be searched
102 * for eigenvalues. Not referenced if RANGE = 'A' or 'I'.
103 *
104 * VU (global input) DOUBLE PRECISION
105 * If RANGE='V', the upper bound of the interval to be searched
106 * for eigenvalues. Not referenced if RANGE = 'A' or 'I'.
107 *
108 * IL (global input) INTEGER
109 * If RANGE='I', the index (from smallest to largest) of the
110 * smallest eigenvalue to be returned. IL >= 1.
111 * Not referenced if RANGE = 'A'.
112 *
113 * IU (global input) INTEGER
114 * If RANGE='I', the index (from smallest to largest) of the
115 * largest eigenvalue to be returned. min(IL,N) <= IU <= N.
116 * Not referenced if RANGE = 'A'.
117 *
118 * M (global output) INTEGER
119 * Total number of eigenvalues found. 0 <= M <= N.
120 *
121 * NZ (global output) INTEGER
122 * Total number of eigenvectors computed. 0 <= NZ <= M.
123 * The number of columns of Z that are filled.
124 * If JOBZ .NE. 'V', NZ is not referenced.
125 * If JOBZ .EQ. 'V', NZ = M
126 *
127 * W (global output) DOUBLE PRECISION array, dimension (N)
128 * On normal exit, the first M entries contain the selected
129 * eigenvalues in ascending order.
130 *
131 * Z (local output) COMPLEX*16 array,
132 * global dimension (N, N),
133 * local dimension ( LLD_Z, LOCc(JZ+N-1) )
134 * If JOBZ = 'V', then on normal exit the first M columns of Z
135 * contain the orthonormal eigenvectors of the matrix
136 * corresponding to the selected eigenvalues.
137 * If JOBZ = 'N', then Z is not referenced.
138 *
139 * IZ (global input) INTEGER
140 * Z's global row index, which points to the beginning of the
141 * submatrix which is to be operated on.
142 * It should be set to 1 when operating on a full matrix.
143 *
144 * JZ (global input) INTEGER
145 * Z's global column index, which points to the beginning of
146 * the submatrix which is to be operated on.
147 * It should be set to 1 when operating on a full matrix.
148 *
149 * DESCZ (global and local input) INTEGER array of dimension DLEN_.
150 * The array descriptor for the distributed matrix Z.
151 * DESCZ( CTXT_ ) must equal DESCA( CTXT_ )
152 *
153 * WORK (local workspace/output) COMPLEX*16 array,
154 * dimension (LWORK)
155 * WORK(1) returns workspace adequate workspace to allow
156 * optimal performance.
157 *
158 * LWORK (local input) INTEGER
159 * Size of WORK array, must be at least 3.
160 * If only eigenvalues are requested:
161 * LWORK >= N + MAX( NB * ( NP00 + 1 ), NB * 3 )
162 * If eigenvectors are requested:
163 * LWORK >= N + ( NP00 + MQ00 + NB ) * NB
164 * For definitions of NP00 & MQ00, see LRWORK.
165 *
166 * For optimal performance, greater workspace is needed, i.e.
167 * LWORK >= MAX( LWORK, NHETRD_LWORK )
168 * Where LWORK is as defined above, and
169 * NHETRD_LWORK = N + 2*( ANB+1 )*( 4*NPS+2 ) +
170 * ( NPS + 1 ) * NPS
171 *
172 * ICTXT = DESCA( CTXT_ )
173 * ANB = PJLAENV( ICTXT, 3, 'PZHETTRD', 'L', 0, 0, 0, 0 )
174 * SQNPC = SQRT( DBLE( NPROW * NPCOL ) )
175 * NPS = MAX( NUMROC( N, 1, 0, 0, SQNPC ), 2*ANB )
176 *
177 * If LWORK = -1, then LWORK is global input and a workspace
178 * query is assumed; the routine only calculates the
179 * optimal size for all work arrays. Each of these
180 * values is returned in the first entry of the corresponding
181 * work array, and no error message is issued by PXERBLA.
182 * NOTE THAT FOR OPTIMAL PERFORMANCE, LWOPT IS RETURNED
183 * (THE OPTIMUM WORKSPACE) RATHER THAN THE MINIMUM NECESSARY
184 * WORKSPACE LWMIN WHEN A WORKSPACE QUERY IS ISSUED.
185 * FOR VERY SMALL MATRICES, LWOPT >> LWMIN.
186 *
187 * RWORK (local workspace/output) DOUBLE PRECISION array,
188 * dimension (LRWORK)
189 * On return, RWORK(1) contains the optimal amount of
190 * workspace required for efficient execution.
191 * if JOBZ='N' RWORK(1) = optimal amount of workspace
192 * required to compute the eigenvalues.
193 * if JOBZ='V' RWORK(1) = optimal amount of workspace
194 * required to compute eigenvalues and eigenvectors.
195 *
196 * LRWORK (local input) INTEGER
197 * Size of RWORK, must be at least 3.
198 * See below for definitions of variables used to define LRWORK.
199 * If no eigenvectors are requested (JOBZ = 'N') then
200 * LRWORK >= 2 + 5 * N + MAX( 12 * N, NB * ( NP00 + 1 ) )
201 * If eigenvectors are requested (JOBZ = 'V' ) then
202 * the amount of workspace required is:
203 * LRWORK >= 2 + 5 * N + MAX( 18*N, NP00 * MQ00 + 2 * NB * NB ) +
204 * (2 + ICEIL( NEIG, NPROW*NPCOL))*N
205 *
206 * Variable definitions:
207 * NEIG = number of eigenvectors requested
208 * NB = DESCA( MB_ ) = DESCA( NB_ ) =
209 * DESCZ( MB_ ) = DESCZ( NB_ )
210 * NN = MAX( N, NB, 2 )
211 * DESCA( RSRC_ ) = DESCA( NB_ ) = DESCZ( RSRC_ ) =
212 * DESCZ( CSRC_ ) = 0
213 * NP00 = NUMROC( NN, NB, 0, 0, NPROW )
214 * MQ00 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL )
215 * ICEIL( X, Y ) is a ScaLAPACK function returning
216 * ceiling(X/Y)
217 *
218 * If LRWORK = -1, then LRWORK is global input and a workspace
219 * query is assumed; the routine only calculates the size
220 * required for optimal performance for all work arrays. Each of
221 * these values is returned in the first entry of the
222 * corresponding work arrays, and no error message is issued by
223 * PXERBLA.
224 *
225 * IWORK (local workspace) INTEGER array
226 * On return, IWORK(1) contains the amount of integer workspace
227 * required.
228 *
229 * LIWORK (local input) INTEGER
230 * size of IWORK
231 *
232 * Let NNP = MAX( N, NPROW*NPCOL + 1, 4 ). Then:
233 * LIWORK >= 12*NNP + 2*N when the eigenvectors are desired
234 * LIWORK >= 10*NNP + 2*N when only the eigenvalues have to be computed
235 *
236 * If LIWORK = -1, then LIWORK is global input and a workspace
237 * query is assumed; the routine only calculates the minimum
238 * and optimal size for all work arrays. Each of these
239 * values is returned in the first entry of the corresponding
240 * work array, and no error message is issued by PXERBLA.
241 *
242 * INFO (global output) INTEGER
243 * = 0: successful exit
244 * < 0: If the i-th argument is an array and the j-entry had
245 * an illegal value, then INFO = -(i*100+j), if the i-th
246 * argument is a scalar and had an illegal value, then
247 * INFO = -i.
248 *
249 * Notes
250 * =====
251 *
252 * Each global data object is described by an associated description
253 * vector. This vector stores the information required to establish
254 * the mapping between an object element and its corresponding process
255 * and memory location.
256 *
257 * Let A be a generic term for any 2D block cyclicly distributed array.
258 * Such a global array has an associated description vector DESCA.
259 * In the following comments, the character _ should be read as
260 * "of the global array".
261 *
262 * NOTATION STORED IN EXPLANATION
263 * --------------- -------------- --------------------------------------
264 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
265 * DTYPE_A = 1.
266 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
267 * the BLACS process grid A is distribu-
268 * ted over. The context itself is glo-
269 * bal, but the handle (the integer
270 * value) may vary.
271 * M_A (global) DESCA( M_ ) The number of rows in the global
272 * array A.
273 * N_A (global) DESCA( N_ ) The number of columns in the global
274 * array A.
275 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
276 * the rows of the array.
277 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
278 * the columns of the array.
279 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
280 * row of the array A is distributed.
281 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
282 * first column of the array A is
283 * distributed.
284 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
285 * array. LLD_A >= MAX(1,LOCr(M_A)).
286 *
287 * Let K be the number of rows or columns of a distributed matrix,
288 * and assume that its process grid has dimension p x q.
289 * LOCr( K ) denotes the number of elements of K that a process
290 * would receive if K were distributed over the p processes of its
291 * process column.
292 * Similarly, LOCc( K ) denotes the number of elements of K that a
293 * process would receive if K were distributed over the q processes of
294 * its process row.
295 * The values of LOCr() and LOCc() may be determined via a call to the
296 * ScaLAPACK tool function, NUMROC:
297 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
298 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
299 * An upper bound for these quantities may be computed by:
300 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
301 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
302 *
303 * PZHEEVR assumes IEEE 754 standard compliant arithmetic.
304 *
305 * Alignment requirements
306 * ======================
307 *
308 * The distributed submatrices A(IA:*, JA:*) and Z(IZ:IZ+M-1,JZ:JZ+N-1)
309 * must satisfy the following alignment properties:
310 *
311 * 1.Identical (quadratic) dimension:
312 * DESCA(M_) = DESCZ(M_) = DESCA(N_) = DESCZ(N_)
313 * 2.Quadratic conformal blocking:
314 * DESCA(MB_) = DESCA(NB_) = DESCZ(MB_) = DESCZ(NB_)
315 * DESCA(RSRC_) = DESCZ(RSRC_)
316 * 3.MOD( IA-1, MB_A ) = MOD( IZ-1, MB_Z ) = 0
317 * 4.IAROW = IZROW
318 *
319 *
320 * .. Parameters ..
321  INTEGER CTXT_, M_, N_,
322  $ MB_, NB_, RSRC_, CSRC_
323  PARAMETER ( CTXT_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
324  $ rsrc_ = 7, csrc_ = 8 )
325  DOUBLE PRECISION ZERO
326  PARAMETER ( ZERO = 0.0d0 )
327 * ..
328 * .. Local Scalars ..
329  LOGICAL ALLEIG, COLBRT, DOBCST, FINISH, FIRST, INDEIG,
330  $ lower, lquery, valeig, vstart, wantz
331  INTEGER ANB, DOL, DOU, DSTCOL, DSTROW, EIGCNT, FRSTCL,
332  $ I, IAROW, ICTXT, IIL, IINDERR, IINDWLC, IINFO,
333  $ iiu, im, indd, indd2, inde, inde2, inderr,
334  $ indilu, indrtau, indrw, indrwork, indtau,
335  $ indwlc, indwork, ipil, ipiu, iproc, izrow,
336  $ lastcl, lengthi, lengthi2, liwmin, llrwork,
337  $ llwork, lrwmin, lrwopt, lwmin, lwopt, maxcls,
338  $ mq00, mycol, myil, myiu, myproc, myrow, mz, nb,
339  $ ndepth, needil, neediu, nhetrd_lwopt, nnp,
340  $ np00, npcol, nprocs, nprow, nps, nsplit,
341  $ offset, parity, rlengthi, rlengthi2, rstarti,
342  $ size1, size2, sqnpc, srccol, srcrow, starti,
343  $ zoffset
344 
345  DOUBLE PRECISION PIVMIN, SAFMIN, SCALE, VLL, VUU, WL,
346  $ WU
347 *
348 * .. Local Arrays ..
349  INTEGER IDUM1( 4 ), IDUM2( 4 )
350 * ..
351 * .. External Functions ..
352  LOGICAL LSAME
353  INTEGER ICEIL, INDXG2P, NUMROC, PJLAENV
354  DOUBLE PRECISION PDLAMCH
355  EXTERNAL ICEIL, INDXG2P, LSAME, NUMROC, PDLAMCH,
356  $ pjlaenv
357 * ..
358 * .. External Subroutines ..
359  EXTERNAL blacs_gridinfo, chk1mat, dcopy, dgebr2d,
360  $ dgebs2d, dgerv2d, dgesd2d, dlarrc, dlasrt2,
361  $ dstegr2a, dstegr2b, dstegr2, igebr2d,
362  $ igebs2d, igerv2d, igesd2d, igsum2d, pchk1mat,
363  $ pchk2mat, pdlared1d, pxerbla, pzelget,
364  $ pzhentrd, pzlaevswp, pzunmtr
365 * ..
366 * .. Intrinsic Functions ..
367  INTRINSIC abs, dble, dcmplx, ichar, int, max, min, mod,
368  $ sqrt
369 * ..
370 * .. Executable Statements ..
371 *
372 
373  info = 0
374 
375 ***********************************************************************
376 *
377 * Decode character arguments to find out what the code should do
378 *
379 ***********************************************************************
380  wantz = lsame( jobz, 'V' )
381  lower = lsame( uplo, 'L' )
382  alleig = lsame( range, 'A' )
383  valeig = lsame( range, 'V' )
384  indeig = lsame( range, 'I' )
385  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
386 
387 ***********************************************************************
388 *
389 * GET MACHINE PARAMETERS
390 *
391 ***********************************************************************
392  ictxt = desca( ctxt_ )
393  safmin = pdlamch( ictxt, 'Safe minimum' )
394 
395 ***********************************************************************
396 *
397 * Set up pointers into the (complex) WORK array
398 *
399 ***********************************************************************
400  indtau = 1
401  indwork = indtau + n
402  llwork = lwork - indwork + 1
403 
404 ***********************************************************************
405 *
406 * Set up pointers into the RWORK array
407 *
408 ***********************************************************************
409  indrtau = 1
410  indd = indrtau + n
411  inde = indd + n + 1
412  indd2 = inde + n + 1
413  inde2 = indd2 + n
414  indrwork = inde2 + n
415  llrwork = lrwork - indrwork + 1
416 
417 ***********************************************************************
418 *
419 * BLACS PROCESSOR GRID SETUP
420 *
421 ***********************************************************************
422  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
423 
424 
425  nprocs = nprow * npcol
426  myproc = myrow * npcol + mycol
427  IF( nprow.EQ.-1 ) THEN
428  info = -( 800+ctxt_ )
429  ELSE IF( wantz ) THEN
430  IF( ictxt.NE.descz( ctxt_ ) ) THEN
431  info = -( 2100+ctxt_ )
432  END IF
433  END IF
434 
435 ***********************************************************************
436 *
437 * COMPUTE REAL WORKSPACE
438 *
439 ***********************************************************************
440  IF ( alleig ) THEN
441  mz = n
442  ELSE IF ( indeig ) THEN
443  mz = iu - il + 1
444  ELSE
445 * Take upper bound for VALEIG case
446  mz = n
447  END IF
448 *
449  nb = desca( nb_ )
450  np00 = numroc( n, nb, 0, 0, nprow )
451  mq00 = numroc( mz, nb, 0, 0, npcol )
452  IF ( wantz ) THEN
453  indrw = indrwork + max(18*n, np00*mq00 + 2*nb*nb)
454  lrwmin = indrw - 1 + (iceil(mz, nprocs) + 2)*n
455  lwmin = n + max((np00 + mq00 + nb) * nb, 3 * nb)
456  ELSE
457  indrw = indrwork + 12*n
458  lrwmin = indrw - 1
459  lwmin = n + max( nb*( np00 + 1 ), 3 * nb )
460  END IF
461 
462 * The code that validates the input requires 3 workspace entries
463  lrwmin = max(3, lrwmin)
464  lrwopt = lrwmin
465  lwmin = max(3, lwmin)
466  lwopt = lwmin
467 *
468  anb = pjlaenv( ictxt, 3, 'PZHETTRD', 'L', 0, 0, 0, 0 )
469  sqnpc = int( sqrt( dble( nprocs ) ) )
470  nps = max( numroc( n, 1, 0, 0, sqnpc ), 2*anb )
471  nhetrd_lwopt = 2*( anb+1 )*( 4*nps+2 ) + ( nps+4 )*nps
472  lwopt = max( lwopt, n+nhetrd_lwopt )
473 *
474  size1 = indrw - indrwork
475 
476 ***********************************************************************
477 *
478 * COMPUTE INTEGER WORKSPACE
479 *
480 ***********************************************************************
481  nnp = max( n, nprocs+1, 4 )
482  IF ( wantz ) THEN
483  liwmin = 12*nnp + 2*n
484  ELSE
485  liwmin = 10*nnp + 2*n
486  END IF
487 
488 ***********************************************************************
489 *
490 * Set up pointers into the IWORK array
491 *
492 ***********************************************************************
493 * Pointer to eigenpair distribution over processors
494  indilu = liwmin - 2*nprocs + 1
495  size2 = indilu - 2*n
496 
497 
498 ***********************************************************************
499 *
500 * Test the input arguments.
501 *
502 ***********************************************************************
503  IF( info.EQ.0 ) THEN
504  CALL chk1mat( n, 4, n, 4, ia, ja, desca, 8, info )
505  IF( wantz )
506  $ CALL chk1mat( n, 4, n, 4, iz, jz, descz, 21, info )
507 *
508  IF( info.EQ.0 ) THEN
509  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
510  info = -1
511  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
512  info = -2
513  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
514  info = -3
515  ELSE IF( mod( ia-1, desca( mb_ ) ).NE.0 ) THEN
516  info = -6
517  ELSE IF( valeig .AND. n.GT.0 .AND. vu.LE.vl ) THEN
518  info = -10
519  ELSE IF( indeig .AND. ( il.LT.1 .OR. il.GT.max( 1, n ) ) )
520  $ THEN
521  info = -11
522  ELSE IF( indeig .AND. ( iu.LT.min( n, il ) .OR. iu.GT.n ))
523  $ THEN
524  info = -12
525  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
526  info = -21
527  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
528  info = -23
529  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
530  info = -25
531  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
532  info = -( 800+nb_ )
533  END IF
534  IF( wantz ) THEN
535  iarow = indxg2p( 1, desca( nb_ ), myrow,
536  $ desca( rsrc_ ), nprow )
537  izrow = indxg2p( 1, desca( nb_ ), myrow,
538  $ descz( rsrc_ ), nprow )
539  IF( iarow.NE.izrow ) THEN
540  info = -19
541  ELSE IF( mod( ia-1, desca( mb_ ) ).NE.
542  $ mod( iz-1, descz( mb_ ) ) ) THEN
543  info = -19
544  ELSE IF( desca( m_ ).NE.descz( m_ ) ) THEN
545  info = -( 2100+m_ )
546  ELSE IF( desca( n_ ).NE.descz( n_ ) ) THEN
547  info = -( 2100+n_ )
548  ELSE IF( desca( mb_ ).NE.descz( mb_ ) ) THEN
549  info = -( 2100+mb_ )
550  ELSE IF( desca( nb_ ).NE.descz( nb_ ) ) THEN
551  info = -( 2100+nb_ )
552  ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) ) THEN
553  info = -( 2100+rsrc_ )
554  ELSE IF( desca( csrc_ ).NE.descz( csrc_ ) ) THEN
555  info = -( 2100+csrc_ )
556  ELSE IF( ictxt.NE.descz( ctxt_ ) ) THEN
557  info = -( 2100+ctxt_ )
558  END IF
559  END IF
560  END IF
561  idum2( 1 ) = 1
562  IF( lower ) THEN
563  idum1( 2 ) = ichar( 'L' )
564  ELSE
565  idum1( 2 ) = ichar( 'U' )
566  END IF
567  idum2( 2 ) = 2
568  IF( alleig ) THEN
569  idum1( 3 ) = ichar( 'A' )
570  ELSE IF( indeig ) THEN
571  idum1( 3 ) = ichar( 'I' )
572  ELSE
573  idum1( 3 ) = ichar( 'V' )
574  END IF
575  idum2( 3 ) = 3
576  IF( lquery ) THEN
577  idum1( 4 ) = -1
578  ELSE
579  idum1( 4 ) = 1
580  END IF
581  idum2( 4 ) = 4
582  IF( wantz ) THEN
583  idum1( 1 ) = ichar( 'V' )
584  CALL pchk2mat( n, 4, n, 4, ia, ja, desca, 8, n, 4, n, 4,iz,
585  $ jz, descz, 21, 4, idum1, idum2, info )
586  ELSE
587  idum1( 1 ) = ichar( 'N' )
588  CALL pchk1mat( n, 4, n, 4, ia, ja, desca, 8, 4, idum1,
589  $ idum2, info )
590  END IF
591  work( 1 ) = dcmplx( lwopt )
592  rwork( 1 ) = dble( lrwopt )
593  iwork( 1 ) = liwmin
594  END IF
595 *
596  IF( info.NE.0 ) THEN
597  CALL pxerbla( ictxt, 'PZHEEVR', -info )
598  RETURN
599  ELSE IF( lquery ) THEN
600  RETURN
601  END IF
602 
603 ***********************************************************************
604 *
605 * Quick return if possible
606 *
607 ***********************************************************************
608  IF( n.EQ.0 ) THEN
609  IF( wantz ) THEN
610  nz = 0
611  END IF
612  m = 0
613  work( 1 ) = dcmplx( lwopt )
614  rwork( 1 ) = dble( lrwopt )
615  iwork( 1 ) = liwmin
616  RETURN
617  END IF
618 
619  IF( valeig ) THEN
620  vll = vl
621  vuu = vu
622  ELSE
623  vll = zero
624  vuu = zero
625  END IF
626 *
627 * No scaling done here, leave this to MRRR kernel.
628 * Scale tridiagonal rather than full matrix.
629 *
630 ***********************************************************************
631 *
632 * REDUCE MATRIX TO REAL SYMMETRIC TRIDIAGONAL FORM.
633 *
634 ***********************************************************************
635 
636 
637  CALL pzhentrd( uplo, n, a, ia, ja, desca, rwork( indd ),
638  $ rwork( inde ), work( indtau ), work( indwork ),
639  $ llwork, rwork( indrwork ), llrwork,iinfo )
640 
641 
642  IF (iinfo .NE. 0) THEN
643  CALL pxerbla( ictxt, 'PZHENTRD', -iinfo )
644  RETURN
645  END IF
646 
647 ***********************************************************************
648 *
649 * DISTRIBUTE TRIDIAGONAL TO ALL PROCESSORS
650 *
651 ***********************************************************************
652  offset = 0
653  IF( ia.EQ.1 .AND. ja.EQ.1 .AND.
654  $ desca( rsrc_ ).EQ.0 .AND. desca( csrc_ ).EQ.0 )
655  $ THEN
656  CALL pdlared1d( n, ia, ja, desca, rwork( indd ),
657  $ rwork( indd2 ), rwork( indrwork ), llrwork )
658 *
659  CALL pdlared1d( n, ia, ja, desca, rwork( inde ),
660  $ rwork( inde2 ), rwork( indrwork ), llrwork )
661  IF( .NOT.lower )
662  $ offset = 1
663  ELSE
664  DO 10 i = 1, n
665  CALL pzelget( 'A', ' ', work( indwork ), a,
666  $ i+ia-1, i+ja-1, desca )
667  rwork( indd2+i-1 ) = dble( work( indwork ) )
668  10 CONTINUE
669  IF( lsame( uplo, 'U' ) ) THEN
670  DO 20 i = 1, n - 1
671  CALL pzelget( 'A', ' ', work( indwork ), a,
672  $ i+ia-1, i+ja, desca )
673  rwork( inde2+i-1 ) = dble( work( indwork ) )
674  20 CONTINUE
675  ELSE
676  DO 30 i = 1, n - 1
677  CALL pzelget( 'A', ' ', work( indwork ), a,
678  $ i+ia, i+ja-1, desca )
679  rwork( inde2+i-1 ) = dble( work( indwork ) )
680  30 CONTINUE
681  END IF
682  END IF
683 
684 
685 
686 
687 ***********************************************************************
688 *
689 * SET IIL, IIU
690 *
691 ***********************************************************************
692  IF ( alleig ) THEN
693  iil = 1
694  iiu = n
695  ELSE IF ( indeig ) THEN
696  iil = il
697  iiu = iu
698  ELSE IF ( valeig ) THEN
699  CALL dlarrc('T', n, vll, vuu, rwork( indd2 ),
700  $ rwork( inde2 + offset ), safmin, eigcnt, iil, iiu, info)
701 * Refine upper bound N that was taken
702  mz = eigcnt
703  iil = iil + 1
704  ENDIF
705 
706  IF(mz.EQ.0) THEN
707  m = 0
708  IF( wantz ) THEN
709  nz = 0
710  END IF
711  work( 1 ) = dble( lwopt )
712  iwork( 1 ) = liwmin
713  RETURN
714  END IF
715 
716  myil = 0
717  myiu = 0
718  m = 0
719  im = 0
720 
721 ***********************************************************************
722 *
723 * COMPUTE WORK ASSIGNMENTS
724 *
725 ***********************************************************************
726 
727 *
728 * Each processor computes the work assignments for all processors
729 *
730  CALL pmpim2( iil, iiu, nprocs,
731  $ iwork(indilu), iwork(indilu+nprocs) )
732 *
733 * Find local work assignment
734 *
735  myil = iwork(indilu+myproc)
736  myiu = iwork(indilu+nprocs+myproc)
737 
738 
739  zoffset = max(0, myil - iil - 1)
740  first = ( myil .EQ. iil )
741 
742 
743 ***********************************************************************
744 *
745 * CALLS TO MRRR KERNEL
746 *
747 ***********************************************************************
748  IF(.NOT.wantz) THEN
749 *
750 * Compute eigenvalues only.
751 *
752  iinfo = 0
753  IF ( myil.GT.0 ) THEN
754  dol = 1
755  dou = myiu - myil + 1
756  CALL dstegr2( jobz, 'I', n, rwork( indd2 ),
757  $ rwork( inde2+offset ), vll, vuu, myil, myiu,
758  $ im, w( 1 ), rwork( indrw ), n,
759  $ myiu - myil + 1,
760  $ iwork( 1 ), rwork( indrwork ), size1,
761  $ iwork( 2*n+1 ), size2,
762  $ dol, dou, zoffset, iinfo )
763 * DSTEGR2 zeroes out the entire W array, so we can't just give
764 * it the part of W we need. So here we copy the W entries into
765 * their correct location
766  DO 49 i = 1, im
767  w( myil-iil+i ) = w( i )
768  49 CONTINUE
769 * W( MYIL ) is at W( MYIL - IIL + 1 )
770 * W( X ) is at W(X - IIL + 1 )
771  END IF
772  IF (iinfo .NE. 0) THEN
773  CALL pxerbla( ictxt, 'DSTEGR2', -iinfo )
774  RETURN
775  END IF
776  ELSEIF ( wantz .AND. nprocs.EQ.1 ) THEN
777 *
778 * Compute eigenvalues and -vectors, but only on one processor
779 *
780  iinfo = 0
781  IF ( myil.GT.0 ) THEN
782  dol = myil - iil + 1
783  dou = myiu - iil + 1
784  CALL dstegr2( jobz, 'I', n, rwork( indd2 ),
785  $ rwork( inde2+offset ), vll, vuu, iil, iiu,
786  $ im, w( 1 ), rwork( indrw ), n,
787  $ n,
788  $ iwork( 1 ), rwork( indrwork ), size1,
789  $ iwork( 2*n+1 ), size2, dol, dou,
790  $ zoffset, iinfo )
791  ENDIF
792  IF (iinfo .NE. 0) THEN
793  CALL pxerbla( ictxt, 'DSTEGR2', -iinfo )
794  RETURN
795  END IF
796  ELSEIF ( wantz ) THEN
797 * Compute representations in parallel.
798 * Share eigenvalue computation for root between all processors
799 * Then compute the eigenvectors.
800  iinfo = 0
801 * Part 1. compute root representations and root eigenvalues
802  IF ( myil.GT.0 ) THEN
803  dol = myil - iil + 1
804  dou = myiu - iil + 1
805  CALL dstegr2a( jobz, 'I', n, rwork( indd2 ),
806  $ rwork( inde2+offset ), vll, vuu, iil, iiu,
807  $ im, w( 1 ), rwork( indrw ), n,
808  $ n, rwork( indrwork ), size1,
809  $ iwork( 2*n+1 ), size2, dol,
810  $ dou, needil, neediu,
811  $ inderr, nsplit, pivmin, scale, wl, wu,
812  $ iinfo )
813  ENDIF
814  IF (iinfo .NE. 0) THEN
815  CALL pxerbla( ictxt, 'DSTEGR2A', -iinfo )
816  RETURN
817  END IF
818 *
819 * The second part of parallel MRRR, the representation tree
820 * construction begins. Upon successful completion, the
821 * eigenvectors have been computed. This is indicated by
822 * the flag FINISH.
823 *
824  vstart = .true.
825  finish = (myil.LE.0)
826 C Part 2. Share eigenvalues and uncertainties between all processors
827  iinderr = indrwork + inderr - 1
828 
829 *
830 
831 
832 *
833 * There are currently two ways to communicate eigenvalue information
834 * using the BLACS.
835 * 1.) BROADCAST
836 * 2.) POINT2POINT between collaborators (those processors working
837 * jointly on a cluster.
838 * For efficiency, BROADCAST has been disabled.
839 * At a later stage, other more efficient communication algorithms
840 * might be implemented, e. g. group or tree-based communication.
841 
842  dobcst = .false.
843  IF(dobcst) THEN
844 * First gather everything on the first processor.
845 * Then use BROADCAST-based communication
846  DO 45 i = 2, nprocs
847  IF (myproc .EQ. (i - 1)) THEN
848  dstrow = 0
849  dstcol = 0
850  starti = dol
851  iwork(1) = starti
852  IF(myil.GT.0) THEN
853  lengthi = myiu - myil + 1
854  ELSE
855  lengthi = 0
856  ENDIF
857  iwork(2) = lengthi
858  CALL igesd2d( ictxt, 2, 1, iwork, 2,
859  $ dstrow, dstcol )
860  IF (( starti.GE.1 ) .AND. ( lengthi.GE.1 )) THEN
861  lengthi2 = 2*lengthi
862 * Copy eigenvalues into communication buffer
863  CALL dcopy(lengthi,w( starti ),1,
864  $ rwork( indd ), 1)
865 * Copy uncertainties into communication buffer
866  CALL dcopy(lengthi,rwork(iinderr+starti-1),1,
867  $ rwork( indd+lengthi ), 1)
868 * send buffer
869  CALL dgesd2d( ictxt, lengthi2,
870  $ 1, rwork( indd ), lengthi2,
871  $ dstrow, dstcol )
872  END IF
873  ELSE IF (myproc .EQ. 0) THEN
874  srcrow = (i-1) / npcol
875  srccol = mod(i-1, npcol)
876  CALL igerv2d( ictxt, 2, 1, iwork, 2,
877  $ srcrow, srccol )
878  starti = iwork(1)
879  lengthi = iwork(2)
880  IF (( starti.GE.1 ) .AND. ( lengthi.GE.1 )) THEN
881  lengthi2 = 2*lengthi
882 * receive buffer
883  CALL dgerv2d( ictxt, lengthi2, 1,
884  $ rwork(indd), lengthi2, srcrow, srccol )
885 * copy eigenvalues from communication buffer
886  CALL dcopy( lengthi, rwork(indd), 1,
887  $ w( starti ), 1)
888 * copy uncertainties (errors) from communication buffer
889  CALL dcopy(lengthi,rwork(indd+lengthi),1,
890  $ rwork( iinderr+starti-1 ), 1)
891  END IF
892  END IF
893  45 CONTINUE
894  lengthi = iiu - iil + 1
895  lengthi2 = lengthi * 2
896  IF (myproc .EQ. 0) THEN
897 * Broadcast eigenvalues and errors to all processors
898  CALL dcopy(lengthi,w ,1, rwork( indd ), 1)
899  CALL dcopy(lengthi,rwork( iinderr ),1,
900  $ rwork( indd+lengthi ), 1)
901  CALL dgebs2d( ictxt, 'A', ' ', lengthi2, 1,
902  $ rwork(indd), lengthi2 )
903  ELSE
904  srcrow = 0
905  srccol = 0
906  CALL dgebr2d( ictxt, 'A', ' ', lengthi2, 1,
907  $ rwork(indd), lengthi2, srcrow, srccol )
908  CALL dcopy( lengthi, rwork(indd), 1, w, 1)
909  CALL dcopy(lengthi,rwork(indd+lengthi),1,
910  $ rwork( iinderr ), 1)
911  END IF
912  ELSE
913 * Enable point2point communication between collaborators
914 
915 * Find collaborators of MYPROC
916  IF( (nprocs.GT.1).AND.(myil.GT.0) ) THEN
917  CALL pmpcol( myproc, nprocs, iil, needil, neediu,
918  $ iwork(indilu), iwork(indilu+nprocs),
919  $ colbrt, frstcl, lastcl )
920  ELSE
921  colbrt = .false.
922  ENDIF
923 
924  IF(colbrt) THEN
925 * If the processor collaborates with others,
926 * communicate information.
927  DO 47 iproc = frstcl, lastcl
928  IF (myproc .EQ. iproc) THEN
929  starti = dol
930  iwork(1) = starti
931  lengthi = myiu - myil + 1
932  iwork(2) = lengthi
933 
934  IF ((starti.GE.1) .AND. (lengthi.GE.1)) THEN
935 * Copy eigenvalues into communication buffer
936  CALL dcopy(lengthi,w( starti ),1,
937  $ rwork(indd), 1)
938 * Copy uncertainties into communication buffer
939  CALL dcopy(lengthi,
940  $ rwork( iinderr+starti-1 ),1,
941  $ rwork(indd+lengthi), 1)
942  ENDIF
943 
944  DO 46 i = frstcl, lastcl
945  IF(i.EQ.myproc) GOTO 46
946  dstrow = i/ npcol
947  dstcol = mod(i, npcol)
948  CALL igesd2d( ictxt, 2, 1, iwork, 2,
949  $ dstrow, dstcol )
950  IF ((starti.GE.1) .AND. (lengthi.GE.1)) THEN
951  lengthi2 = 2*lengthi
952 * send buffer
953  CALL dgesd2d( ictxt, lengthi2,
954  $ 1, rwork(indd), lengthi2,
955  $ dstrow, dstcol )
956  END IF
957  46 CONTINUE
958  ELSE
959  srcrow = iproc / npcol
960  srccol = mod(iproc, npcol)
961  CALL igerv2d( ictxt, 2, 1, iwork, 2,
962  $ srcrow, srccol )
963  rstarti = iwork(1)
964  rlengthi = iwork(2)
965  IF ((rstarti.GE.1 ) .AND. (rlengthi.GE.1 )) THEN
966  rlengthi2 = 2*rlengthi
967  CALL dgerv2d( ictxt, rlengthi2, 1,
968  $ rwork(inde), rlengthi2,
969  $ srcrow, srccol )
970 * copy eigenvalues from communication buffer
971  CALL dcopy( rlengthi,rwork(inde), 1,
972  $ w( rstarti ), 1)
973 * copy uncertainties (errors) from communication buffer
974  CALL dcopy(rlengthi,rwork(inde+rlengthi),1,
975  $ rwork( iinderr+rstarti-1 ), 1)
976  END IF
977  END IF
978  47 CONTINUE
979  ENDIF
980  ENDIF
981 
982 * Part 3. Compute representation tree and eigenvectors.
983 * What follows is a loop in which the tree
984 * is constructed in parallel from top to bottom,
985 * on level at a time, until all eigenvectors
986 * have been computed.
987 *
988  100 CONTINUE
989  IF ( myil.GT.0 ) THEN
990  CALL dstegr2b( jobz, n, rwork( indd2 ),
991  $ rwork( inde2+offset ),
992  $ im, w( 1 ), rwork( indrw ), n, n,
993  $ iwork( 1 ), rwork( indrwork ), size1,
994  $ iwork( 2*n+1 ), size2, dol,
995  $ dou, needil, neediu, indwlc,
996  $ pivmin, scale, wl, wu,
997  $ vstart, finish,
998  $ maxcls, ndepth, parity, zoffset, iinfo )
999  iindwlc = indrwork + indwlc - 1
1000  IF(.NOT.finish) THEN
1001  IF((needil.LT.dol).OR.(neediu.GT.dou)) THEN
1002  CALL pmpcol( myproc, nprocs, iil, needil, neediu,
1003  $ iwork(indilu), iwork(indilu+nprocs),
1004  $ colbrt, frstcl, lastcl )
1005  ELSE
1006  colbrt = .false.
1007  frstcl = myproc
1008  lastcl = myproc
1009  ENDIF
1010 *
1011 * Check if this processor collaborates, i.e.
1012 * communication is needed.
1013 *
1014  IF(colbrt) THEN
1015  DO 147 iproc = frstcl, lastcl
1016  IF (myproc .EQ. iproc) THEN
1017  starti = dol
1018  iwork(1) = starti
1019  IF(myil.GT.0) THEN
1020  lengthi = myiu - myil + 1
1021  ELSE
1022  lengthi = 0
1023  ENDIF
1024  iwork(2) = lengthi
1025  IF ((starti.GE.1).AND.(lengthi.GE.1)) THEN
1026 * Copy eigenvalues into communication buffer
1027  CALL dcopy(lengthi,
1028  $ rwork( iindwlc+starti-1 ),1,
1029  $ rwork(indd), 1)
1030 * Copy uncertainties into communication buffer
1031  CALL dcopy(lengthi,
1032  $ rwork( iinderr+starti-1 ),1,
1033  $ rwork(indd+lengthi), 1)
1034  ENDIF
1035 
1036  DO 146 i = frstcl, lastcl
1037  IF(i.EQ.myproc) GOTO 146
1038  dstrow = i/ npcol
1039  dstcol = mod(i, npcol)
1040  CALL igesd2d( ictxt, 2, 1, iwork, 2,
1041  $ dstrow, dstcol )
1042  IF ((starti.GE.1).AND.(lengthi.GE.1)) THEN
1043  lengthi2 = 2*lengthi
1044 * send buffer
1045  CALL dgesd2d( ictxt, lengthi2,
1046  $ 1, rwork(indd), lengthi2,
1047  $ dstrow, dstcol )
1048  END IF
1049  146 CONTINUE
1050  ELSE
1051  srcrow = iproc / npcol
1052  srccol = mod(iproc, npcol)
1053  CALL igerv2d( ictxt, 2, 1, iwork, 2,
1054  $ srcrow, srccol )
1055  rstarti = iwork(1)
1056  rlengthi = iwork(2)
1057  IF ((rstarti.GE.1).AND.(rlengthi.GE.1)) THEN
1058  rlengthi2 = 2*rlengthi
1059  CALL dgerv2d( ictxt,rlengthi2, 1,
1060  $ rwork(inde),rlengthi2,
1061  $ srcrow, srccol )
1062 * copy eigenvalues from communication buffer
1063  CALL dcopy(rlengthi,rwork(inde), 1,
1064  $ rwork( iindwlc+rstarti-1 ), 1)
1065 * copy uncertainties (errors) from communication buffer
1066  CALL dcopy(rlengthi,rwork(inde+rlengthi),
1067  $ 1,rwork( iinderr+rstarti-1 ), 1)
1068  END IF
1069  END IF
1070  147 CONTINUE
1071  ENDIF
1072  GOTO 100
1073  ENDIF
1074  ENDIF
1075  IF (iinfo .NE. 0) THEN
1076  CALL pxerbla( ictxt, 'DSTEGR2B', -iinfo )
1077  RETURN
1078  END IF
1079 *
1080  ENDIF
1081 
1082 *
1083 ***********************************************************************
1084 *
1085 * MAIN PART ENDS HERE
1086 *
1087 ***********************************************************************
1088 *
1089 
1090 ***********************************************************************
1091 *
1092 * ALLGATHER: EACH PROCESSOR SENDS ITS EIGENVALUES TO THE FIRST ONE,
1093 * THEN THE FIRST PROCESSOR BROADCASTS ALL EIGENVALUES
1094 *
1095 ***********************************************************************
1096 
1097  DO 50 i = 2, nprocs
1098  IF (myproc .EQ. (i - 1)) THEN
1099  dstrow = 0
1100  dstcol = 0
1101  starti = myil - iil + 1
1102  iwork(1) = starti
1103  IF(myil.GT.0) THEN
1104  lengthi = myiu - myil + 1
1105  ELSE
1106  lengthi = 0
1107  ENDIF
1108  iwork(2) = lengthi
1109  CALL igesd2d( ictxt, 2, 1, iwork, 2,
1110  $ dstrow, dstcol )
1111  IF ((starti.GE.1).AND.(lengthi.GE.1)) THEN
1112  CALL dgesd2d( ictxt, lengthi,
1113  $ 1, w( starti ), lengthi,
1114  $ dstrow, dstcol )
1115  ENDIF
1116  ELSE IF (myproc .EQ. 0) THEN
1117  srcrow = (i-1) / npcol
1118  srccol = mod(i-1, npcol)
1119  CALL igerv2d( ictxt, 2, 1, iwork, 2,
1120  $ srcrow, srccol )
1121  starti = iwork(1)
1122  lengthi = iwork(2)
1123  IF ((starti.GE.1).AND.(lengthi.GE.1)) THEN
1124  CALL dgerv2d( ictxt, lengthi, 1,
1125  $ w( starti ), lengthi, srcrow, srccol )
1126  ENDIF
1127  ENDIF
1128  50 CONTINUE
1129 
1130 * Accumulate M from all processors
1131  m = im
1132  CALL igsum2d( ictxt, 'A', ' ', 1, 1, m, 1, -1, -1 )
1133 
1134 * Broadcast eigenvalues to all processors
1135  IF (myproc .EQ. 0) THEN
1136 * Send eigenvalues
1137  CALL dgebs2d( ictxt, 'A', ' ', m, 1, w, m )
1138  ELSE
1139  srcrow = 0
1140  srccol = 0
1141  CALL dgebr2d( ictxt, 'A', ' ', m, 1,
1142  $ w, m, srcrow, srccol )
1143  END IF
1144 *
1145 * Sort the eigenvalues and keep permutation in IWORK to
1146 * sort the eigenvectors accordingly
1147 *
1148  DO 160 i = 1, m
1149  iwork( nprocs+1+i ) = i
1150  160 CONTINUE
1151  CALL dlasrt2( 'I', m, w, iwork( nprocs+2 ), iinfo )
1152  IF (iinfo.NE.0) THEN
1153  CALL pxerbla( ictxt, 'DLASRT2', -iinfo )
1154  RETURN
1155  END IF
1156 
1157 ***********************************************************************
1158 *
1159 * TRANSFORM Z FROM 1D WORKSPACE INTO 2D BLOCKCYCLIC STORAGE
1160 *
1161 ***********************************************************************
1162  IF ( wantz ) THEN
1163  DO 170 i = 1, m
1164  iwork( m+nprocs+1+iwork( nprocs+1+i ) ) = i
1165  170 CONTINUE
1166 * Store NVS in IWORK(1:NPROCS+1) for PZLAEVSWP
1167  iwork( 1 ) = 0
1168  DO 180 i = 1, nprocs
1169 * Find IL and IU for processor i-1
1170 * Has already been computed by PMPIM2 and stored
1171  ipil = iwork(indilu+i-1)
1172  ipiu = iwork(indilu+nprocs+i-1)
1173  IF (ipil .EQ. 0) THEN
1174  iwork( i + 1 ) = iwork( i )
1175  ELSE
1176  iwork( i + 1 ) = iwork( i ) + ipiu - ipil + 1
1177  ENDIF
1178  180 CONTINUE
1179 
1180  IF ( first ) THEN
1181  CALL pzlaevswp(n, rwork( indrw ), n, z, iz, jz,
1182  $ descz, iwork( 1 ), iwork( nprocs+m+2 ), rwork( indrwork ),
1183  $ size1 )
1184  ELSE
1185  CALL pzlaevswp(n, rwork( indrw + n ), n, z, iz, jz,
1186  $ descz, iwork( 1 ), iwork( nprocs+m+2 ), rwork( indrwork ),
1187  $ size1 )
1188  END IF
1189 *
1190  nz = m
1191 *
1192 
1193 ***********************************************************************
1194 *
1195 * Compute eigenvectors of A from eigenvectors of T
1196 *
1197 ***********************************************************************
1198  IF( nz.GT.0 ) THEN
1199  CALL pzunmtr( 'L', uplo, 'N', n, nz, a, ia, ja, desca,
1200  $ work( indtau ), z, iz, jz, descz,
1201  $ work( indwork ), llwork, iinfo )
1202  END IF
1203  IF (iinfo.NE.0) THEN
1204  CALL pxerbla( ictxt, 'PZUNMTR', -iinfo )
1205  RETURN
1206  END IF
1207 *
1208 
1209  END IF
1210 *
1211  work( 1 ) = dcmplx( lwopt )
1212  rwork( 1 ) = dble( lrwopt )
1213  iwork( 1 ) = liwmin
1214 
1215  RETURN
1216 *
1217 * End of PZHEEVR
1218 *
1219  END
pzelget
subroutine pzelget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pzelget.f:2
max
#define max(A, B)
Definition: pcgemr.c:180
pzlaevswp
subroutine pzlaevswp(N, ZIN, LDZI, Z, IZ, JZ, DESCZ, NVS, KEY, RWORK, LRWORK)
Definition: pzlaevswp.f:5
dstegr2b
subroutine dstegr2b(JOBZ, N, D, E, M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK, LIWORK, DOL, DOU, NEEDIL, NEEDIU, INDWLC, PIVMIN, SCALE, WL, WU, VSTART, FINISH, MAXCLS, NDEPTH, PARITY, ZOFFSET, INFO)
Definition: dstegr2b.f:7
pzheevr
subroutine pzheevr(JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, VL, VU, IL, IU, M, NZ, W, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
Definition: pzheevr.f:6
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
pchk2mat
subroutine pchk2mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, MB, MBPOS0, NB, NBPOS0, IB, JB, DESCB, DESCBPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:175
pmpim2
subroutine pmpim2(IL, IU, NPROCS, PMYILS, PMYIUS)
Definition: pmpim2.f:7
pmpcol
subroutine pmpcol(MYPROC, NPROCS, IIL, NEEDIL, NEEDIU, PMYILS, PMYIUS, COLBRT, FRSTCL, LASTCL)
Definition: pmpcol.f:9
pzunmtr
subroutine pzunmtr(SIDE, UPLO, TRANS, M, N, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pzunmtr.f:3
pzhentrd
subroutine pzhentrd(UPLO, N, A, IA, JA, DESCA, D, E, TAU, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pzhentrd.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
dstegr2
subroutine dstegr2(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, WORK, LWORK, IWORK, LIWORK, DOL, DOU, ZOFFSET, INFO)
Definition: dstegr2.f:4
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdlared1d
subroutine pdlared1d(N, IA, JA, DESC, BYCOL, BYALL, WORK, LWORK)
Definition: pdlared1d.f:2
dstegr2a
subroutine dstegr2a(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, WORK, LWORK, IWORK, LIWORK, DOL, DOU, NEEDIL, NEEDIU, INDERR, NSPLIT, PIVMIN, SCALE, WL, WU, INFO)
Definition: dstegr2a.f:6
min
#define min(A, B)
Definition: pcgemr.c:181