PDSYGVX(3) ScaLAPACK routine of NEC Numeric Library Collection PDSYGVX(3)
NAME
PDSYGVX computes all the eigenvalues, and optionally, the eigenvectors
of a real generalized SY-definite eigenproblem.
SYNOPSIS
SUBROUTINE PDSYGVX( IBTYPE, JOBZ, RANGE, UPLO, N, A, IA, JA,
DESCA, B, IB, JB, DESCB, VL, VU, IL, IU, ABSTOL, M,
NZ, W, ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK,
LIWORK, IFAIL, ICLUSTR, GAP, INFO )
CHARACTER JOBZ, RANGE, UPLO
INTEGER IA, IB, IBTYPE, IL, INFO, IU, IZ, JA, JB, JZ,
LIWORK, LWORK, M, N, NZ
DOUBLE PRECISION ABSTOL, ORFAC, VL, VU
INTEGER DESCA( * ), DESCB( * ), DESCZ( * ), ICLUSTR( * ),
IFAIL( * ), IWORK( * )
DOUBLE PRECISION A( * ), B( * ), GAP( * ), W( * ), WORK( *
), Z( * )
PURPOSE
PDSYGVX computes all the eigenvalues, and optionally, the eigenvectors
of a real generalized SY-definite eigenproblem, of the form sub( A
)*x=(lambda)*sub( B )*x, sub( A )*sub( B )x=(lambda)*x, or sub( B
)*sub( A )*x=(lambda)*x. Here sub( A ) denoting A( IA:IA+N-1,
JA:JA+N-1 ) is assumed to be SY, and sub( B ) denoting B( IB:IB+N-1,
JB:JB+N-1 ) is assumed to be symmetric positive definite.
Notes
=====
Each global data object is described by an associated description vec-
tor. This vector stores the information required to establish the map-
ping 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 col-
umn.
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
ARGUMENTS
IBTYPE (global input) INTEGER
Specifies the problem type to be solved:
= 1: sub( A )*x = (lambda)*sub( B )*x
= 2: sub( A )*sub( B )*x = (lambda)*x
= 3: sub( B )*sub( A )*x = (lambda)*x
JOBZ (global input) CHARACTER*1
= '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
= 'U': Upper triangles of sub( A ) and sub( B ) are stored;
= 'L': Lower triangles of sub( A ) and sub( B ) are stored.
N (global input) INTEGER
The order of the matrices sub( A ) and sub( B ). N >= 0.
A (local input/local output) DOUBLE PRECISION pointer into the
local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
On entry, this array contains the local pieces of the N-by-N
symmetric distributed matrix sub( A ). If UPLO = 'U', the lead-
ing N-by-N upper triangular part of sub( A ) contains the upper
triangular part of the matrix. If UPLO = 'L', the leading N-
by-N lower triangular part of sub( A ) contains the lower tri-
angular part of the matrix.
On exit, if JOBZ = 'V', then if INFO = 0, sub( A ) contains the
distributed matrix Z of eigenvectors. The eigenvectors are
normalized as follows: if IBTYPE = 1 or 2, Z**T*sub( B )*Z = I;
if IBTYPE = 3, Z**T*inv( sub( B ) )*Z = I. If JOBZ = 'N', then
on exit the upper triangle (if UPLO='U') or the lower triangle
(if UPLO='L') of sub( A ), including the diagonal, is
destroyed.
IA (global input) INTEGER
The row index in the global array A indicating the first row of
sub( A ).
JA (global input) INTEGER
The column index in the global array A indicating the first
column of sub( A ).
DESCA (global and local input) INTEGER array of dimension DLEN_.
The array descriptor for the distributed matrix A. If DESCA(
CTXT_ ) is incorrect, PDSYGVX cannot guarantee correct error
reporting.
B (local input/local output) DOUBLE PRECISION pointer into the
local memory to an array of dimension (LLD_B, LOCc(JB+N-1)).
On entry, this array contains the local pieces of the N-by-N
symmetric distributed matrix sub( B ). If UPLO = 'U', the lead-
ing N-by-N upper triangular part of sub( B ) contains the upper
triangular part of the matrix. If UPLO = 'L', the leading N-
by-N lower triangular part of sub( B ) contains the lower tri-
angular part of the matrix.
On exit, if INFO <= N, the part of sub( B ) containing thema-
trix is overwritten by the triangular factor U or L from the
Cholesky factorization sub( B ) = U**T*U or sub( B ) = L*L**T.
IB (global input) INTEGER
The row index in the global array B indicating the first row of
sub( B ).
JB (global input) INTEGER
The column index in the global array B indicating the first
column of sub( B ).
DESCB (global and local input) INTEGER array of dimension DLEN_.
The array descriptor for the distributed matrix B. DESCB(
CTXT_ ) must equal DESCA( CTXT_ )
VL (global input) DOUBLE PRECISION
If RANGE='V', the lower bound of the interval to be searched
for eigenvalues. Not referenced if RANGE = 'A' or 'I'.
VU (global input) DOUBLE PRECISION
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) DOUBLE PRECISION
If JOBZ='V', setting ABSTOL to PDLAMCH( CONTEXT, 'U') yields
the most orthogonal eigenvectors.
The absolute error tolerance for the eigenvalues. An approxi-
mate 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*PDLAMCH('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*PDLAMCH('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 num-
ber 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 PDSYGVX 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.) PDSYGVX is always able to detect insuffi-
cient space without computation unless RANGE .EQ. 'V'.
W (global output) DOUBLE PRECISION array, dimension (N)
On normal exit, the first M entries contain the selected eigen-
values in ascending order.
ORFAC (global input) DOUBLE PRECISION
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.
ORFAC should be identical on all processes.
Z (local output) DOUBLE PRECISION 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 correspond-
ing to the selected eigenvalues. If an eigenvector fails to
converge, then that column of Z contains the latest approxima-
tion 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
The row index in the global array Z indicating the first row of
sub( Z ).
JZ (global input) INTEGER
The column index in the global array Z indicating the first
column of sub( Z ).
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) DOUBLE PRECISION array,
dimension max(3,LWORK)
if JOBZ='N' WORK(1) = optimal amount of workspace
required to compute eigenvalues efficiently
if JOBZ='V' WORK(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.
LWORK (local input) INTEGER
See below for definitions of variables used to define LWORK.
If no eigenvectors are requested (JOBZ = 'N') then
LWORK >= 5 * N + MAX( 5 * NN, NB * ( NP0 + 1 ) )
If eigenvectors are requested (JOBZ = 'V' ) then
the amount of workspace required to guarantee that all
eigenvectors are computed is:
LWORK >= 5 * N + MAX( 5*NN, NP0 * MQ0 + 2 * NB * NB ) +
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 LWORK:
(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 LWORK is too small:
If LWORK is too small to guarantee orthogonality,
PDSYGVX attempts to maintain orthogonality in
the clusters with the smallest
spacing between the eigenvalues.
If LWORK is too small to compute all the eigenvectors
requested, no computation is performed and INFO=-23
is returned. Note that when RANGE='V', PDSYGVX does
not know how many eigenvectors are requested until
the eigenvalues are computed. Therefore, when RANGE='V'
and as long as LWORK is large enough to allow PDSYGVX to
compute the eigenvalues, PDSYGVX will compute the
eigenvalues and as many eigenvectors as it can.
Relationship between workspace, orthogonality & performance:
Greater performance can be achieved if adequate workspace
is provided. On the other hand, in some situations,
performance can decrease as the workspace provided
increases above the workspace amount shown below:
For optimal performance, greater workspace may be
needed, i.e.
LWORK >= MAX( LWORK, 5 * N + NSYTRD_LWOPT,
NSYGST_LWOPT )
Where:
LWORK, as defined previously, depends upon the number
of eigenvectors requested, and
NSYTRD_LWOPT = N + 2*( ANB+1 )*( 4*NPS+2 ) +
( NPS + 3 ) * NPS
NSYGST_LWOPT = 2*NP0*NB + NQ0*NB + NB*NB
ANB = PJLAENV( DESCA( CTXT_), 3, 'PDSYTTRD', 'L',
0, 0, 0, 0)
SQNPC = INT( SQRT( DBLE( NPROW * NPCOL ) ) )
NPS = MAX( NUMROC( N, 1, 0, 0, SQNPC ), 2*ANB )
NB = DESCA( MB_ )
NP0 = NUMROC( N, NB, 0, 0, NPROW )
NQ0 = NUMROC( N, NB, 0, 0, NPCOL )
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.
For large N, no extra workspace is needed, however the
biggest boost in performance comes for small N, so it
is wise to provide the extra workspace (typically less
than a Megabyte per process).
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)
PDSTEIN will perform no better than DSTEIN 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 LWORK = -1, then LWORK is global input and a workspace query
is assumed; the routine only calculates the size required for
optimal performance on 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.
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 (output) INTEGER array, dimension (N)
IFAIL provides additional information when INFO .NE. 0 If
(MOD(INFO/16,2).NE.0) then IFAIL(1) indicates the order of the
smallest minor which is not positive definite. If
(MOD(INFO,2).NE.0) on exit, then IFAIL contains the indices of
the eigenvectors that failed to converge.
If neither of the above error conditions hold and JOBZ = 'V',
then the first M elements of IFAIL are set to zero.
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
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) DOUBLE PRECISION array, dimension
(NPROW*NPCOL)
This array contains the gap between eigenvalues whose eigenvec-
tors 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 cor-
respoding 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. 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
PDSYGVX 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 PDSTEBZ failed to
compute eigenvalues.
Send e-mail to scalapack@cs.utk.edu
if (MOD(INFO/16,2).NE.0), then B was not positive
definite. IFAIL(1) indicates the order of
the smallest minor which is not positive definite.
ALIGNMENT REQUIREMENTS
The distributed submatrices A(IA:*, JA:*), C(IC:IC+M-1,JC:JC+N-1), and
B( IB:IB+N-1, JB:JB+N-1 ) must verify some alignment properties, namely
the following expressions should be true:
DESCA(MB_) = DESCA(NB_)
IA = IB = IZ
JA = IB = JZ
DESCA(M_) = DESCB(M_) =DESCZ(M_)
DESCA(N_) = DESCB(N_)= DESCZ(N_)
DESCA(MB_) = DESCB(MB_) = DESCZ(MB_)
DESCA(NB_) = DESCB(NB_) = DESCZ(NB_)
DESCA(RSRC_) = DESCB(RSRC_) = DESCZ(RSRC_)
DESCA(CSRC_) = DESCB(CSRC_) = DESCZ(CSRC_)
MOD( IA-1, DESCA( MB_ ) ) = 0
MOD( JA-1, DESCA( NB_ ) ) = 0
MOD( IB-1, DESCB( MB_ ) ) = 0
MOD( JB-1, DESCB( NB_ ) ) = 0
ScaLAPACK routine 31 October 2017 PDSYGVX(3)