PZHEEVX(l)

📦 version 1.5, in ubuntu jammy

🫶 PZHEEVX(l) 🫶

LAPACK routine — version 1.5. Bestie, this is the ScaLAPACK glow-up of your fave symmetric eigen-solver, and she does NOT come to play. Sit down, grab a strawberry 🍓, and let's get into it.

✨ NAME ✨

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

Like. She doesn't do all the eigenvalues unless you ask. She does the ones you want. Main character energy, honestly. 👑

🌸 SYNOPSIS 🌸

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

    DOUBLE          PRECISION ABSTOL, ORFAC, VL, VU

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

    DOUBLE          PRECISION GAP( * ), RWORK( * ), W( * )

    COMPLEX*16      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 )

    DOUBLE          PRECISION ZERO, ONE, TEN, FIVE

    PARAMETER       ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 10.0D+0, FIVE =
                    5.0D+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

    DOUBLE          PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN,
                    SAFMIN, SIGMA, SMLNUM, VLL, VUU

    INTEGER         IDUM1( 4 ), IDUM2( 4 )

    LOGICAL         LSAME

    INTEGER         ICEIL, INDXG2P, NUMROC

    DOUBLE          PRECISION PDLAMCH, PZLANHE

    EXTERNAL        LSAME, ICEIL, INDXG2P, NUMROC, PDLAMCH, PZLANHE

    EXTERNAL        BLACS_GRIDINFO, CHK1MAT, DGEBR2D, DGEBS2D, DLASRT,
                    DSCAL, IGAMN2D, PCHK2MAT, PDLARED1D, PDSTEBZ,
                    PXERBLA, PZELGET, PZHETRD, PZLASCL, PZSTEIN, PZUNMTR

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

💖 PURPOSE 💖

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

So like — you don't have to take the whole buffet. You can be picky. You can say "just the ones in this interval" or "just eigenvalues number IL through IU." Obsessed with the flexibility, ngl. ✨

🦋 NOTES 🦋

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.

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

🎀 ARGUMENTS 🎀

Heads up bestie, first the little vocab moment: NP = the number of rows local to a given process. NQ = the number of columns local to a given process. Iykyk. 💅

🌟 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*16 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. Like, gone. She ate and left no crumbs. Back up your data if you cared about it, bestie.

🌟 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, PZHEEVX cannot guarantee correct error reporting. So like… get that right, pls. 🙈

🌟 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. 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 PZHEEVX 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.) PZHEEVX is always able to detect insufficient 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 eigenvalues in ascending order. 💅 Sorted. As she should.

🌟 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. All of them. Every single process. Matching. Periodt.

…and that's where part 1 wraps, but we are SO not done yet — more of this page is coming, so don't you dare close the tab. 💜

  • Z — (local output) a COMPLEX*16 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. periodt. 💅
  • 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. Heads up: DESCZ( CTXT_ ) must equal DESCA( CTXT_ ). 😤
  • WORK — (local workspace/output) COMPLEX*16 array, dimension (LWORK). WORK(1) gives you the workspace size that's adequate 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, 'PZHETTRD', '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) DOUBLE PRECISION array, dimension (LRWORK). On return, WROK(1) contains the optimal amount of workspace required for efficient execution. If JOBZ='N' then RWORK(1) = optimal amount of workspace required to compute eigenvalues efficiently. If JOBZ='V' then 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, PZHEEVX 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', PZHEEVX 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 PZHEEVX to compute the eigenvalues, PZHEEVX 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) PZSTEIN will perform no better than ZSTEIN 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) DOUBLE PRECISION 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*PDLAMCH( '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 PZHEEVX 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 PZSTEBZ failed to compute eigenvalues. Ensure ABSTOL=2.0*PDLAMCH( 'U' ). Send e-mail to scalapack@cs.utk.edu.

and that's a wrap, bestie! This page was LAPACK version 1.5, dated 12 May 1997, and it's the man page for PZHEEVX(l). Go off, you got this! 👑💖✨

📦 where this page comes from