PDSYEV(l)

📦 version 1.5, in ubuntu jammy

💖 PDSYEV(l) — LAPACK routine (version 1.3) 💖

okay bestie, gather round, because PDSYEV(l) is here and she is that girl. This is a LAPACK routine, version 1.3, and she computes all eigenvalues and, if you want them, eigenvectors of a real symmetric matrix A by calling the recommended sequence of ScaLAPACK routines. No cap. 💅

✨ NAME ✨

PDSYEV — computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A by calling the recommended sequence of ScaLAPACK routines. 🦋

🎀 SYNOPSIS 🎀

here's the exact signature + all the declarations, copied with love and zero edits:

       SUBROUTINE PDSYEV( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ,
                          WORK, LWORK, INFO )

           CHARACTER      JOBZ, UPLO

           INTEGER        IA, INFO, IZ, JA, JZ, LWORK, N

           INTEGER        DESCA( * ), DESCZ( * )

           DOUBLE         PRECISION A( * ), 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 ZERO, ONE

           PARAMETER      ( ZERO = 0.0D+0, ONE = 1.0D+0 )

           INTEGER        ITHVAL

           PARAMETER      ( ITHVAL = 10 )

           LOGICAL        LOWER, WANTZ

           INTEGER        CONTEXTC, CSRC_A, I, IACOL, IAROW, ICOFFA, IINFO,
                          INDD, INDD2, INDE, INDE2, INDTAU, INDWORK, INDWORK2,
                          IROFFA, IROFFZ, ISCALE, IZROW, J, K, LCM, LCMQ, LDC,
                          LLWORK, LWMIN, MB_A, MB_Z, MYCOL, MYPCOLC, MYPROWC,
                          MYROW, NB, NB_A, NB_Z, NN, NP, NPCOL, NPCOLC, NPROCS,
                          NPROW, NPROWC, NQ, NRC, QRMEM, RSRC_A, RSRC_Z,
                          SIZEMQRLEFT, SIZEMQRRIGHT

           DOUBLE         PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
                          SIGMA, SMLNUM

           INTEGER        DESCQR( 10 ), IDUM1( 3 ), IDUM2( 3 )

           LOGICAL        LSAME

           INTEGER        ILCM, INDXG2P, NUMROC, SL_GRIDRESHAPE

           DOUBLE         PRECISION PDLAMCH, PDLANSY

           EXTERNAL       LSAME, ILCM, INDXG2P, NUMROC, SL_GRIDRESHAPE,
                          PDLAMCH, PDLANSY

           EXTERNAL       BLACS_GRIDEXIT, BLACS_GRIDINFO, CHK1MAT, DCOPY,
                          DESCINIT, DGAMN2D, DGAMX2D, DSCAL, DSTEQR2, PCHK2MAT,
                          PDELGET, PDGEMR2D, PDLASCL, PDLASET, PDORMTR,
                          PDSYTRD, PXERBLA

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

🌸 PURPOSE 🌸

okay so PDSYEV computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A by calling the recommended sequence of ScaLAPACK routines. That's the whole point, and she eats every time. 👑

⚠️ real talk warning: in its present form, PDSYEV assumes a homogeneous system and makes no checks for consistency of the eigenvalues or eigenvectors across the different processes. Because of this, it is possible that a heterogeneous system may return incorrect results without any error messages. so like, be careful bestie, because she will not warn you. 💀

🌈 NOTES 🌈

A description vector is associated with each 2D block-cyclicly distributed matrix. This vector stores the information required to establish the mapping between a matrix entry and its corresponding process and memory location. it's giving admin, it's giving mapping, it's giving necessary. ✨

In the following comments, the character _ should be read as "of the distributed matrix". Let A be a generic term for any 2D block cyclicly distributed matrix. Its description vector is DESCA:

  • DTYPE_A (global) — stored in DESCA( DTYPE_) — the descriptor type.
  • 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 distributed matrix A.
  • N_A (global) — stored in DESCA( N_ ) — the number of columns in the distributed matrix A.
  • MB_A (global) — stored in DESCA( MB_ ) — the blocking factor used to distribute the rows of A.
  • NB_A (global) — stored in DESCA( NB_ ) — the blocking factor used to distribute the columns of A.
  • RSRC_A (global) — stored in DESCA( RSRC_ ) — the process row over which the first row of the matrix A is distributed.
  • CSRC_A (global) — stored in DESCA( CSRC_ ) — the process column over which the first column of A is distributed.
  • LLD_A (local) — stored in DESCA( LLD_ ) — the leading dimension of the local array storing the local blocks of the distributed matrix A. 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 ).

✨ ARGUMENTS ✨

okay so this is the part where we meet every single argument, no cap. each one is doing its own lil thing and you need to know all of them, periodt. 💖

  • 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.
  • UPLO (global input) CHARACTER*1 — specifies whether the upper or lower triangular part of the symmetric 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. so if N is negative, that's a no from the routine, bestie.
  • A (local input/workspace) block cyclic DOUBLE PRECISION array, global dimension (N, N), local dimension ( LLD_A, LOCc(JA+N-1) )
    • On entry, the symmetric matrix A. If UPLO = 'U', only the upper triangular part of A is used to define the elements of the symmetric matrix. If UPLO = 'L', only the lower triangular part of A is used to define the elements of the symmetric 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. don't say we didn't warn you 😭
  • 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, PDSYEV cannot guarantee correct error reporting. so pls get your context right, iykyk.
  • W (global output) DOUBLE PRECISION array, dimension (N) — On normal exit, the first M entries contain the selected eigenvalues in ascending order. ascending, like your faves in a ranking. slay.
  • 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 JOBZ = 'N', then Z is not referenced. literally not even looked at. ghosted.
  • 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_ ) no exceptions, bestie. they gotta match.
  • WORK (local workspace/output) DOUBLE PRECISION array, dimension (LWORK)
    • Version 1.0: on output, WORK(1) returns the workspace needed to guarantee completion.
    • If the input parameters are incorrect, WORK(1) may also be incorrect.
    • If JOBZ='N' WORK(1) = minimal=optimal amount of workspace
    • If JOBZ='V' WORK(1) = minimal workspace required to generate all the eigenvectors.
  • 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 + SIZESYTRD + 1

where

SIZESYTRD = The workspace requirement for PDSYTRD
            and is MAX( NB * ( NP +1 ), 3 * NB )

If eigenvectors are requested (JOBZ = 'V' ) then the amount of workspace required to guarantee that all eigenvectors are computed is:

QRMEM = 2*N-2
LWMIN = 5*N + N*LDC + MAX( SIZEMQRLEFT, QRMEM ) + 1

Variable definitions:

NB = DESCA( MB_ ) = DESCA( NB_ ) =
     DESCZ( MB_ ) = DESCZ( NB_ )
NN = MAX( N, NB, 2 )
DESCA( RSRC_ ) = DESCA( RSRC_ ) = DESCZ( RSRC_ ) =
                 DESCZ( CSRC_ ) = 0
NP = NUMROC( NN, NB, 0, 0, NPROW )
NQ = NUMROC( MAX( N, NB, 2 ), NB, 0, 0, NPCOL )
NRC = NUMROC( N, NB, MYPROWC, 0, NPROCS)
LDC = MAX( 1, NRC )
SIZEMQRLEFT = The workspace requirement for PDORMTR
              when it's SIDE argument is 'L'.

With MYPROWC defined when a new context is created as:

CALL BLACS_GET( DESCA( CTXT_ ), 0, CONTEXTC )
CALL BLACS_GRIDINIT( CONTEXTC, 'R', NPROCS, 1 )
CALL BLACS_GRIDINFO( CONTEXTC, NPROWC, NPCOLC, MYPROWC, MYPCOLC )

If LWORK = -1, the LWORK is global input and a workspace query is assumed; the routine only calculates the minimum size for the WORK array. The required workspace is returned as the first element of WORK and no error message is issued by PXERBLA. so if you're just asking "how much space do i need", set it to -1 and vibe.

  • INFO (global output) INTEGER
    • = 0: successful exit we love a zero, it's giving success 💅
    • < 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 INFO = 1 through N, the i(th) eigenvalue did not converge in DSTEQR2 after a total of 30*N iterations. If INFO = N+1, then PDSYEV has detected heterogeneity by finding that eigenvalues were not identical across the process grid. In this case, the accuracy of the results from PDSYEV cannot be guaranteed. so if you see N+1, bestie, the accuracy is lowkey not guaranteed. be careful.

and that's the whole argument list, no cap. every single one accounted for, periodt. 💖

this iconic page is from LAPACK version 1.3, dated 12 May 1997. so vintage, so slay, mwah 💋✨

📦 where this page comes from