PSPOSVX(l)

📦 version 1.5, in ubuntu jammy

💖 PSPOSVX 💖

okay so this is a LAPACK routine, and the version we're looking at is version 1.5. The manual name is PSPOSVX(l) — that's the one, that's our girl. She solves linear systems by Cholesky and then tells you how much to trust the answer. Main character energy, fr fr. 🎀

🌸 NAME 🌸

PSPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to compute the solution to a real system of linear equations

A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),

like, that's the whole job bestie. 💅

🎀 SYNOPSIS 🎀

SUBROUTINE PSPOSVX( FACT, UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF,
                    DESCAF, EQUED, SR, SC, B, IB, JB, DESCB, X, IX, JX,
                    DESCX, RCOND, FERR, BERR, WORK, LWORK, IWORK,
                    LIWORK, INFO )

    CHARACTER       EQUED, FACT, UPLO

    INTEGER         IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
                    LWORK, N, NRHS

    REAL            RCOND

    INTEGER         DESCA( * ), DESCAF( * ), DESCB( * ), DESCX( * ),
                    IWORK( * )

    REAL            A( * ), AF( * ), B( * ), BERR( * ), FERR( * ), SC( *
                    ), SR( * ), WORK( * ), X( * )

look at all those arguments — she's high maintenance but she's worth it, iykyk. ✨

🦋 PURPOSE 🦋

So PSPOSVX takes the Cholesky factorization A = U**T*U or A = L*L**T and uses it to compute the solution to a real system of linear equations

A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),

where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix, and X and B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices. 🫶

She also hands you error bounds on the solution and a condition estimate. Not just solving it, she's grading it. Iconic.

Also, in all the comments that follow, Y denotes Y(IY:IY+M-1,JY:JY+K-1), a M-by-K matrix, where Y can be A, AF, B or X. 🌈

📝 Notes 📝

Every global data object has an associated description vector bestie. That vector stores all the info you need to map an element of the object to its process and its memory location. It's giving ✨directory✨.

Let A be a generic term for any 2D block cyclicly distributed array. That global array has an associated description vector DESCA. In the comments below, read the character _ 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)).

Now, let K be the number of rows or columns of a distributed matrix, and assume its process grid has dimension p x q.

LOCr( K ) is 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 ) is the number of elements of K that a process would receive if K were distributed over the q processes of its process row. 🦋

You can get LOCr() and LOCc() by calling 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 ).

And if you want an upper bound for these quantities, that's:

LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

math but make it cute 💀✨

👑 DESCRIPTION 👑

Here's the whole routine, step by step, no notes skipped, bestie:

  1. If FACT = 'E', real scaling factors are computed to equilibrate the system:

diag(SR) * A * diag(SC) * inv(diag(SC)) * X = diag(SR) * B

Whether or not the system will be equilibrated depends on the scaling of the matrix A, but if equilibration is used, A is overwritten by diag(SR)*A*diag(SC) and B by diag(SR)*B. 😳 (yes she overwrites your stuff, that's the deal)

  1. If FACT = 'N' or 'E', the Cholesky decomposition is used to factor the matrix A (after equilibration if FACT = 'E') as
A = U**T* U,  if UPLO = 'U', or
A = L * L**T,  if UPLO = 'L',

where U is an upper triangular matrix and L is a lower triangular matrix.

  1. The factored form of A is used to estimate the condition number of the matrix A. If the reciprocal of the condition number is less than machine precision, steps 4-6 are skipped. 🫢
  2. The system of equations is solved for X using the factored form of A.
  3. Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. ✨ glow up ✨
  4. If equilibration was used, the matrix X is premultiplied by diag(SR) so that it solves the original system before equilibration. 🎀

more of this page is coming bestie, we're not done yet 👑

💅 ARGUMENTS 💅

okay bestie, deep breath, this is the long part but it's literally just the cast list. every single one of these is a character and they all matter, no skipping 🫶

  • FACT (global input) CHARACTER — tells us whether the factored form of matrix A is already being handed to us on entry, and if it's not, whether we should equilibrate A before factoring it. the options are:
    • 'F': on entry, AF holds the factored form of A. if EQUED = 'Y', then A has been equilibrated with the scaling factors given by S. A and AF will not be modified. untouched, periodt.
    • 'N': A gets copied to AF and then factored. we're doing the work ourselves.
    • 'E': A gets equilibrated if it needs to be, then copied to AF and factored.
  • UPLO (global input) CHARACTER — which triangle of A are we even looking at:
    • 'U': upper triangle of A is stored.
    • 'L': lower triangle of A is stored.
  • N (global input) INTEGER — the number of rows and columns we're operating on, i.e. the order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1). constraint: N >= 0.
  • NRHS (global input) INTEGER — the number of right hand sides, i.e. the number of columns of the distributed submatrices B and X. constraint: NRHS >= 0.
  • A (local input/local output) REAL pointer into the local memory to an array of local dimension ( LLD_A, LOCc(JA+N-1) ). on entry, this is the symmetric matrix A — except if FACT = 'F' and EQUED = 'Y', in which case A must contain the equilibrated matrix diag(SR)*A*diag(SC). if UPLO = 'U', the leading N-by-N upper triangular part of A holds the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. if UPLO = 'L', the leading N-by-N lower triangular part of A holds the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.

on exit, if FACT = 'E' and EQUED = 'Y', A gets overwritten by diag(SR)*A*diag(SC). glowed up and scaled 💅

  • IA (global input) INTEGER — the row index in the global array A marking the first row of sub( A ).
  • JA (global input) INTEGER — the column index in the global array A marking the first column of sub( A ).
  • DESCA (global and local input) INTEGER array of dimension DLEN_ — the array descriptor for the distributed matrix A.
  • AF (local input or local output) REAL pointer into the local memory to an array of local dimension ( LLD_AF, LOCc(JA+N-1)). here's the tea depending on FACT:
    • if FACT = 'F', then AF is an input argument, and on entry it contains the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T, in the same storage format as A. if EQUED .ne. 'N', then AF is the factored form of the equilibrated matrix diag(SR)*A*diag(SC).
    • if FACT = 'N', then AF is an output argument, and on exit it returns the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T of the original matrix A.
    • if FACT = 'E', then AF is an output argument, and on exit it returns the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T of the equilibrated matrix A (peep the description of A for what that equilibrated matrix looks like).
  • IAF (global input) INTEGER — the row index in the global array AF marking the first row of sub( AF ).
  • JAF (global input) INTEGER — the column index in the global array AF marking the first column of sub( AF ).
  • DESCAF (global and local input) INTEGER array of dimension DLEN_ — the array descriptor for the distributed matrix AF.
  • EQUED (global input/global output) CHARACTER — says what form of equilibration went down:
    • 'N': no equilibration (always true if FACT = 'N').
    • 'Y': equilibration was done, i.e. A has been replaced by diag(SR) * A * diag(SC).

EQUED is an input variable if FACT = 'F'; otherwise it's an output variable. she switches roles, iconic.

  • SR (local input/local output) REAL array, dimension (LLD_A). the scale factors for A distributed across process rows; not accessed if EQUED = 'N'. SR is an input variable if FACT = 'F'; otherwise SR is an output variable. if FACT = 'F' and EQUED = 'Y', every element of SR must be positive.
  • SC (local input/local output) REAL array, dimension (LOC(N_A)). the scale factors for A distributed across process columns; not accessed if EQUED = 'N'. SC is an input variable if FACT = 'F'; otherwise SC is an output variable. if FACT = 'F' and EQUED = 'Y', every element of SC must be positive.
  • B (local input/local output) REAL pointer into the local memory to an array of local dimension ( LLD_B, LOCc(JB+NRHS-1) ). on entry, the N-by-NRHS right-hand side matrix B. on exit: if EQUED = 'N', B is not modified; if TRANS = 'N' and EQUED = 'R' or 'B', B gets overwritten by diag(R)*B; if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B gets overwritten by diag(C)*B.
  • IB (global input) INTEGER — the row index in the global array B marking the first row of sub( B ).
  • JB (global input) INTEGER — the column index in the global array B marking the first column of sub( B ).
  • DESCB (global and local input) INTEGER array of dimension DLEN_ — the array descriptor for the distributed matrix B.
  • X (local input/local output) REAL pointer into the local memory to an array of local dimension ( LLD_X, LOCc(JX+NRHS-1) ). if INFO = 0, this is the N-by-NRHS solution matrix X to the original system of equations. heads up: A and B are modified on exit if EQUED .ne. 'N', and the solution to the equilibrated system is inv(diag(SC))*X if TRANS = 'N' and EQUED = 'C' or or 'B'. (yes, that "or or" is in there, i'm just the messenger 🫠)
  • IX (global input) INTEGER — the row index in the global array X marking the first row of sub( X ).
  • JX (global input) INTEGER — the column index in the global array X marking the first column of sub( X ).
  • DESCX (global and local input) INTEGER array of dimension DLEN_ — the array descriptor for the distributed matrix X.
  • RCOND (global output) REAL — the estimate of the reciprocal condition number of the matrix A after equilibration (if that happened). if RCOND is less than the machine precision (especially if RCOND = 0), the matrix is singular to working precision. that situation gets flagged by a return code of INFO > 0, and the solution and error bounds are not computed. so like. nothing gets computed. it's giving "we tried" 😭
  • FERR (local output) REAL array, dimension (LOC(N_B)) — the estimated forward error bounds for each solution vector X(j) (the j-th column of the solution matrix X). if XTRUE is the true solution, FERR(j) bounds the magnitude of the largest entry in (X(j) - XTRUE) divided by the magnitude of the largest entry in X(j). the quality of that error bound depends on the quality of the estimate of norm(inv(A)) computed in the code; if the estimate of norm(inv(A)) is accurate, the error bound is guaranteed. accuracy is load-bearing, fr fr.
  • BERR (local output) REAL array, dimension (LOC(N_B)) — the componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any entry of A or B that makes X(j) an exact solution).
  • WORK (local workspace/local output) REAL array, dimension (LWORK). on exit, WORK(1) returns the minimal and optimal LWORK.
  • LWORK (local or global input) INTEGER — the dimension of the array WORK. LWORK is local input and must be at least
LWORK = MAX( PSPOCON( LWORK ), PSPORFS( LWORK ) ) + LOCr( N_A ).
LWORK = 3*DESCA( LLD_ )

if LWORK = -1, then LWORK is global input and a workspace query is assumed; the routine only calculates the minimum and optimal size for all work arrays. each of those values comes back in the first entry of the corresponding work array, and no error message is issued by PXERBLA. sneaky little query mode, we love her 🧚

  • IWORK (local workspace/local output) INTEGER array, dimension (LIWORK). on exit, IWORK(1) returns the minimal and optimal LIWORK.
  • LIWORK (local or global input) INTEGER — the dimension of the array IWORK. LIWORK is local input and must be at least
LIWORK = DESCA( LLD_ )
LIWORK = LOCr(N_A).

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 those values comes back in the first entry of the corresponding work array, and no error message is issued by PXERBLA.

  • INFO (global output) INTEGER — the vibe report:
    • = 0: successful exit. slay. ✨
    • < 0: if INFO = -i, the i-th argument had an illegal value.
    • > 0: if INFO = i, and i is
      • <= N: if INFO = i, the leading minor of order i of A is not positive definite, so the factorization could not be completed, and the solution and error bounds could not be computed.
      • = N+1: RCOND is less than machine precision. the factorization has been completed, but the matrix is singular to working precision, and the solution and error bounds have not been computed.
okay that was a LOT of arguments but you ate that, bestie. every single one of them matters, no cap 💖

this whole page came from LAPACK version 1.5, dated 12 May 1997 — yes she's an icon from the 90s and she's still running the math world, we stan a vintage queen 👑🌸

📦 where this page comes from