NAME¶

PCHEEVX - compute selected eigenvalues and, optionally, eigenvectors of a complex hermitian matrix A by calling the recommended sequence of ScaLAPACK routines

SYNOPSIS¶

SUBROUTINE PCHEEVX(

JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, VL, VU, IL, IU, ABSTOL, M, NZ, W, ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO )

CHARACTER JOBZ, RANGE, UPLO INTEGER IA, IL, INFO, IU, IZ, JA, JZ, LIWORK, LRWORK, LWORK, M, N, NZ REAL ABSTOL, ORFAC, VL, VU INTEGER DESCA( * ), DESCZ( * ), ICLUSTR( * ), IFAIL( * ), IWORK( * ) REAL GAP( * ), RWORK( * ), W( * ) COMPLEX A( * ), WORK( * ), Z( * ) INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, MB_, NB_, RSRC_, CSRC_, LLD_ PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) REAL ZERO, ONE, TEN, FIVE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 10.0E+0, FIVE = 5.0E+0 ) INTEGER IERREIN, IERRCLS, IERRSPC, IERREBZ PARAMETER ( IERREIN = 1, IERRCLS = 2, IERRSPC = 4, IERREBZ = 8 ) LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, QUICKRETURN, VALEIG, WANTZ CHARACTER ORDER INTEGER CSRC_A, I, IACOL, IAROW, ICOFFA, IINFO, INDD, INDD2, INDE, INDE2, INDIBL, INDISP, INDRWORK, INDTAU, INDWORK, IROFFA, IROFFZ, ISCALE, ISIZESTEBZ, ISIZESTEIN, IZROW, LALLWORK, LIWMIN, LLRWORK, LLWORK, LRWMIN, LWMIN, MAXEIGS, MB_A, MB_Z, MQ0, MYCOL, MYROW, NB, NB_A, NB_Z, NEIG, NN, NNP, NP0, NPCOL, NPROCS, NPROW, NQ0, NSPLIT, NZZ, OFFSET, RSRC_A, RSRC_Z, SIZEHEEVX, SIZEORMTR, SIZESTEIN REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, VLL, VUU INTEGER IDUM1( 4 ), IDUM2( 4 ) LOGICAL LSAME INTEGER ICEIL, INDXG2P, NUMROC REAL PCLANHE, PSLAMCH EXTERNAL LSAME, ICEIL, INDXG2P, NUMROC, PCLANHE, PSLAMCH EXTERNAL BLACS_GRIDINFO, CHK1MAT, IGAMN2D, PCELGET, PCHETRD, PCHK2MAT, PCLASCL, PCSTEIN, PCUNMTR, PSLARED1D, PSSTEBZ, PXERBLA, SGEBR2D, SGEBS2D, SLASRT, SSCAL INTRINSIC ABS, CMPLX, ICHAR, MAX, MIN, MOD, REAL, SQRT

PURPOSE¶

PCHEEVX computes selected eigenvalues and, optionally, eigenvectors of a complex hermitian matrix A by calling the recommended sequence of ScaLAPACK routines. Eigenvalues/vectors can be selected by specifying a range of values or a range of indices for the desired eigenvalues.

NOTES¶

Each global data object is described by an associated description vector. This vector stores the information required to establish the mapping between an object element and its corresponding process and memory location.

Let A be a generic term for any 2D block cyclicly distributed array. Such a global array has an associated description vector DESCA. In the following comments, the character _ should be read as "of the global array".

NOTATION STORED IN EXPLANATION
--------------- -------------- --------------------------------------
DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
DTYPE_A = 1.
CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
the BLACS process grid A is distribu-
ted over. The context itself is glo-
bal, but the handle (the integer
value) may vary.
M_A (global) DESCA( M_ ) The number of rows in the global
array A.
N_A (global) DESCA( N_ ) The number of columns in the global
array A.
MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
the rows of the array.
NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
the columns of the array.
RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
row of the array A is distributed.
CSRC_A (global) DESCA( CSRC_ ) The process column over which the
first column of the array A is
distributed.
LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
array. LLD_A >= MAX(1,LOCr(M_A)).

Let K be the number of rows or columns of a distributed matrix, and assume that its process grid has dimension p x q. LOCr( K ) denotes the number of elements of K that a process would receive if K were distributed over the p processes of its process column.S
Similarly, LOCc( K ) denotes the number of elements of K that a process would receive if K were distributed over the q processes of its process row.
The values of LOCr() and LOCc() may be determined via a call to the ScaLAPACK tool function, NUMROC:
LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
An upper bound for these quantities may be computed by:
LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

PCHEEVX assumes IEEE 754 standard compliant arithmetic. To port to a system which does not have IEEE 754 arithmetic, modify the appropriate SLmake.inc file to include the compiler switch -DNO_IEEE. This switch only affects the compilation of pslaiect.c.

ARGUMENTS¶

NP = the number of rows local to a given process.
NQ = the number of columns local to a given process.

JOBZ (global input) CHARACTER*1
Specifies whether or not to compute the eigenvectors:
= 'N': Compute eigenvalues only.
= 'V': Compute eigenvalues and eigenvectors.

RANGE (global input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the interval [VL,VU] will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.

UPLO (global input) CHARACTER*1
Specifies whether the upper or lower triangular part of the
Hermitian matrix A is stored:
= 'U': Upper triangular
= 'L': Lower triangular

N (global input) INTEGER
The number of rows and columns of the matrix A. N >= 0.

A (local input/workspace) block cyclic COMPLEX array,
global dimension (N, N),
local dimension ( LLD_A, LOCc(JA+N-1) )

On entry, the Hermitian matrix A. If UPLO = 'U', only the
upper triangular part of A is used to define the elements of
the Hermitian matrix. If UPLO = 'L', only the lower
triangular part of A is used to define the elements of the
Hermitian matrix.

On exit, the lower triangle (if UPLO='L') or the upper
triangle (if UPLO='U') of A, including the diagonal, is
destroyed.

IA (global input) INTEGER
A's global row index, which points to the beginning of the
submatrix which is to be operated on.

JA (global input) INTEGER
A's global column index, which points to the beginning of
the submatrix which is to be operated on.

DESCA (global and local input) INTEGER array of dimension DLEN_.
The array descriptor for the distributed matrix A.
If DESCA( CTXT_ ) is incorrect, PCHEEVX cannot guarantee
correct error reporting.

VL (global input) REAL
If RANGE='V', the lower bound of the interval to be searched
for eigenvalues. Not referenced if RANGE = 'A' or 'I'.

VU (global input) REAL
If RANGE='V', the upper bound of the interval to be searched
for eigenvalues. Not referenced if RANGE = 'A' or 'I'.

IL (global input) INTEGER
If RANGE='I', the index (from smallest to largest) of the
smallest eigenvalue to be returned. IL >= 1.
Not referenced if RANGE = 'A' or 'V'.

IU (global input) INTEGER
If RANGE='I', the index (from smallest to largest) of the
largest eigenvalue to be returned. min(IL,N) <= IU <= N.
Not referenced if RANGE = 'A' or 'V'.

ABSTOL (global input) REAL
If JOBZ='V', setting ABSTOL to PSLAMCH( CONTEXT, 'U') yields
the most orthogonal eigenvectors.

The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to

ABSTOL + EPS * max( |a|,|b| ) ,

where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*norm(T) will be used in its place,
where norm(T) is the 1-norm of the tridiagonal matrix
obtained by reducing A to tridiagonal form.

Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*PSLAMCH('S') not zero.
If this routine returns with ((MOD(INFO,2).NE.0) .OR.
(MOD(INFO/8,2).NE.0)), indicating that some eigenvalues or
eigenvectors did not converge, try setting ABSTOL to
2*PSLAMCH('S').

See "Computing Small Singular Values of Bidiagonal Matrices
with Guaranteed High Relative Accuracy," by Demmel and
Kahan, LAPACK Working Note #3.

See "On the correctness of Parallel Bisection in Floating
Point" by Demmel, Dhillon and Ren, LAPACK Working Note #70

M (global output) INTEGER
Total number of eigenvalues found. 0 <= M <= N.

NZ (global output) INTEGER
Total number of eigenvectors computed. 0 <= NZ <= M.
The number of columns of Z that are filled.
If JOBZ .NE. 'V', NZ is not referenced.
If JOBZ .EQ. 'V', NZ = M unless the user supplies
insufficient space and PCHEEVX is not able to detect this
before beginning computation. To get all the eigenvectors
requested, the user must supply both sufficient
space to hold the eigenvectors in Z (M .LE. DESCZ(N_))
and sufficient workspace to compute them. (See LWORK below.)
PCHEEVX is always able to detect insufficient space without
computation unless RANGE .EQ. 'V'.

W (global output) REAL array, dimension (N)
On normal exit, the first M entries contain the selected
eigenvalues in ascending order.

ORFAC (global input) REAL
Specifies which eigenvectors should be reorthogonalized.
Eigenvectors that correspond to eigenvalues which are within
tol=ORFAC*norm(A) of each other are to be reorthogonalized.
However, if the workspace is insufficient (see LWORK),
tol may be decreased until all eigenvectors to be
reorthogonalized can be stored in one process.
No reorthogonalization will be done if ORFAC equals zero.
A default value of 10^-3 is used if ORFAC is negative.
ORFAC should be identical on all processes.

Z (local output) COMPLEX array,
global dimension (N, N),
local dimension ( LLD_Z, LOCc(JZ+N-1) )
If JOBZ = 'V', then on normal exit the first M columns of Z
contain the orthonormal eigenvectors of the matrix
corresponding to the selected eigenvalues. If an eigenvector
fails to converge, then that column of Z contains the latest
approximation to the eigenvector, and the index of the
eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.

IZ (global input) INTEGER
Z's global row index, which points to the beginning of the
submatrix which is to be operated on.

JZ (global input) INTEGER
Z's global column index, which points to the beginning of
the submatrix which is to be operated on.

DESCZ (global and local input) INTEGER array of dimension DLEN_.
The array descriptor for the distributed matrix Z.
DESCZ( CTXT_ ) must equal DESCA( CTXT_ )

WORK (local workspace/output) COMPLEX array,
dimension (LWORK)
WORK(1) returns workspace adequate workspace to allow
optimal performance.

LWORK (local input) INTEGER
Size of WORK array. If only eigenvalues are requested:
LWORK >= N + MAX( NB * ( NP0 + 1 ), 3 )
If eigenvectors are requested:
LWORK >= N + ( NP0 + MQ0 + NB ) * NB
with NQ0 = NUMROC( NN, NB, 0, 0, NPCOL ).

For optimal performance, greater workspace is needed, i.e.
LWORK >= MAX( LWORK, NHETRD_LWORK )
Where LWORK is as defined above, and
NHETRD_LWORK = N + 2*( ANB+1 )*( 4*NPS+2 ) +
( NPS + 1 ) * NPS

ICTXT = DESCA( CTXT_ )
ANB = PJLAENV( ICTXT, 3, 'PCHETTRD', 'L', 0, 0, 0, 0 )
SQNPC = SQRT( DBLE( NPROW * NPCOL ) )
NPS = MAX( NUMROC( N, 1, 0, 0, SQNPC ), 2*ANB )

NUMROC is a ScaLAPACK tool functions;
PJLAENV is a ScaLAPACK envionmental inquiry function
MYROW, MYCOL, NPROW and NPCOL can be determined by calling
the subroutine BLACS_GRIDINFO.

If LWORK = -1, then LWORK is global input and a workspace
query is assumed; the routine only calculates the
optimal size for all work arrays. Each of these
values is returned in the first entry of the corresponding
work array, and no error message is issued by PXERBLA.

RWORK (local workspace/output) REAL array,
dimension (LRWORK)
On return, WROK(1) contains the optimal amount of
workspace required for efficient execution.
if JOBZ='N' RWORK(1) = optimal amount of workspace
required to compute eigenvalues efficiently
if JOBZ='V' RWORK(1) = optimal amount of workspace
required to compute eigenvalues and eigenvectors
efficiently with no guarantee on orthogonality.
If RANGE='V', it is assumed that all eigenvectors
may be required.

LRWORK (local input) INTEGER
Size of RWORK
See below for definitions of variables used to define LRWORK.
If no eigenvectors are requested (JOBZ = 'N') then
LRWORK >= 5 * NN + 4 * N
If eigenvectors are requested (JOBZ = 'V' ) then
the amount of workspace required to guarantee that all
eigenvectors are computed is:
LRWORK >= 4*N + MAX( 5*NN, NP0 * MQ0 ) +
ICEIL( NEIG, NPROW*NPCOL)*NN

The computed eigenvectors may not be orthogonal if the
minimal workspace is supplied and ORFAC is too small.
If you want to guarantee orthogonality (at the cost
of potentially poor performance) you should add
the following to LRWORK:
(CLUSTERSIZE-1)*N
where CLUSTERSIZE is the number of eigenvalues in the
largest cluster, where a cluster is defined as a set of
close eigenvalues: { W(K),...,W(K+CLUSTERSIZE-1) |
W(J+1) <= W(J) + ORFAC*2*norm(A) }
Variable definitions:
NEIG = number of eigenvectors requested
NB = DESCA( MB_ ) = DESCA( NB_ ) =
DESCZ( MB_ ) = DESCZ( NB_ )
NN = MAX( N, NB, 2 )
DESCA( RSRC_ ) = DESCA( NB_ ) = DESCZ( RSRC_ ) =
DESCZ( CSRC_ ) = 0
NP0 = NUMROC( NN, NB, 0, 0, NPROW )
MQ0 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL )
ICEIL( X, Y ) is a ScaLAPACK function returning
ceiling(X/Y)

When LRWORK is too small:
If LRWORK is too small to guarantee orthogonality,
PCHEEVX attempts to maintain orthogonality in
the clusters with the smallest
spacing between the eigenvalues.
If LRWORK is too small to compute all the eigenvectors
requested, no computation is performed and INFO=-25
is returned. Note that when RANGE='V', PCHEEVX does
not know how many eigenvectors are requested until
the eigenvalues are computed. Therefore, when RANGE='V'
and as long as LRWORK is large enough to allow PCHEEVX to
compute the eigenvalues, PCHEEVX will compute the
eigenvalues and as many eigenvectors as it can.

Relationship between workspace, orthogonality & performance:
If CLUSTERSIZE >= N/SQRT(NPROW*NPCOL), then providing
enough space to compute all the eigenvectors
orthogonally will cause serious degradation in
performance. In the limit (i.e. CLUSTERSIZE = N-1)
PCSTEIN will perform no better than CSTEIN on 1
processor.
For CLUSTERSIZE = N/SQRT(NPROW*NPCOL) reorthogonalizing
all eigenvectors will increase the total execution time
by a factor of 2 or more.
For CLUSTERSIZE > N/SQRT(NPROW*NPCOL) execution time will
grow as the square of the cluster size, all other factors
remaining equal and assuming enough workspace. Less
workspace means less reorthogonalization but faster
execution.

If LRWORK = -1, then LRWORK is global input and a workspace
query is assumed; the routine only calculates the size
required for optimal performance for all work arrays. Each of
these values is returned in the first entry of the
corresponding work arrays, and no error message is issued by
PXERBLA.

IWORK (local workspace) INTEGER array
On return, IWORK(1) contains the amount of integer workspace
required.

LIWORK (local input) INTEGER
size of IWORK
LIWORK >= 6 * NNP
Where:
NNP = MAX( N, NPROW*NPCOL + 1, 4 )
If LIWORK = -1, then LIWORK is global input and a workspace
query is assumed; the routine only calculates the minimum
and optimal size for all work arrays. Each of these
values is returned in the first entry of the corresponding
work array, and no error message is issued by PXERBLA.

IFAIL (global output) INTEGER array, dimension (N)
If JOBZ = 'V', then on normal exit, the first M elements of
IFAIL are zero. If (MOD(INFO,2).NE.0) on exit, then
IFAIL contains the
indices of the eigenvectors that failed to converge.
If JOBZ = 'N', then IFAIL is not referenced.

ICLUSTR (global output) integer array, dimension (2*NPROW*NPCOL)
This array contains indices of eigenvectors corresponding to
a cluster of eigenvalues that could not be reorthogonalized
due to insufficient workspace (see LWORK, ORFAC and INFO).
Eigenvectors corresponding to clusters of eigenvalues indexed
ICLUSTR(2*I-1) to ICLUSTR(2*I), could not be
reorthogonalized due to lack of workspace. Hence the
eigenvectors corresponding to these clusters may not be
orthogonal. ICLUSTR() is a zero terminated array.
(ICLUSTR(2*K).NE.0 .AND. ICLUSTR(2*K+1).EQ.0) if and only if
K is the number of clusters
ICLUSTR is not referenced if JOBZ = 'N'

GAP (global output) REAL array,
dimension (NPROW*NPCOL)
This array contains the gap between eigenvalues whose
eigenvectors could not be reorthogonalized. The output
values in this array correspond to the clusters indicated
by the array ICLUSTR. As a result, the dot product between
eigenvectors correspoding to the I^th cluster may be as high
as ( C * n ) / GAP(I) where C is a small constant.

INFO (global output) INTEGER
= 0: successful exit
< 0: If the i-th argument is an array and the j-entry had
an illegal value, then INFO = -(i*100+j), if the i-th
argument is a scalar and had an illegal value, then
INFO = -i.
> 0: if (MOD(INFO,2).NE.0), then one or more eigenvectors
failed to converge. Their indices are stored
in IFAIL. Ensure ABSTOL=2.0*PSLAMCH( 'U' )
Send e-mail to scalapack@cs.utk.edu
if (MOD(INFO/2,2).NE.0),then eigenvectors corresponding
to one or more clusters of eigenvalues could not be
reorthogonalized because of insufficient workspace.
The indices of the clusters are stored in the array
ICLUSTR.
if (MOD(INFO/4,2).NE.0), then space limit prevented
PCHEEVX from computing all of the eigenvectors
between VL and VU. The number of eigenvectors
computed is returned in NZ.
if (MOD(INFO/8,2).NE.0), then PCSTEBZ failed to compute
eigenvalues. Ensure ABSTOL=2.0*PSLAMCH( 'U' )
Send e-mail to scalapack@cs.utk.edu