PSSYGVX(l)

📦 version 1.5, in ubuntu jammy

PSSYGVX 💖 the LAPACK routine, version 1.5, manual page PSSYGVX(l)

okay bestie, sit down, hydrate, because we're about to get into a whole generalized symmetric-definite eigenvalue situation and it's honestly giving main character energy 🎀

✨ NAME ✨

PSSYGVX — computes all the eigenvalues, and optionally, the eigenvectors of a real generalized SY-definite eigenproblem. no cap. 💅

🌸 SYNOPSIS 🌸

SUBROUTINE PSSYGVX( 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

    REAL            ABSTOL, ORFAC, VL, VU

    INTEGER         DESCA( * ), DESCB( * ), DESCZ( * ), ICLUSTR( * ),
                    IFAIL( * ), IWORK( * )

    REAL            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 )

    REAL            ONE

    PARAMETER       ( ONE = 1.0E+0 )

    REAL            FIVE, ZERO

    PARAMETER       ( FIVE = 5.0E+0, ZERO = 0.0E+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

    REAL            EPS, SCALE

    INTEGER         IDUM1( 5 ), IDUM2( 5 )

    LOGICAL         LSAME

    INTEGER         ICEIL, INDXG2P, NUMROC

    REAL            PSLAMCH

    EXTERNAL        LSAME, ICEIL, INDXG2P, NUMROC, PSLAMCH

    EXTERNAL        BLACS_GRIDINFO, CHK1MAT, PCHK1MAT, PCHK2MAT,
                    PSPOTRF, PSSYEVX, PSSYGST, PSTRMM, PSTRSM, PXERBLA,
                    SGEBR2D, SGEBS2D, SSCAL

    INTRINSIC       ABS, ICHAR, MAX, MIN, MOD, REAL

💜 PURPOSE 💜

PSSYGVX computes all the eigenvalues, and optionally, the eigenvectors of a real generalized SY-definite eigenproblem, and the forms it handles are, and I quote, fr fr:

  • sub( A )*x=(lambda)*sub( B )*x
  • sub( A )*sub( B )x=(lambda)*x
  • sub( B )*sub( A )*x=(lambda)*x

Here sub( A ) — that's A( IA:IA+N-1, JA:JA+N-1 ) — is assumed to be SY, and sub( B ) — that's B( IB:IB+N-1, JB:JB+N-1 ) — is assumed to be symmetric positive definite. periodt. 💅

🦋 NOTES 🦋

okay so every global data object gets described by an associated description vector. that vector is literally holding the info you need to map an element of the object onto its process and its memory location. iconic behavior. ✨

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". (yes, the underscore is doing a little cosplay. let it. 🎀)

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

Now 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 ) = how many elements of K a process would receive if K were distributed over the p processes of its process column.
  • LOCc( K ) = how many elements of K 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, and it's literally a vibe selector:
    • = 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. Do we want the eigenvectors or not, bestie?
    • = 'N': Compute eigenvalues only;
    • = 'V': Compute eigenvalues and eigenvectors.
  • RANGE — (global input) CHARACTER*1. Which eigenvalues are we hunting? 👑
    • = '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. Which triangle are we storing, babes? 🌸
    • = '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 ). Constraint: N >= 0.
  • A — (local input/local output) REAL 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. 😭 gone. reduced to atoms.
  • 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. Heads up bestie: if DESCA( CTXT_ ) is incorrect, PSSYGVX cannot guarantee correct error reporting. 💀
  • B — (local input/local output) REAL 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. Hard requirement, no negotiating: DESCB( CTXT_ ) must equal DESCA( CTXT_ ). 💅
  • 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'. she's off duty in those cases. 🙈
  • 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. Constraint: 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. Constraint: min(IL,N) <= IU <= N. Not referenced if RANGE = 'A' or 'V'.
  • ABSTOL — (global input) REAL. Listen, this one matters so pay attention bestie 🫶
    • If JOBZ='V', setting ABSTOL to PSLAMCH( CONTEXT, 'U') yields the most orthogonal eigenvectors. slay.
    • 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. 📚
part 1 out of 3, we're not done yet, keep scrolling bestie 💖
📎 See "On the correctness of Parallel Bisection in Floating Point" by Demmel, Dhillon and Ren, LAPACK Working Note #70

🎀 M — the eigenvalue headcount

Type: (global output) INTEGER

Total number of eigenvalues found. Range check bestie: 0 <= M <= N. That's the whole vibe, no notes. 💅

🦋 NZ — how many eigenvectors actually got built

Type: (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', then NZ is not referenced — don't even look at her.
  • If JOBZ .EQ. 'V', then NZ = M — unless you supplied insufficient space and PSSYGVX is not able to detect that before beginning computation.

To get all the eigenvectors you asked for, you must supply both of these:

  • sufficient space to hold the eigenvectors in Z (M .LE. DESCZ(N_)), and
  • sufficient workspace to compute them (peep LWORK below 👇).

PSSYGVX is always able to detect insufficient space without computation — unless RANGE .EQ. 'V'. Then it can't clock it in advance and that's just how it is. 😭

🌸 W — the actual eigenvalues

Type: (global output) REAL array, dimension (N)

On normal exit, the first M entries contain the selected eigenvalues in ascending order. Sorted, slay, iconic. ✨

✨ ORFAC — the reorthogonalization vibe check

Type: (global input) REAL

Specifies which eigenvectors should be reorthogonalized. Eigenvectors that correspond to eigenvalues sitting within tol=ORFAC*norm(A) of each other get reorthogonalized.

However, if the workspace is insufficient (see LWORK), tol may be decreased until all the eigenvectors to be reorthogonalized can be stored in one process. No cap.

  • If ORFAC equals zero, no reorthogonalization will be done. Chill.
  • If ORFAC is negative, a default value of 10^-3 is used.

ORFAC should be identical on all processes. Periodt.

🫶 Z — the eigenvector storage

Type: (local output) REAL 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 that eigenvector is returned in IFAIL.
  • If JOBZ = 'N', then Z is not referenced. She's not here, don't text her. 🙈

🧚 IZ

Type: (global input) INTEGER

The row index in the global array Z indicating the first row of sub( Z ).

🧚 JZ

Type: (global input) INTEGER

The column index in the global array Z indicating the first column of sub( Z ).

💜 DESCZ

Type: (global and local input) INTEGER array of dimension DLEN_.

The array descriptor for the distributed matrix Z. And heads up: DESCZ( CTXT_ ) must equal DESCA( CTXT_ ). Non-negotiable, fr fr.

🌈 WORK

Type: (local workspace/output) REAL 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 — the big one

Type: (local input) INTEGER

See below for definitions of the variables used to define LWORK. Grab a snack, bestie. 🍓

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. Not a drill.

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, PSSYGVX 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', PSSYGVX 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 PSSYGVX to compute the eigenvalues, PSSYGVX will compute the eigenvalues and as many eigenvectors as it can. She's doing her best, okay. 😭

🔥 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, 'PSSYTTRD', '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 function; PJLAENV is a ScaLAPACK environmental 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). Tiny cost, huge glow-up. ✨

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), PSSTEIN will perform no better than SSTEIN 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. Free real estate, we love a quiet exit. 🤫

🎀 IWORK

Type: (local workspace) INTEGER array

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

🎀 LIWORK

Type: (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 — the receipts

Type: (output) INTEGER array, dimension (N)

IFAIL provides additional information when INFO .NE. 0. The tea:

  • 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 — the clusters that didn't get reorthogonalized

Type: (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 — how far apart they were

Type: (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 corresponding to the I^th cluster may be as high as ( C * n ) / GAP(I) where C is a small constant.

  • INFO (global output) INTEGER — bestie, this is the ✨vibe check✨ on how the whole thing went. Read it like a group chat receipt 💅
    • = 0 — slay!! successful exit, we ate and left no crumbs 👑
    • < 0 — uh oh. 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. (yes, the numbers are doing math at you. that's the point. 🧮)
    • > 0 — ok so this one is the drama one, it's giving layered mystery 😭 the bits stack, so you gotta decode them one by one:
      • if (MOD(INFO,2).NE.0), then one or more eigenvectors failed to converge. their indices are stored in IFAIL. and like, actually — send e-mail to scalapack@cs.utk.edu if this happens, they wanna know 🫶
      • 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 PSSYGVX 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 PSSTEBZ failed to compute eigenvalues. send e-mail to scalapack@cs.utk.edu — yes again, they really do want the tea ☕
      • 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. 🌸
lowkey iconic that a single integer is out here carrying this much lore. read it carefully, don't skim it, iykyk 🦋

and that's the whole page, bestie!! PSSYGVX(l) is living in LAPACK version 1.5, dated 12 May 1997 — yeah she's vintage, she's a classic, she's been serving since before most of us could walk 💖✨ go forth and factor responsibly 🎀

📦 where this page comes from