PDSYGVX(l)

📦 version 1.5, in ubuntu jammy

hey bestie! welcome to the PDSYGVX(l) man page — that's a LAPACK routine, version 1.5, and she's literally the moment. no cap. let's get into it 💅✨

✨ NAME ✨

PDSYGVX — compute all the eigenvalues, and optionally, the eigenvectors of a real generalized SY-definite eigenproblem. periodt. 💖

🌸 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( * )

           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 )

           DOUBLE          PRECISION ONE

           PARAMETER       ( ONE = 1.0D+0 )

           DOUBLE          PRECISION FIVE, ZERO

           PARAMETER       ( FIVE = 5.0D+0, ZERO = 0.0D+0 )

           INTEGER         IERRNPD

           PARAMETER       ( IERRNPD = 16 )

           LOGICAL         ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ

           CHARACTER       TRANS

           INTEGER         IACOL, IAROW, IBCOL, IBROW, ICOFFA, ICOFFB, ICTXT,
                           IROFFA, IROFFB, LIWMIN, LWMIN, MQ0, MYCOL, MYROW,
                           NB, NEIG, NN, NP0, NPCOL, NPROW

           DOUBLE          PRECISION EPS, SCALE

           INTEGER         IDUM1( 5 ), IDUM2( 5 )

           LOGICAL         LSAME

           INTEGER         ICEIL, INDXG2P, NUMROC

           DOUBLE          PRECISION PDLAMCH

           EXTERNAL        LSAME, ICEIL, INDXG2P, NUMROC, PDLAMCH

           EXTERNAL        BLACS_GRIDINFO, CHK1MAT, DGEBR2D, DGEBS2D, DSCAL,
                           PCHK1MAT, PCHK2MAT, PDPOTRF, PDSYEVX, PDSYGST,
                           PDTRMM, PDTRSM, PXERBLA

           INTRINSIC       ABS, DBLE, ICHAR, MAX, MIN, MOD

🎀 PURPOSE 🎀

okay so PDSYGVX computes all the eigenvalues, and optionally, the eigenvectors, of a real generalized SY-definite eigenproblem. it's giving three possible forms:

  • sub( A )*x=(lambda)*sub( B )*x
  • sub( A )*sub( B )x=(lambda)*x
  • 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. slay. 👑

💅 NOTES 💅

so every 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. iykyk.

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".

  • DTYPE_A (global) — stored in DESCA( DTYPE_ ) — the descriptor type. in this case, DTYPE_A = 1.
  • CTXT_A (global) — stored in DESCA( CTXT_ ) — the BLACS context handle, indicating the BLACS process grid A is distributed over. the context itself is global, but the handle (the integer value) may vary.
  • M_A (global) — stored in DESCA( M_ ) — the number of rows in the global array A.
  • N_A (global) — stored in DESCA( N_ ) — the number of columns in the global array A.
  • MB_A (global) — stored in DESCA( MB_ ) — the blocking factor used to distribute the rows of the array.
  • NB_A (global) — stored in DESCA( NB_ ) — the blocking factor used to distribute the columns of the array.
  • RSRC_A (global) — stored in DESCA( RSRC_ ) — the process row over which the first row of the array A is distributed.
  • CSRC_A (global) — stored in DESCA( CSRC_ ) — the process column over which the first column of the array A is distributed.
  • LLD_A (local) — stored in 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. 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 👑

okay bestie, here come the arguments. pay attention, because they're all important. 💖

  • 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 leading 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 triangular 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 leading 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 triangular part of the matrix.

On exit, if INFO <= N, the part of sub( B ) containing the matrix 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 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*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. It's literally the number of columns of Z that are filled. If JOBZ .NE. 'V', NZ is not referenced. If JOBZ .EQ. 'V', NZ = M unless you supply insufficient space and PDSYGVX can't detect this before beginning computation. To get all the eigenvectors you requested, you 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 insufficient space without computation unless RANGE .EQ. 'V'. No cap. ✨
  • W — (global output) DOUBLE PRECISION array, dimension (N). On normal exit, the first M entries contain the selected eigenvalues 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. Bestie, keep it consistent! 💅
  • 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 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. 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 (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. This one is a whole saga, so buckle up:
    • 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. (We keep the typo "envionmental" because that's what the docs say, lol.) 🙈

  • 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 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) — DOUBLE PRECISION array, dimension (NPROW*NPCOL) 💖

Okay so this one is like the emotional support array for eigenvectors that just could NOT get their act together. It holds the gap between eigenvalues whose eigenvectors couldn't be reorthogonalized. The values you get back line up one-for-one with the clusters that ICLUSTR is pointing at, iykyk.

And here's the tea 🫖 — because of that, the dot product between eigenvectors corresponding to the I^th cluster can be as big as:

( C * n ) / GAP(I)

where C is a small constant. So yeah, tiny gap = big ol' overlap. It's giving "they're basically the same vector at this point" 💀

  • INFO (global output) — INTEGER ✨

The final boss of the output. Read it carefully bestie, it's telling you everything.

  • = 0: successful exit 🎉 slay, we ate, no notes.
  • < 0: an argument was illegal. 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. Math homework, but for bugs 🙈
  • > 0: bestie… something went sideways. And it's bitmask-coded, so get your calculator out 🧮 Each of these is a separate "oh no":
    • 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. Give it more memory next time, bestie 💅
    • 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. So like… partial slay, I guess 🥲
    • 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. Not the vibe, not the moment, absolutely not the matrix 🚫
fr fr though — that e-mail address is literally in the spec, so if INFO is throwing those bits at you, that's your move. No cap.

and that's the whole page bestie 🫶 go forth and PDSYGVX responsibly 💜

this iconic lil number comes to you from LAPACK version 1.5, dated 12 May 1997 — yeah she's vintage, she's archival, and she still runs circles around your fave 💅✨

📦 where this page comes from