ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzheev.f
Go to the documentation of this file.
1  SUBROUTINE pzheev( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ,
2  $ DESCZ, WORK, LWORK, RWORK, LRWORK, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * August 14, 2001
8 *
9 * .. Scalar Arguments ..
10  CHARACTER JOBZ, UPLO
11  INTEGER IA, INFO, IZ, JA, JZ, LRWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * ), DESCZ( * )
15  DOUBLE PRECISION RWORK( * ), W( * )
16  COMPLEX*16 A( * ), WORK( * ), Z( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PZHEEV computes selected eigenvalues and, optionally, eigenvectors
23 * of a complex Hermitian matrix A by calling the recommended sequence
24 * of ScaLAPACK routines.
25 *
26 * In its present form, PZHEEV assumes a homogeneous system and makes
27 * only spot checks of the consistency of the eigenvalues across the
28 * different processes. Because of this, it is possible that a
29 * heterogeneous system may return incorrect results without any error
30 * messages.
31 *
32 * Notes
33 * =====
34 * A description vector is associated with each 2D block-cyclicly dis-
35 * tributed matrix. This vector stores the information required to
36 * establish the mapping between a matrix entry and its corresponding
37 * process and memory location.
38 *
39 * In the following comments, the character _ should be read as
40 * "of the distributed matrix". Let A be a generic term for any 2D
41 * block cyclicly distributed matrix. Its description vector is DESCA:
42 *
43 * NOTATION STORED IN EXPLANATION
44 * --------------- -------------- --------------------------------------
45 * DTYPE_A(global) DESCA( DTYPE_) The descriptor type.
46 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
47 * the BLACS process grid A is distribu-
48 * ted over. The context itself is glo-
49 * bal, but the handle (the integer
50 * value) may vary.
51 * M_A (global) DESCA( M_ ) The number of rows in the distributed
52 * matrix A.
53 * N_A (global) DESCA( N_ ) The number of columns in the distri-
54 * buted matrix A.
55 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
56 * the rows of A.
57 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
58 * the columns of A.
59 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
60 * row of the matrix A is distributed.
61 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
62 * first column of A is distributed.
63 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
64 * array storing the local blocks of the
65 * distributed matrix A.
66 * LLD_A >= MAX(1,LOCr(M_A)).
67 *
68 * Let K be the number of rows or columns of a distributed matrix,
69 * and assume that its process grid has dimension p x q.
70 * LOCr( K ) denotes the number of elements of K that a process
71 * would receive if K were distributed over the p processes of its
72 * process column.
73 * Similarly, LOCc( K ) denotes the number of elements of K that a
74 * process would receive if K were distributed over the q processes of
75 * its process row.
76 * The values of LOCr() and LOCc() may be determined via a call to the
77 * ScaLAPACK tool function, NUMROC:
78 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
79 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
80 *
81 * Arguments
82 * =========
83 *
84 * NP = the number of rows local to a given process.
85 * NQ = the number of columns local to a given process.
86 *
87 * JOBZ (global input) CHARACTER*1
88 * Specifies whether or not to compute the eigenvectors:
89 * = 'N': Compute eigenvalues only.
90 * = 'V': Compute eigenvalues and eigenvectors.
91 *
92 * UPLO (global input) CHARACTER*1
93 * Specifies whether the upper or lower triangular part of the
94 * Hermitian matrix A is stored:
95 * = 'U': Upper triangular
96 * = 'L': Lower triangular
97 *
98 * N (global input) INTEGER
99 * The number of rows and columns of the matrix A. N >= 0.
100 *
101 * A (local input/workspace) block cyclic COMPLEX*16 array,
102 * global dimension (N, N), local dimension ( LLD_A,
103 * LOCc(JA+N-1) )
104 *
105 * On entry, the Hermitian matrix A. If UPLO = 'U', only the
106 * upper triangular part of A is used to define the elements of
107 * the Hermitian matrix. If UPLO = 'L', only the lower
108 * triangular part of A is used to define the elements of the
109 * Hermitian matrix.
110 *
111 * On exit, the lower triangle (if UPLO='L') or the upper
112 * triangle (if UPLO='U') of A, including the diagonal, is
113 * destroyed.
114 *
115 * IA (global input) INTEGER
116 * A's global row index, which points to the beginning of the
117 * submatrix which is to be operated on.
118 *
119 * JA (global input) INTEGER
120 * A's global column index, which points to the beginning of
121 * the submatrix which is to be operated on.
122 *
123 * DESCA (global and local input) INTEGER array of dimension DLEN_.
124 * The array descriptor for the distributed matrix A.
125 * If DESCA( CTXT_ ) is incorrect, PZHEEV cannot guarantee
126 * correct error reporting.
127 *
128 * W (global output) DOUBLE PRECISION array, dimension (N)
129 * If INFO=0, the 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 *
143 * JZ (global input) INTEGER
144 * Z's global column index, which points to the beginning of
145 * the submatrix which is to be operated on.
146 *
147 * DESCZ (global and local input) INTEGER array of dimension DLEN_.
148 * The array descriptor for the distributed matrix Z.
149 * DESCZ( CTXT_ ) must equal DESCA( CTXT_ )
150 *
151 * WORK (local workspace/output) COMPLEX*16 array,
152 * dimension (LWORK)
153 * On output, WORK(1) returns the workspace needed to guarantee
154 * completion. If the input parameters are incorrect, WORK(1)
155 * may also be incorrect.
156 *
157 * If JOBZ='N' WORK(1) = minimal workspace for eigenvalues only.
158 * If JOBZ='V' WORK(1) = minimal workspace required to
159 * generate all the eigenvectors.
160 *
161 *
162 * LWORK (local input) INTEGER
163 * See below for definitions of variables used to define LWORK.
164 * If no eigenvectors are requested (JOBZ = 'N') then
165 * LWORK >= MAX( NB*( NP0+1 ), 3 ) +3*N
166 * If eigenvectors are requested (JOBZ = 'V' ) then
167 * the amount of workspace required:
168 * LWORK >= (NP0 + NQ0 + NB)*NB + 3*N + N^2
169 *
170 * Variable definitions:
171 * NB = DESCA( MB_ ) = DESCA( NB_ ) =
172 * DESCZ( MB_ ) = DESCZ( NB_ )
173 * NP0 = NUMROC( NN, NB, 0, 0, NPROW )
174 * NQ0 = NUMROC( MAX( N, NB, 2 ), NB, 0, 0, NPCOL )
175 *
176 * If LWORK = -1, the LWORK is global input and a workspace
177 * query is assumed; the routine only calculates the minimum
178 * size for the WORK array. The required workspace is returned
179 * as the first element of WORK and no error message is issued
180 * by PXERBLA.
181 *
182 * RWORK (local workspace/output) COMPLEX*16 array,
183 * dimension (LRWORK)
184 * On output RWORK(1) returns the
185 * DOUBLE PRECISION workspace needed to
186 * guarantee completion. If the input parameters are incorrect,
187 * RWORK(1) may also be incorrect.
188 *
189 * LRWORK (local input) INTEGER
190 * Size of RWORK array.
191 * If eigenvectors are desired (JOBZ = 'V') then
192 * LRWORK >= 2*N + 2*N-2
193 * If eigenvectors are not desired (JOBZ = 'N') then
194 * LRWORK >= 2*N
195 *
196 * If LRWORK = -1, the LRWORK is global input and a workspace
197 * query is assumed; the routine only calculates the minimum
198 * size for the RWORK array. The required workspace is returned
199 * as the first element of RWORK and no error message is issued
200 * by PXERBLA.
201 *
202 * INFO (global output) INTEGER
203 * = 0: successful exit
204 * < 0: If the i-th argument is an array and the j-entry had
205 * an illegal value, then INFO = -(i*100+j), if the i-th
206 * argument is a scalar and had an illegal value, then
207 * INFO = -i.
208 * > 0: If INFO = 1 through N, the i(th) eigenvalue did not
209 * converge in ZSTEQR2 after a total of 30*N iterations.
210 * If INFO = N+1, then PZHEEV has detected heterogeneity
211 * by finding that eigenvalues were not identical across
212 * the process grid. In this case, the accuracy of
213 * the results from PZHEEV cannot be guaranteed.
214 *
215 * Alignment requirements
216 * ======================
217 *
218 * The distributed submatrices A(IA:*, JA:*) and C(IC:IC+M-1,JC:JC+N-1)
219 * must verify some alignment properties, namely the following
220 * expressions should be true:
221 *
222 * ( MB_A.EQ.NB_A.EQ.MB_Z .AND. IROFFA.EQ.IROFFZ .AND. IROFFA.EQ.0 .AND.
223 * IAROW.EQ.IZROW )
224 * where
225 * IROFFA = MOD( IA-1, MB_A ) and ICOFFA = MOD( JA-1, NB_A ).
226 *
227 * =====================================================================
228 *
229 * Version 1.4 limitations:
230 * DESCA(MB_) = DESCA(NB_)
231 * DESCA(M_) = DESCZ(M_)
232 * DESCA(N_) = DESCZ(N_)
233 * DESCA(MB_) = DESCZ(MB_)
234 * DESCA(NB_) = DESCZ(NB_)
235 * DESCA(RSRC_) = DESCZ(RSRC_)
236 *
237 * .. Parameters ..
238  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
239  $ mb_, nb_, rsrc_, csrc_, lld_
240  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
241  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
242  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
243  DOUBLE PRECISION ZERO, ONE
244  parameter( zero = 0.0d+0, one = 1.0d+0 )
245  COMPLEX*16 CZERO, CONE
246  parameter( czero = ( 0.0d+0, 0.0d+0 ),
247  $ cone = ( 1.0d+0, 0.0d+0 ) )
248  INTEGER ITHVAL
249  parameter( ithval = 10 )
250 * ..
251 * .. Local Scalars ..
252  LOGICAL LOWER, WANTZ
253  INTEGER CONTEXTC, CSRC_A, I, IACOL, IAROW, ICOFFA,
254  $ iinfo, indd, inde, indrd, indre, indrwork,
255  $ indtau, indwork, indwork2, iroffa, iroffz,
256  $ iscale, izrow, j, k, ldc, llrwork, llwork,
257  $ lrmin, lrwmin, lwmin, mb_a, mb_z, mycol,
258  $ mypcolc, myprowc, myrow, nb, nb_a, nb_z, np0,
259  $ npcol, npcolc, nprocs, nprow, nprowc, nq0, nrc,
260  $ rsizezsteqr2, rsrc_a, rsrc_z, sizepzhetrd,
261  $ sizepzunmtr, sizezsteqr2
262  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
263  $ smlnum
264 * ..
265 * .. Local Arrays ..
266  INTEGER DESCQR( 10 ), IDUM1( 3 ), IDUM2( 3 )
267 * ..
268 * .. External Functions ..
269  LOGICAL LSAME
270  INTEGER INDXG2P, NUMROC, SL_GRIDRESHAPE
271  DOUBLE PRECISION PDLAMCH, PZLANHE
272  EXTERNAL lsame, indxg2p, numroc, sl_gridreshape,
273  $ pdlamch, pzlanhe
274 * ..
275 * .. External Subroutines ..
276  EXTERNAL blacs_gridexit, blacs_gridinfo, chk1mat, dcopy,
277  $ descinit, dgamn2d, dgamx2d, dscal, pchk1mat,
278  $ pchk2mat, pxerbla, pzelget, pzgemr2d, pzhetrd,
279  $ pzlascl, pzlaset, pzunmtr, zsteqr2
280 * ..
281 * .. Intrinsic Functions ..
282  INTRINSIC abs, dble, dcmplx, ichar, int, max, min, mod,
283  $ sqrt
284 * ..
285 * .. Executable Statements ..
286 * This is just to keep ftnchek and toolpack/1 happy
287  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
288  $ rsrc_.LT.0 )RETURN
289 *
290 * Quick return
291 *
292  IF( n.EQ.0 )
293  $ RETURN
294 *
295 * Test the input arguments.
296 *
297  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
298  info = 0
299 *
300 * Initialize pointer to some safe value
301 *
302  indtau = 1
303  indd = 1
304  inde = 1
305  indwork = 1
306  indwork2 = 1
307 *
308  indre = 1
309  indrd = 1
310  indrwork = 1
311 *
312  wantz = lsame( jobz, 'V' )
313  IF( nprow.EQ.-1 ) THEN
314  info = -( 700+ctxt_ )
315  ELSE IF( wantz ) THEN
316  IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
317  info = -( 1200+ctxt_ )
318  END IF
319  END IF
320  IF( info.EQ.0 ) THEN
321  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
322  IF( wantz )
323  $ CALL chk1mat( n, 3, n, 3, iz, jz, descz, 12, info )
324 *
325  IF( info.EQ.0 ) THEN
326 *
327 * Get machine constants.
328 *
329  safmin = pdlamch( desca( ctxt_ ), 'Safe minimum' )
330  eps = pdlamch( desca( ctxt_ ), 'Precision' )
331  smlnum = safmin / eps
332  bignum = one / smlnum
333  rmin = sqrt( smlnum )
334  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
335 *
336  nprocs = nprow*npcol
337  nb_a = desca( nb_ )
338  mb_a = desca( mb_ )
339  nb = nb_a
340  lower = lsame( uplo, 'L' )
341 *
342  rsrc_a = desca( rsrc_ )
343  csrc_a = desca( csrc_ )
344  iroffa = mod( ia-1, mb_a )
345  icoffa = mod( ja-1, nb_a )
346  iarow = indxg2p( 1, nb_a, myrow, rsrc_a, nprow )
347  iacol = indxg2p( 1, mb_a, mycol, csrc_a, npcol )
348  np0 = numroc( n+iroffa, nb, myrow, iarow, nprow )
349  nq0 = numroc( n+icoffa, nb, mycol, iacol, npcol )
350  IF( wantz ) THEN
351  nb_z = descz( nb_ )
352  mb_z = descz( mb_ )
353  rsrc_z = descz( rsrc_ )
354  iroffz = mod( iz-1, mb_a )
355  izrow = indxg2p( 1, nb_a, myrow, rsrc_z, nprow )
356  ELSE
357  iroffz = 0
358  izrow = 0
359  END IF
360 *
361 * COMPLEX*16 work space for PZHETRD
362 *
363  CALL pzhetrd( uplo, n, a, ia, ja, desca, rwork( indd ),
364  $ rwork( inde ), work( indtau ),
365  $ work( indwork ), -1, iinfo )
366  sizepzhetrd = int( abs( work( 1 ) ) )
367 *
368 * COMPLEX*16 work space for PZUNMTR
369 *
370  IF( wantz ) THEN
371  CALL pzunmtr( 'L', uplo, 'N', n, n, a, ia, ja, desca,
372  $ work( indtau ), z, iz, jz, descz,
373  $ work( indwork ), -1, iinfo )
374  sizepzunmtr = int( abs( work( 1 ) ) )
375  ELSE
376  sizepzunmtr = 0
377  END IF
378 *
379 * DOUBLE PRECISION work space for ZSTEQR2
380 *
381  IF( wantz ) THEN
382  rsizezsteqr2 = min( 1, 2*n-2 )
383  ELSE
384  rsizezsteqr2 = 0
385  END IF
386 *
387 * Initialize the context of the single column distributed
388 * matrix required by ZSTEQR2. This specific distribution
389 * allows each process to do 1/pth of the work updating matrix
390 * Q during ZSTEQR2 and achieve some parallelization to an
391 * otherwise serial subroutine.
392 *
393  ldc = 0
394  IF( wantz ) THEN
395  contextc = sl_gridreshape( desca( ctxt_ ), 0, 1, 1,
396  $ nprocs, 1 )
397  CALL blacs_gridinfo( contextc, nprowc, npcolc, myprowc,
398  $ mypcolc )
399  nrc = numroc( n, nb_a, myprowc, 0, nprocs )
400  ldc = max( 1, nrc )
401  CALL descinit( descqr, n, n, nb, nb, 0, 0, contextc, ldc,
402  $ info )
403  END IF
404 *
405 * COMPLEX*16 work space for ZSTEQR2
406 *
407  IF( wantz ) THEN
408  sizezsteqr2 = n*ldc
409  ELSE
410  sizezsteqr2 = 0
411  END IF
412 *
413 * Set up pointers into the WORK array
414 *
415  indtau = 1
416  indd = indtau + n
417  inde = indd + n
418  indwork = inde + n
419  indwork2 = indwork + n*ldc
420  llwork = lwork - indwork + 1
421 *
422 * Set up pointers into the RWORK array
423 *
424  indre = 1
425  indrd = indre + n
426  indrwork = indrd + n
427  llrwork = lrwork - indrwork + 1
428 *
429 * Compute the total amount of space needed
430 *
431  lrwmin = 2*n + rsizezsteqr2
432  lwmin = 3*n + max( sizepzhetrd, sizepzunmtr, sizezsteqr2 )
433 *
434  END IF
435  IF( info.EQ.0 ) THEN
436  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
437  info = -1
438  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
439  info = -2
440  ELSE IF( lwork.LT.lwmin .AND. lwork.NE.-1 ) THEN
441  info = -14
442  ELSE IF( lrwork.LT.lrwmin .AND. lrwork.NE.-1 ) THEN
443  info = -16
444  ELSE IF( iroffa.NE.0 ) THEN
445  info = -5
446  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
447  info = -( 700+nb_ )
448  END IF
449  IF( wantz ) THEN
450  IF( iroffa.NE.iroffz ) THEN
451  info = -10
452  ELSE IF( iarow.NE.izrow ) THEN
453  info = -10
454  ELSE IF( desca( m_ ).NE.descz( m_ ) ) THEN
455  info = -( 1200+m_ )
456  ELSE IF( desca( n_ ).NE.descz( n_ ) ) THEN
457  info = -( 1200+n_ )
458  ELSE IF( desca( mb_ ).NE.descz( mb_ ) ) THEN
459  info = -( 1200+mb_ )
460  ELSE IF( desca( nb_ ).NE.descz( nb_ ) ) THEN
461  info = -( 1200+nb_ )
462  ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) ) THEN
463  info = -( 1200+rsrc_ )
464  ELSE IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
465  info = -( 1200+ctxt_ )
466  END IF
467  END IF
468  END IF
469  IF( wantz ) THEN
470  idum1( 1 ) = ichar( 'V' )
471  ELSE
472  idum1( 1 ) = ichar( 'N' )
473  END IF
474  idum2( 1 ) = 1
475  IF( lower ) THEN
476  idum1( 2 ) = ichar( 'L' )
477  ELSE
478  idum1( 2 ) = ichar( 'U' )
479  END IF
480  idum2( 2 ) = 2
481  IF( lwork.EQ.-1 ) THEN
482  idum1( 3 ) = -1
483  ELSE
484  idum1( 3 ) = 1
485  END IF
486  idum2( 3 ) = 3
487  IF( wantz ) THEN
488  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 7, n, 3, n, 3, iz,
489  $ jz, descz, 12, 3, idum1, idum2, info )
490  ELSE
491  CALL pchk1mat( n, 3, n, 3, ia, ja, desca, 7, 3, idum1,
492  $ idum2, info )
493  END IF
494  work( 1 ) = dcmplx( lwmin )
495  rwork( 1 ) = dble( lrwmin )
496  END IF
497 *
498  IF( info.NE.0 ) THEN
499  CALL pxerbla( desca( ctxt_ ), 'PZHEEV', -info )
500  IF( wantz )
501  $ CALL blacs_gridexit( contextc )
502  RETURN
503  ELSE IF( lwork.EQ.-1 .OR. lrwork.EQ.-1 ) THEN
504  IF( wantz )
505  $ CALL blacs_gridexit( contextc )
506  RETURN
507  END IF
508 *
509 * Scale matrix to allowable range, if necessary.
510 *
511  iscale = 0
512 *
513  anrm = pzlanhe( 'M', uplo, n, a, ia, ja, desca,
514  $ rwork( indrwork ) )
515 *
516 *
517  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
518  iscale = 1
519  sigma = rmin / anrm
520  ELSE IF( anrm.GT.rmax ) THEN
521  iscale = 1
522  sigma = rmax / anrm
523  END IF
524 *
525  IF( iscale.EQ.1 ) THEN
526  CALL pzlascl( uplo, one, sigma, n, n, a, ia, ja, desca, iinfo )
527  END IF
528 *
529 * Reduce Hermitian matrix to tridiagonal form.
530 *
531  CALL pzhetrd( uplo, n, a, ia, ja, desca, rwork( indrd ),
532  $ rwork( indre ), work( indtau ), work( indwork ),
533  $ llwork, iinfo )
534 *
535 * Copy the values of D, E to all processes.
536 *
537  DO 10 i = 1, n
538  CALL pzelget( 'A', ' ', work( indd+i-1 ), a, i+ia-1, i+ja-1,
539  $ desca )
540  rwork( indrd+i-1 ) = dble( work( indd+i-1 ) )
541  10 CONTINUE
542  IF( lsame( uplo, 'U' ) ) THEN
543  DO 20 i = 1, n - 1
544  CALL pzelget( 'A', ' ', work( inde+i-1 ), a, i+ia-1, i+ja,
545  $ desca )
546  rwork( indre+i-1 ) = dble( work( inde+i-1 ) )
547  20 CONTINUE
548  ELSE
549  DO 30 i = 1, n - 1
550  CALL pzelget( 'A', ' ', work( inde+i-1 ), a, i+ia, i+ja-1,
551  $ desca )
552  rwork( indre+i-1 ) = dble( work( inde+i-1 ) )
553  30 CONTINUE
554  END IF
555 *
556  IF( wantz ) THEN
557 *
558  CALL pzlaset( 'Full', n, n, czero, cone, work( indwork ), 1, 1,
559  $ descqr )
560 *
561 * ZSTEQR2 is a modified version of LAPACK's CSTEQR. The
562 * modifications allow each process to perform partial updates
563 * to matrix Q.
564 *
565  CALL zsteqr2( 'I', n, rwork( indrd ), rwork( indre ),
566  $ work( indwork ), ldc, nrc, rwork( indrwork ),
567  $ info )
568 *
569  CALL pzgemr2d( n, n, work( indwork ), 1, 1, descqr, z, ia, ja,
570  $ descz, contextc )
571 *
572  CALL pzunmtr( 'L', uplo, 'N', n, n, a, ia, ja, desca,
573  $ work( indtau ), z, iz, jz, descz,
574  $ work( indwork ), llwork, iinfo )
575 *
576  ELSE
577 *
578  CALL zsteqr2( 'N', n, rwork( indrd ), rwork( indre ),
579  $ work( indwork ), 1, 1, rwork( indrwork ), info )
580  END IF
581 *
582 * Copy eigenvalues from workspace to output array
583 *
584  CALL dcopy( n, rwork( indd ), 1, w, 1 )
585 *
586 * If matrix was scaled, then rescale eigenvalues appropriately.
587 *
588  IF( iscale.EQ.1 ) THEN
589  CALL dscal( n, one / sigma, w, 1 )
590  END IF
591 *
592  work( 1 ) = dble( lwmin )
593 *
594 * Free up resources
595 *
596  IF( wantz ) THEN
597  CALL blacs_gridexit( contextc )
598  END IF
599 *
600 * Compare every ith eigenvalue, or all if there are only a few,
601 * across the process grid to check for heterogeneity.
602 *
603  IF( n.LE.ithval ) THEN
604  j = n
605  k = 1
606  ELSE
607  j = n / ithval
608  k = ithval
609  END IF
610 *
611  lrmin = int( rwork( 1 ) )
612  indtau = 0
613  inde = indtau + j
614  DO 40 i = 1, j
615  rwork( i+indtau ) = w( ( i-1 )*k+1 )
616  rwork( i+inde ) = w( ( i-1 )*k+1 )
617  40 CONTINUE
618 *
619  CALL dgamn2d( desca( ctxt_ ), 'All', ' ', j, 1, rwork( 1+indtau ),
620  $ j, 1, 1, -1, -1, 0 )
621  CALL dgamx2d( desca( ctxt_ ), 'All', ' ', j, 1, rwork( 1+inde ),
622  $ j, 1, 1, -1, -1, 0 )
623 *
624  DO 50 i = 1, j
625  IF( info.EQ.0 .AND. ( rwork( i+indtau )-rwork( i+inde ).NE.
626  $ zero ) ) THEN
627  info = n + 1
628  END IF
629  50 CONTINUE
630  rwork( 1 ) = lrmin
631 *
632  RETURN
633 *
634 * End of PZHEEV
635 *
636  END
pzelget
subroutine pzelget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pzelget.f:2
max
#define max(A, B)
Definition: pcgemr.c:180
zsteqr2
subroutine zsteqr2(COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO)
Definition: zsteqr2.f:2
pzlascl
subroutine pzlascl(TYPE, CFROM, CTO, M, N, A, IA, JA, DESCA, INFO)
Definition: pzlascl.f:3
pzlaset
subroutine pzlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pzblastst.f:7509
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
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pzheev
subroutine pzheev(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pzheev.f:3
pzunmtr
subroutine pzunmtr(SIDE, UPLO, TRANS, M, N, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pzunmtr.f:3
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
pzhetrd
subroutine pzhetrd(UPLO, N, A, IA, JA, DESCA, D, E, TAU, WORK, LWORK, INFO)
Definition: pzhetrd.f:3
min
#define min(A, B)
Definition: pcgemr.c:181