ppsvx(3)

💅 ppsvx(3) — ppsvx: factor and solve, expert ✨

okay so bestie, welcome to ppsvx(3) — this is the LAPACK page for ppsvx: factor and solve, expert. it's giving "I know what I'm doing and I brought receipts" energy. let's get into it 🫶

---

🌸 NAME 🌸

ppsvx — ppsvx: factor and solve, expert 🎀

---

✨ SYNOPSIS ✨

Functions 🦋

• subroutine cppsvx (fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info) — CPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
• subroutine dppsvx (fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info) — DPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
• subroutine sppsvx (fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info) — SPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
• subroutine zppsvx (fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info) — ZPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

notice the pattern? c = complex, d = double, s = single, z = double-complex. the real ones pass rwork, the real ones that are single/double pass iwork. the complex ones pass rwork. it's giving "family reunion, everyone shows up" 💀

---

🌈 Detailed Description 🌈

this page is where the actual explanations live. buckle up, it's long but iconic 👑

👑 Function Documentation 👑

subroutine cppsvx (character fact, character uplo, integer n, integer nrhs, complex, dimension( * ) ap, complex, dimension( * ) afp, character equed, real, dimension( * ) s, complex, dimension( ldb, * ) b, integer ldb, complex, dimension( ldx, * ) x, integer ldx, real rcond, real, dimension( * ) ferr, real, dimension( * ) berr, complex, dimension( * ) work, real, dimension( * ) rwork, integer info)

CPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices 💖

🎀 Purpose:

!>
!> CPPSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!> compute the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N Hermitian positive definite matrix stored in
!> packed format and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>

lowkey that's the whole pitch: it solves A * X = B with a Hermitian positive definite A, packed format, Cholesky on deck, and it also hands you error bounds plus a condition estimate. main character behavior, honestly 💅

🔥 Description:

!>
!> The following steps are performed:
!>

1. Equilibrate (maybe) — if FACT = 'E', real scaling factors are computed to equilibrate the system:

diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B

whether or not the system actually gets equilibrated depends on the scaling of the matrix A. but if equilibration is used, A is overwritten by diag(S)*A*diag(S) and B by diag(S)*B. 😭 (yes, they get modified, we'll get to that)

2. Factor — if FACT = 'N' or 'E', the Cholesky decomposition factors the matrix A (after equilibration if FACT = 'E') as:

A = U**H * U ,  if UPLO = 'U', or
A = L * L**H,  if UPLO = 'L',

where U is upper triangular, L is lower triangular, and **H means conjugate transpose. periodt.

3. Check the minors — if the leading principal minor of order i is not positive, the routine returns with INFO = i. otherwise, the factored form of A is used to estimate the condition number of A. if the reciprocal of the condition number is less than machine precision, INFO = N+1 is returned as a warning — but like, the routine still goes on to solve for X and compute error bounds as described below. it's not a quitter 🫶
4. Solve — the system of equations is solved for X using the factored form of A.
5. Refine — iterative refinement is applied to improve the computed solution matrix, and to calculate error bounds and backward error estimates for it. glow-up round ✨
6. Un-equilibrate — if equilibration was used, the matrix X is premultiplied by diag(S) so that it solves the original system before equilibration. we put it back how we found it, bestie 🙈

---

🧸 Parameters 🧸

FACT

!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AFP contains the factored form of A.
!>                  If EQUED = 'Y', the matrix A has been equilibrated
!>                  with scaling factors given by S.  AP and AFP will not
!>                  be modified.
!>          = 'N':  The matrix A will be copied to AFP and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AFP and factored.
!>

so FACT is the "did you already do the work or not" switch:

• = 'F': you're bringing the factored form of A in AFP on entry. if EQUED = 'Y', A has already been equilibrated with scaling factors from S. AP and AFP will not be modified 🎀
• = 'N': A gets copied to AFP and factored
• = 'E': A gets equilibrated if necessary, then copied to AFP and factored

UPLO

!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>

which triangle of A are we storing? that one. 'U' = upper, 'L' = lower. choose wisely 🦋

N

!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>

the number of linear equations, i.e. the order of the matrix A. and it's N >= 0, no negative sizes, that's not a vibe.

NRHS

!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>

the number of right hand sides = the number of columns of B and X. also NRHS >= 0. bestie said no negatives in this house 🚫

AP

!>          AP is COMPLEX array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array, except if FACT = 'F'
!>          and EQUED = 'Y', then A must contain the equilibrated matrix
!>          diag(S)*A*diag(S).  The j-th column of A is stored in the
!>          array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.  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 is overwritten by
!>          diag(S)*A*diag(S).
!>

okay this one's the ✨packed storage✨ moment. AP is COMPLEX with dimension (N*(N+1)/2) — that's the whole triangle, no zero-padding, we're saving memory like it's our job (it is).

on entry it holds the upper or lower triangle of the Hermitian matrix A, packed columnwise in a linear array. except — and this is the iconic exception — if FACT = 'F' and EQUED = 'Y', then A must contain the equilibrated matrix diag(S)*A*diag(S).

the j-th column of A lives in AP like this:

if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

see below for further details (there's literally a whole section below, patience bestie 😭). 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 is overwritten by diag(S)*A*diag(S). slay but also watch out 💀

AFP

!>          AFP is COMPLEX array, dimension (N*(N+1)/2)
!>          If FACT = 'F', then AFP is an input argument and on entry
!>          contains the triangular factor U or L from the Cholesky
!>          factorization A = U**H*U or A = L*L**H, in the same storage
!>          format as A.  If EQUED .ne. 'N', then AFP is the factored
!>          form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AFP is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**H * U or A = L * L**H of the original
!>          matrix A.
!>
!>          If FACT = 'E', then AFP is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**H*U or A = L*L**H of the equilibrated
!>          matrix A (see the description of AP for the form of the
!>          equilibrated matrix).
!>

AFP is the factor-holding twin, also COMPLEX, also dimension (N*(N+1)/2). whether it's input or output is literally just vibes based on FACT:

• FACT = 'F' → input argument. on entry it contains the triangular factor U or L from the Cholesky factorization A = U**H*U or A = L*L**H, in the same storage format as A. if EQUED .ne. 'N', then AFP is the factored form of the equilibrated matrix A 🫶
• FACT = 'N' → output argument. on exit it returns the triangular factor U or L from the Cholesky factorization A = U**H * U or A = L * L**H of the original matrix A
• FACT = 'E' → output argument. on exit it returns the triangular factor U or L from the Cholesky factorization A = U**H*U or A = L*L**H of the equilibrated matrix A (see the description of AP for the form of the equilibrated matrix)

EQUED

!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'Y':  Equilibration was done, i.e., A has been replaced by
!>                  diag(S) * A * diag(S).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>

EQUED is CHARACTER*1 and it tells you what equilibration happened:

• = 'N': no equilibration (always true if FACT = 'N')
• = 'Y': equilibration was done, i.e. A has been replaced by diag(S) * A * diag(S)

and here's the gossip: EQUED is an input argument if FACT = 'F'; otherwise it's an output argument. it's giving "depends on the situation" 💅

S

!>          S is REAL array, dimension (N)
!>          The scale factors for A; not accessed if EQUED = 'N'.  S is
!>          an input argument if FACT = 'F'; otherwise, S is an output
!>          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
!>          must be positive.
!>

S is the REAL array of dimension (N) holding the scale factors for A. not accessed if EQUED = 'N' (if we're not equilibrating, we don't touch it, boundaries bestie 🦋).

S is an input argument if FACT = 'F'; otherwise it's an output argument. and important: if FACT = 'F' and EQUED = 'Y', each element of S must be positive. no zeroes, no negatives, positive only, periodt ⭐

B

!>          B is COMPLEX array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
!>          B is overwritten by diag(S) * B.
!>

B is COMPLEX, dimension (LDB,NRHS). on entry: the N-by-NRHS right hand side matrix B. on exit: if EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten by diag(S) * B. so like, check EQUED before you act surprised, ok 😭

LDB

!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>

LDB is INTEGER — the leading dimension of the array B. and LDB >= max(1,N). don't test it 🌸

X

!>          X is COMPLEX array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
!>          the original system of equations.  Note that if EQUED = 'Y',
!>          A and B are modified on exit, and the solution to the
!>          equilibrated system is inv(diag(S))*X.
!>

X is COMPLEX, dimension (LDX,NRHS). if INFO = 0 or INFO = N+1, this is the N-by-NRHS solution matrix X to the original system of equations 🎀

⚠️ note that if EQUED = 'Y', A and B are modified on exit, and the solution to the equilibrated system is inv(diag(S))*X. keep that in your notes bestie, it matters.

LDX

!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>

LDX is INTEGER, the leading dimension of the array X, and LDX >= max(1,N). same rules, different array 💜

RCOND

!>          RCOND is REAL
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>

RCOND is REAL: the estimate of the reciprocal condition number of the matrix A after equilibration (if done). if RCOND is less than machine precision (in particular, if RCOND = 0), the matrix is singular to working precision — and that condition is flagged by a return code of INFO > 0. so if you see that, it's not a glitch, it's a warning 🔥

FERR

!>          FERR is REAL array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>

FERR is a REAL array, dimension (NRHS) — the estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X).

if XTRUE is the true solution corresponding to X(j), then FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). the estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error. it's humble-bragging, honestly 🫶

BERR

!>          BERR is REAL array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>

BERR is a REAL array, dimension (NRHS) — the componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution). tiny change energy ✨

WORK

!>          WORK is COMPLEX array, dimension (2*N)
!>

WORK is a COMPLEX array, dimension (2*N). that's it, that's the whole parameter. it works so WORK can work 😌

RWORK

!>          RWORK is REAL array, dimension (N)
!>

RWORK is a REAL array, dimension (N). the real one next to the complex one. iconic duo 🎀

INFO

• INFO — INTEGER, the receipts on how it went 👀
  • = 0: successful exit ✨ slay, we're done
  • < 0: if INFO = -i, the i-th argument had an illegal value 🙈 (you gave it nonsense, bestie)
  • > 0: if INFO = i, and i is…
    • <= N: the leading principal minor of order i of A is not positive, so the factorization could not be completed, and the solution has not been computed. RCOND = 0 is returned 💀
    • = N+1: U is nonsingular, but RCOND is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest 🫶 (like, don't panic, it's still giving useful answers)

---

🧚 AUTHOR

• Univ. of Tennessee
• Univ. of California Berkeley
• Univ. of Colorado Denver
• NAG Ltd.

---

🎀 FURTHER DETAILS:

The packed storage scheme is illustrated by the following example when N = 4, UPLO = 'U' 👑

Two-dimensional storage of the Hermitian matrix A:

   a11 a12 a13 a14
       a22 a23 a24
           a33 a34     (aij = conjg(aji))
               a44

Packed storage of the upper triangle of A:

AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]

Definition at line 309 of file cppsvx.f.

---

💜 subroutine dppsvx (character fact, character uplo, integer n, integer nrhs, double precision, dimension( * ) ap, double precision, dimension( * ) afp, character equed, double precision, dimension( * ) s, double precision, dimension( ldb, * ) b, integer ldb, double precision, dimension( ldx, * ) x, integer ldx, double precision rcond, double precision, dimension( * ) ferr, double precision, dimension( * ) berr, double precision, dimension( * ) work, integer, dimension( * ) iwork, integer info)

subroutine dppsvx (character fact, character uplo, integer n, integer nrhs,
    double precision, dimension( * ) ap, double precision, dimension( * )
    afp, character equed, double precision, dimension( * ) s, double
    precision, dimension( ldb, * ) b, integer ldb, double precision,
    dimension( ldx, * ) x, integer ldx, double precision rcond, double
    precision, dimension( * ) ferr, double precision, dimension( * ) berr,
    double precision, dimension( * ) work, integer, dimension( * ) iwork,
    integer info)

DPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices 🌸

✨ Purpose:

DPPSVX 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 * X = B,

where A is an N-by-N symmetric positive definite matrix stored in packed format and X and B are N-by-NRHS matrices 🦋

Error bounds on the solution and a condition estimate are also provided. Yes, we're doing the whole glow-up: solution, error bounds, condition estimate, all of it 💅

🌈 Description:

The following steps are performed:

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

   diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * 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(S)*A*diag(S) and B by diag(S)*B.

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

3. If the leading principal minor of order i is not positive, then the routine returns with INFO = i. Otherwise, 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, INFO = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below 🫶
4. The system of equations is solved for X using the factored form of A.
5. Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. Bestie it literally polishes the answer, iconic 🔥
6. If equilibration was used, the matrix X is premultiplied by diag(S) so that it solves the original system before equilibration.

🧁 Parameters

• FACT — FACT is CHARACTER*1 — specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored.
  • = 'F': On entry, AFP contains the factored form of A. If EQUED = 'Y', the matrix A has been equilibrated with scaling factors given by S. AP and AFP will not be modified.
  • = 'N': The matrix A will be copied to AFP and factored.
  • = 'E': The matrix A will be equilibrated if necessary, then copied to AFP and factored.
• UPLO — UPLO is CHARACTER*1
  • = 'U': Upper triangle of A is stored;
  • = 'L': Lower triangle of A is stored.
• N — N is INTEGER — the number of linear equations, i.e., the order of the matrix A. N >= 0.
• NRHS — NRHS is INTEGER — the number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0.
• AP — AP is DOUBLE PRECISION array, dimension (N*(N+1)/2) — on entry, the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array, except if FACT = 'F' and EQUED = 'Y', then A must contain the equilibrated matrix diag(S)*A*diag(S). The j-th column of A is stored in the array AP as follows:

if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

See below for further details. 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 is overwritten by diag(S)*A*diag(S) 😭 she changes in place, main character energy

• AFP — AFP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
  • If FACT = 'F', then AFP is an input argument and on entry 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 AFP is the factored form of the equilibrated matrix A.
  • If FACT = 'N', then AFP is an output argument and on exit 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 AFP is an output argument and on exit 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 (see the description of AP for the form of the equilibrated matrix).
• EQUED — EQUED is CHARACTER*1 — specifies the form of equilibration that was done.
  • = 'N': No equilibration (always true if FACT = 'N').
  • = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).

EQUED is an input argument if FACT = 'F'; otherwise, it is an output argument 🎀 (yes, it switches sides, we love a flexible icon)

• S — S is DOUBLE PRECISION array, dimension (N) — the scale factors for A; not accessed if EQUED = 'N'. S is an input argument if FACT = 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED = 'Y', each element of S must be positive 🌸
• B — B is DOUBLE PRECISION array, dimension (LDB,NRHS) — on entry, the N-by-NRHS right hand side matrix B. On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten by diag(S) * B.
• LDB — LDB is INTEGER — the leading dimension of the array B. LDB >= max(1,N).
• X — X is DOUBLE PRECISION array, dimension (LDX,NRHS) — if INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to the original system of equations. Note that if EQUED = 'Y', A and B are modified on exit, and the solution to the equilibrated system is inv(diag(S))*X 👑
• LDX — LDX is INTEGER — the leading dimension of the array X. LDX >= max(1,N).
• RCOND — RCOND is DOUBLE PRECISION — the estimate of the reciprocal condition number of the matrix A after equilibration (if done). If RCOND is less than the machine precision (in particular, if RCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code of INFO > 0 💀
• FERR — FERR is DOUBLE PRECISION array, dimension (NRHS) — the estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error 🫶 (she's generous, not dramatic)
• BERR — BERR is DOUBLE PRECISION array, dimension (NRHS) — the componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).
• WORK
• WORK — a DOUBLE PRECISION array, dimension (3*N). That's it. That's the whole description. No notes, just vibes 💅
• IWORK — an INTEGER array, dimension (N). Same energy: name, type, size, done. Iconic minimalism honestly ✨
• INFO — INTEGER, and this one's the receipts bestie, so listen up 👑
  • = 0: successful exit 🎉
  • < 0: if INFO = -i, the i-th argument had an illegal value 🙈
  • > 0: if INFO = i, and i is:
    • <= N: the leading principal minor of order i of A is not positive, so the factorization could not be completed, and the solution has not been computed. RCOND = 0 is returned 💀
    • = N+1: U is nonsingular, but RCOND is less than machine precision, meaning the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest 🫶

Author

• Univ. of Tennessee
• Univ. of California Berkeley
• Univ. of Colorado Denver
• NAG Ltd.

Further Details:

The packed storage scheme is illustrated by the following example when N = 4, UPLO = 'U':

Two-dimensional storage of the symmetric matrix A:

   a11 a12 a13 a14
       a22 a23 a24
           a33 a34     (aij = conjg(aji))
               a44

Packed storage of the upper triangle of A:

AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]

Definition at line 309 of file dppsvx.f.

---

🌸 subroutine sppsvx

subroutine sppsvx (character fact, character uplo, integer n, integer nrhs,
    real, dimension( * ) ap, real, dimension( * ) afp, character equed,
    real, dimension( * ) s, real, dimension( ldb, * ) b, integer ldb, real,
    dimension( ldx, * ) x, integer ldx, real rcond, real, dimension( * )
    ferr, real, dimension( * ) berr, real, dimension( * ) work, integer,
    dimension( * ) iwork, integer info)

SPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices 🦋

🎯 Purpose

SPPSVX 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 * X = B, where A is an N-by-N symmetric positive definite matrix stored in packed format and X and B are N-by-NRHS matrices. Error bounds on the solution and a condition estimate are also provided. The whole package, no cap 🎁

🌈 Description

The following steps are performed, in this exact order, main character energy only:

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

   diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * 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(S)*A*diag(S) and B by diag(S)*B.

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

3. If the leading principal minor of order i is not positive, then the routine returns with INFO = i. Otherwise, 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, INFO = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. She doesn't give up, she just files a complaint 💅
4. The system of equations is solved for X using the factored form of A.
5. Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. Glow-up step ✨
6. If equilibration was used, the matrix X is premultiplied by diag(S) so that it solves the original system before equilibration. Undo the makeover so it matches the original pic 🎀

Parameters

• FACT — CHARACTER*1. Specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored.
  • = 'F': on entry, AFP contains the factored form of A. If EQUED = 'Y', the matrix A has been equilibrated with scaling factors given by S. AP and AFP will not be modified.
  • = 'N': the matrix A will be copied to AFP and factored.
  • = 'E': the matrix A will be equilibrated if necessary, then copied to AFP and factored.
• UPLO — CHARACTER*1.
  • = 'U': upper triangle of A is stored;
  • = 'L': lower triangle of A is stored.
• N — INTEGER. The number of linear equations, i.e., the order of the matrix A. N >= 0.
• NRHS — INTEGER. The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0.
• AP — REAL array, dimension (N*(N+1)/2). On entry, the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array, except if FACT = 'F' and EQUED = 'Y', then A must contain the equilibrated matrix diag(S)*A*diag(S). The j-th column of A is stored in the array AP as follows:
  • if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  • if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

See below for further details. 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 is overwritten by diag(S)*A*diag(S).

• AFP — REAL array, dimension (N*(N+1)/2).
  • If FACT = 'F', then AFP is an input argument and on entry 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 AFP is the factored form of the equilibrated matrix A.
  • If FACT = 'N', then AFP is an output argument and on exit 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 AFP is an output argument and on exit 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 (see the description of AP for the form of the equilibrated matrix).
• EQUED — CHARACTER*1. Specifies the form of equilibration that was done.
  • = 'N': no equilibration (always true if FACT = 'N').
  • = 'Y': equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).

EQUED is an input argument if FACT = 'F'; otherwise, it is an output argument. Sneaky little switcharoo, love that for her 🙈

• S — REAL array, dimension (N). The scale factors for A; not accessed if EQUED = 'N'. S is an input argument if FACT = 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED = 'Y', each element of S must be positive.
• B — REAL array, dimension (LDB,NRHS). On entry, the N-by-NRHS right hand side matrix B. On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten by diag(S) * B.
• LDB — INTEGER. The leading dimension of the array B. LDB >= max(1,N).
• X — REAL array, dimension (LDX,NRHS). If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to the original system of equations. Note that if EQUED = 'Y', A and B are modified on exit, and the solution to the equilibrated system is inv(diag(S))*X 🫶
• LDX — INTEGER. The leading dimension of the array X. LDX >= max(1,N).
• RCOND — REAL. The estimate of the reciprocal condition number of the matrix A after equilibration (if done). If RCOND is less than the machine precision (in particular, if RCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code of INFO > 0 🚨
• FERR — REAL array, dimension (NRHS). The estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error. Humble brag but make it math 💅
• BERR — REAL array, dimension (NRHS). The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).
• WORK — REAL array, dimension (3*N) 🌸

IWORK

IWORK is INTEGER array, dimension (N) — and that's genuinely the entire description, bestie. The docs give us nothing more, so I'm not gonna invent lore 💀 she's just an integer workspace of size N, living her truth.

INFO

INFO is INTEGER, and here's the whole gossip chain:

• = 0: successful exit ✨ everything ate, no crumbs.
• < 0: if INFO = -i, the i-th argument had an illegal value 🙈 (you passed something cursed).
• > 0: if INFO = i, and i is…
  • <= N: the leading principal minor of order i of A is not positive, so the factorization could not be completed, and the solution has not been computed. RCOND = 0 is returned. She said nope 💀
  • = N+1: U is nonsingular, but RCOND is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed, because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest. So like… it's giving ✨technically singular✨ but we ball anyway 🫶

---

⭐ Author ⭐

• Univ. of Tennessee
• Univ. of California Berkeley
• Univ. of Colorado Denver
• NAG Ltd.

🦋 Further Details:

Okay so picture this: the packed storage scheme, illustrated with N = 4, UPLO = 'U' — we're storing the upper triangle, no waste, very demure, very mindful 🎀

Two-dimensional storage of the symmetric matrix A:

   a11 a12 a13 a14
       a22 a23 a24
           a33 a34     (aij = conjg(aji))
               a44

And now the same thing, packed 💅 — packed storage of the upper triangle of A:

AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]

Notice how every element is accounted for, character for character? That's the icon behavior, no cap 👑

Definition at line 309 of file sppsvx.f.

---

🌈 subroutine zppsvx

subroutine zppsvx (character fact, character uplo, integer n, integer nrhs,
    complex*16, dimension( * ) ap, complex*16, dimension( * ) afp, character
    equed, double precision, dimension( * ) s, complex*16, dimension( ldb, *
    ) b, integer ldb, complex*16, dimension( ldx, * ) x, integer ldx, double
    precision rcond, double precision, dimension( * ) ferr, double
    precision, dimension( * ) berr, complex*16, dimension( * ) work, double
    precision, dimension( * ) rwork, integer info)

ZPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

🌸 Purpose:

ZPPSVX uses the Cholesky factorization A = U**H * U or A = L * L**H to compute the solution to a complex system of linear equations A * X = B, where A is an N-by-N Hermitian positive definite matrix stored in packed format and X and B are N-by-NRHS matrices. Error bounds on the solution and a condition estimate are also provided. The full package of care and attention, honestly 💖

💜 Description:

The following steps are performed, bestie, in order, no skipping:

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

diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * 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(S)*A*diag(S) and B by diag(S)*B.

2. If FACT = 'N' or 'E', the Cholesky decomposition is used to factor the matrix A (after equilibration if FACT = 'E') as

A = U**H * U ,  if UPLO = 'U', or
A = L * L**H,  if UPLO = 'L',

where U is an upper triangular matrix, L is a lower triangular matrix, and **H indicates conjugate transpose. Slay 💅

3. If the leading principal minor of order i is not positive, then the routine returns with INFO = i. Otherwise, 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, INFO = N+1 is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. (Warning but still functioning — main character energy, honestly.)
4. The system of equations is solved for X using the factored form of A.
5. Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
6. If equilibration was used, the matrix X is premultiplied by diag(S) so that it solves the original system before equilibration. She comes back home to the original coordinates 🧚

🎀 Parameters

FACT

FACT is CHARACTER*1 — specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored.

• = 'F': On entry, AFP contains the factored form of A. If EQUED = 'Y', the matrix A has been equilibrated with scaling factors given by S. AP and AFP will not be modified.
• = 'N': The matrix A will be copied to AFP and factored.
• = 'E': The matrix A will be equilibrated if necessary, then copied to AFP and factored.

UPLO

UPLO is CHARACTER*1:

• = 'U': Upper triangle of A is stored;
• = 'L': Lower triangle of A is stored.

N

N is INTEGER — the number of linear equations, i.e., the order of the matrix A. N >= 0.

NRHS

NRHS is INTEGER — the number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0.

AP

AP is COMPLEX*16 array, dimension (N*(N+1)/2).

On entry, the upper or lower triangle of the Hermitian matrix A, packed columnwise in a linear array, except if FACT = 'F' and EQUED = 'Y', then A must contain the equilibrated matrix diag(S)*A*diag(S). The j-th column of A is stored in the array AP as follows:

• if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
• if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

See below for further details. A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. (A stays unbothered, as she should 💅)

On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by diag(S)*A*diag(S).

AFP

AFP is COMPLEX*16 array, dimension (N*(N+1)/2).

• If FACT = 'F', then AFP is an input argument and on entry contains the triangular factor U or L from the Cholesky factorization A = U**H*U or A = L*L**H, in the same storage format as A. If EQUED .ne. 'N', then AFP is the factored form of the equilibrated matrix A.
• If FACT = 'N', then AFP is an output argument and on exit returns the triangular factor U or L from the Cholesky factorization A = U**H * U or A = L * L**H of the original matrix A.
• If FACT = 'E', then AFP is an output argument and on exit returns the triangular factor U or L from the Cholesky factorization A = U**H * U or A = L * L**H of the equilibrated matrix A (see the description of AP for the form of the equilibrated matrix).

EQUED

EQUED is CHARACTER*1 — specifies the form of equilibration that was done.

• = 'N': No equilibration (always true if FACT = 'N').
• = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).

EQUED is an input argument if FACT = 'F'; otherwise, it is an output argument. Know which one you're dealing with, bestie 🫶

S

S is DOUBLE PRECISION array, dimension (N) — the scale factors for A; not accessed if EQUED = 'N'.

S is an input argument if FACT = 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED = 'Y', each element of S must be positive. Positive only, we don't do negative scale factors in this house 🙅‍♀️

B

B is COMPLEX*16 array, dimension (LDB,NRHS).

On entry, the N-by-NRHS right hand side matrix B.

On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten by diag(S) * B.

LDB

LDB is INTEGER — the leading dimension of the array B. LDB >= max(1,N).

X

X is COMPLEX*16 array, dimension (LDX,NRHS).

If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to the original system of equations. Note that if EQUED = 'Y', A and B are modified on exit, and the solution to the equilibrated system is inv(diag(S))*X. Fold that piece of lore into your brain, it matters 😭

LDX

LDX is INTEGER — the leading dimension of the array X. LDX >= max(1,N).

RCOND

RCOND is DOUBLE PRECISION — the estimate of the reciprocal condition number of the matrix A after equilibration (if done).

If RCOND is less than the machine precision (in particular, if RCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code of INFO > 0. So if you see that, she's telling you something 🚨

FERR

FERR is DOUBLE PRECISION array, dimension (NRHS) — the estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X).

If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error. Overcautious queen, but we respect her 💅

BERR

• BERR — it’s a DOUBLE PRECISION array, dimension (NRHS). The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution). 💅
• WORK — a COMPLEX*16 array, dimension (2*N). ✨
• RWORK — a DOUBLE PRECISION array, dimension (N). 🌸
• INFO — an INTEGER. Here’s the whole mood:
  • = 0: successful exit
  • < 0: if INFO = -i, the i-th argument had an illegal value
  • > 0: if INFO = i, and i is:
    • <= N: the leading principal minor of order i of A is not positive, so the factorization could not be completed, and the solution has not been computed. RCOND = 0 is returned.
    • = N+1: U is nonsingular, but RCOND is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest.

---

✨ Author ✨

No cap, these legends wrote it:

• Univ. of Tennessee
• Univ. of California Berkeley
• Univ. of Colorado Denver
• NAG Ltd.

🌸 Further Details 🌸

The packed storage scheme is illustrated by the following example when N = 4, UPLO = 'U':

Two-dimensional storage of the Hermitian matrix A:

    a11 a12 a13 a14
        a22 a23 a24
            a33 a34     (aij = conjg(aji))
                a44

Packed storage of the upper triangle of A:

AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]

Definition at line 309 of file zppsvx.f. 📌

🤖 Author (Doxygen’s receipt) 🤖

Generated automatically by Doxygen for LAPACK from the source code. 🫶

---

Okay bestie, that’s the whole page! Footer check: LAPACK Version 3.12.0 and this iconic lil page is ppsvx(3) — that’s the source/version flex and page name, signed off with love. Periodt. 💖