realOTHERsolve(3)

📦
version 3.11.0 · 2 versions
, in debian bookworm

💖 NAME 💖

realOTHERsolve — the real Other Solve Routines, bestie. Yes this is the realOTHERsolve(3) page, straight out of LAPACK, and we are about to get so linear-algebra-pilled together. Grab a boba, hydrate, and let's go. 🦋

✨ SYNOPSIS ✨

🎀 Functions 🎀

  • subroutine sgglse (M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, INFO) — SGGLSE solves overdetermined or underdetermined systems for OTHER matrices
  • subroutine spbsv (UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO) — SPBSV computes the solution to system of linear equations A * X = B for OTHER matrices
  • subroutine spbsvx (FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO) — SPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
  • subroutine sppsv (UPLO, N, NRHS, AP, B, LDB, INFO) — SPPSV 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 sspsv (UPLO, N, NRHS, AP, IPIV, B, LDB, INFO) — SSPSV computes the solution to system of linear equations A * X = B for OTHER matrices
  • subroutine sspsvx (FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO) — SSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

🌈 Detailed Description 🌈

okay so fr fr, this is the group of real Other Solve routines. That's the whole tea. It's a group, it's real, it's Other Solve, and it's serving exactly what it says on the tin. 🫶

👑 Function Documentation 👑

📐 subroutine sgglse 📐

subroutine sgglse (integer M, integer N, integer P, real, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, real, dimension( * ) C, real, dimension( * ) D, real, dimension( * ) X, real, dimension( * ) WORK, integer LWORK, integer INFO)

SGGLSE solves overdetermined or underdetermined systems for OTHER matrices. main character energy, honestly. 💅

🎯 Purpose

SGGLSE solves the linear equality-constrained least squares (LSE) problem:

minimize || c - A*x ||_2   subject to   B*x = d

where A is an M-by-N matrix, B is a P-by-N matrix, c is a given M-vector, and d is a given P-vector. It is assumed that P <= N <= M+P, and

                     rank(B) = P and  rank( (A) ) = N.
                                          ( (B) )

These conditions ensure that the LSE problem has a unique solution, which is obtained using a generalized RQ factorization of the matrices (B, A) given by

B = (0 R)*Q,   A = Z*T*Q.

🧾 Parameters

  • M — M is INTEGER. The number of rows of the matrix A. M >= 0.
  • N — N is INTEGER. The number of columns of the matrices A and B. N >= 0.
  • P — P is INTEGER. The number of rows of the matrix B. 0 <= P <= N <= M+P.
  • A — A is REAL array, dimension (LDA,N)
    • On entry, the M-by-N matrix A.
    • On exit, the elements on and above the diagonal of the array contain the min(M,N)-by-N upper trapezoidal matrix T.
  • LDA — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • B — B is REAL array, dimension (LDB,N)
    • On entry, the P-by-N matrix B.
    • On exit, the upper triangle of the subarray B(1:P,N-P+1:N) contains the P-by-P upper triangular matrix R.
  • LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,P).
  • C — C is REAL array, dimension (M)
    • On entry, C contains the right hand side vector for the least squares part of the LSE problem.
    • On exit, the residual sum of squares for the solution is given by the sum of squares of elements N-P+1 to M of vector C.
  • D — D is REAL array, dimension (P)
    • On entry, D contains the right hand side vector for the constrained equation.
    • On exit, D is destroyed. 😭 (gone, reduced to atoms)
  • X — X is REAL array, dimension (N). On exit, X is the solution of the LSE problem.
  • WORK — WORK is REAL array, dimension (MAX(1,LWORK)). On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  • LWORK — LWORK is INTEGER. The dimension of the array WORK. LWORK >= max(1,M+N+P). For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB, where NB is an upper bound for the optimal blocksizes for SGEQRF, SGERQF, SORMQR and SORMRQ.
    • If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.
  • INFO — INFO is INTEGER
    • = 0: successful exit.
    • < 0: if INFO = -i, the i-th argument had an illegal value.
    • = 1: the upper triangular factor R associated with B in the generalized RQ factorization of the pair (B, A) is singular, so that rank(B) < P; the least squares solution could not be computed.
    • = 2: the (N-P) by (N-P) part of the upper trapezoidal factor T associated with A in the generalized RQ factorization of the pair (B, A) is singular, so that
rank( (A) ) < N; the least squares solution could not
     ( (B) )
be computed.

✍️ Author

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

📐 subroutine spbsv 📐

subroutine spbsv (character UPLO, integer N, integer KD, integer NRHS, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldb, * ) B, integer LDB, integer INFO)

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

🎯 Purpose

SPBSV computes the solution to a real system of linear equations

A * X = B,

where A is an N-by-N symmetric positive definite band matrix and X and B are N-by-NRHS matrices.

The Cholesky decomposition is used to factor A as

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

where U is an upper triangular band matrix, and L is a lower triangular band matrix, with the same number of superdiagonals or subdiagonals as A. The factored form of A is then used to solve the system of equations A * X = B.

🧾 Parameters

  • 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.
  • KD — KD is INTEGER. The number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0.
  • NRHS — NRHS is INTEGER. The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
  • AB — AB is REAL array, dimension (LDAB,N)
    • On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows:
      • if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
      • if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD).
    • See below for further details.
    • On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T of the band matrix A, in the same storage format as A.
  • LDAB — LDAB is INTEGER. The leading dimension of the array AB. LDAB >= KD+1.
  • B — B is REAL array, dimension (LDB,NRHS)
    • On entry, the N-by-NRHS right hand side matrix B.
    • On exit, if INFO = 0, the N-by-NRHS solution matrix X.
  • LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
  • INFO — INFO is INTEGER
    • = 0: successful exit
    • < 0: if INFO = -i, the i-th argument had an illegal value
    • > 0: 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 has not been computed.

✍️ Author

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

🔍 Further Details

The band storage scheme is illustrated by the following example, when N = 6, KD = 2, and UPLO = 'U':

On entry:                       On exit:

    *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
    *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
   a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66

Similarly, if UPLO = 'L' the format of A is as follows:

On entry:                       On exit:

   a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
   a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
   a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *

Array elements marked * are not used by the routine.

📐 subroutine spbsvx 📐

subroutine spbsvx (character FACT, character UPLO, integer N, integer KD, integer NRHS, real, dimension( ldab, * ) AB, integer LDAB, real, dimension( ldafb, * ) AFB, integer LDAFB, 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)

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

🎯 Purpose

SPBSVX 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 band matrix and X and B are N-by-NRHS matrices.

Error bounds on the solution and a condition estimate are also provided.

okay bestie we are NOT done — this page keeps going, and so do we. part 2 incoming, don't touch that dial. 💖

✨ Description ✨

okay bestie, here's the whole glow-up sequence this routine runs through, step by step, no skips 💅

  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.

  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 band matrix, and L is a lower triangular band matrix.

  1. If the leading i-by-i principal minor is not positive definite, 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. it's giving "warning but make it work" 😭
  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(S) so that it solves the original system before equilibration.

🎀 Parameters 🎀

FACT

FACT is CHARACTER*1 — tells us whether the factored form of the matrix A is supplied on entry, and if it's not, whether A should be equilibrated before it gets factored.

  • = 'F': On entry, AFB contains the factored form of A. If EQUED = 'Y', the matrix A has been equilibrated with scaling factors given by S. AB and AFB will not be modified.
  • = 'N': The matrix A will be copied to AFB and factored.
  • = 'E': The matrix A will be equilibrated if necessary, then copied to AFB 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. bestie zero is allowed, it's just vibes at that point 💀

KD

KD is INTEGER — the number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 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.

AB

AB is REAL array, dimension (LDAB,N).

On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the 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 j-th column of the array AB as follows:

if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).

See below for further details 👀

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

LDAB

LDAB is INTEGER — the leading dimension of the array A. LDAB >= KD+1.

AFB

AFB is REAL array, dimension (LDAFB,N).

  • If FACT = 'F', then AFB 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 of the band matrix A, in the same storage format as A (see AB). If EQUED = 'Y', then AFB is the factored form of the equilibrated matrix A.
  • If FACT = 'N', then AFB 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.
  • If FACT = 'E', then AFB 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 A for the form of the equilibrated matrix).

LDAFB

LDAFB is INTEGER — the leading dimension of the array AFB. LDAFB >= KD+1.

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.

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

B

B is 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

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

X

X is 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

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

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.

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. lowkey humble of it 🫶

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

WORK

WORK is REAL array, dimension (3*N)

IWORK

IWORK is INTEGER array, dimension (N)

INFO

INFO is INTEGER

  • = 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 minor of order i of A is not positive definite, 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 👑

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

🌈 Further Details 🌈

The band storage scheme is illustrated by the following example, when N = 6, KD = 2, and UPLO = 'U':

Two-dimensional storage of the symmetric matrix A:

   a11  a12  a13
        a22  a23  a24
             a33  a34  a35
                  a44  a45  a46
                       a55  a56
   (aij=conjg(aji))         a66

Band storage of the upper triangle of A:

    *    *   a13  a24  a35  a46
    *   a12  a23  a34  a45  a56
   a11  a22  a33  a44  a55  a66

Similarly, if UPLO = 'L' the format of A is as follows:

   a11  a22  a33  a44  a55  a66
   a21  a32  a43  a54  a65   *
   a31  a42  a53  a64   *    *

Array elements marked * are not used by the routine. they're just there for the aesthetic 🙈

🦋 subroutine sppsv (character UPLO, integer N, integer NRHS, real, dimension( * ) AP, real, dimension( ldb, * ) B, integer LDB, integer INFO) 🦋

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

Purpose

SPPSV computes 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.

The Cholesky decomposition is used to factor A 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. The factored form of A is then used to solve the system of equations A * X = B.

Parameters

UPLO

UPLO is CHARACTER*1

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

N

  • N — N is INTEGER. The number of linear equations, aka the order of the matrix A. So like, N >= 0, periodt.
  • NRHS — NRHS is INTEGER. The number of right hand sides, aka the number of columns of the matrix B. NRHS >= 0. 💅
  • AP — AP is 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. The j-th column of A lives in the array 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. 👀

On exit, if INFO = 0, this holds the factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T, in the same storage format as A. Iconic behavior honestly. ✨

  • B — B is REAL array, dimension (LDB,NRHS).

On entry: the N-by-NRHS right hand side matrix B. On exit, if INFO = 0, it's the N-by-NRHS solution matrix X. glow up fr fr 🌸

  • LDB — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
  • INFO — INFO is INTEGER.
    • = 0: successful exit 🎉
    • < 0: if INFO = -i, the i-th argument had an illegal value 😭
    • > 0: 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 has not been computed. 💀

👑 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 ]

like, look at her go 🦋

🫶 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. she's thorough, bestie 💖

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

  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. If the leading i-by-i principal minor is not positive definite, 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.
  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.
  4. 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 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 — AFP is 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 — 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.

  • 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. 🍓
  • B — B is 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 — LDB is INTEGER. The leading dimension of the array B. LDB >= max(1,N).
  • X — X is 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 — LDX is INTEGER. The leading dimension of the array X. LDX >= max(1,N).
  • 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. 🚨
  • 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.
  • 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).
  • WORK — WORK is REAL array, dimension (3*N)
  • IWORK — IWORK is INTEGER array, dimension (N)
  • INFO
  • INFO is an INTEGER — think of it as the ✨vibe report✨ for the whole call:
    • = 0: successful exit 🎉 everything slayed, no notes.
    • < 0: if INFO = -i, the i-th argument had an illegal value 🙈 — like, argument number i is the problem child.
    • > 0: if INFO = i, and i is...
      • <= N: the leading minor of order i of A is not positive definite 😭, so the factorization could not be completed, and the solution has not been computed. RCOND = 0 is returned. (Translation: A said "no thank you" and we got nothing back.)
      • = 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. (Bad vibes, but the answer might still eat. iykyk.)

👑 Author 👑

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

(The group project that actually carried. 💅)

🎀 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 ]

Like — the whole upper triangle gets folded up into one lil linear array. Iconic. 🦋

🌸 subroutine sspsv 🌸

subroutine sspsv (character UPLO, integer N, integer NRHS, real, dimension(
* ) AP, integer, dimension( * ) IPIV, real, dimension( ldb, * ) B,
integer LDB, integer INFO)

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

✨ Purpose ✨

SSPSV computes the solution to a real system of linear equations

   A * X = B,

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

The diagonal pivoting method is used to factor A as

   A = U * D * U**T,  if UPLO = 'U', or
   A = L * D * L**T,  if UPLO = 'L',

where U (or L) is a product of permutation and unit upper (lower) triangular matrices, D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B. Bestie, that's the whole plot. 🫶

🎀 Parameters 🎀

  • 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 matrix B. NRHS >= 0.
  • AP
    • AP is 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. 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. (The receipts are in "Further Details", stay tuned 👀)
    • On exit, the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as a packed triangular matrix in the same storage format as A.
  • IPIV
    • IPIV is INTEGER array, dimension (N)
    • Details of the interchanges and the block structure of D, as determined by SSPTRF. If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  • B
    • B is REAL array, dimension (LDB,NRHS)
    • On entry, the N-by-NRHS right hand side matrix B.
    • On exit, if INFO = 0, the N-by-NRHS solution matrix X. (There it is. The moment. 👑)
  • LDB
    • LDB is INTEGER
    • The leading dimension of the array B. LDB >= max(1,N).
  • INFO
    • INFO is INTEGER
    • = 0: successful exit
    • < 0: if INFO = -i, the i-th argument had an illegal value
    • > 0: if INFO = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed. 😭 (D really said "0" and meant it.)

👑 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 = aji)
               a44

Packed storage of the upper triangle of A:

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

Note the (aij = aji) here — real symmetric, no conjugation, we keep it simple. 💅

💖 subroutine sspsvx 💖

subroutine sspsvx (character FACT, character UPLO, integer N, integer NRHS,
real, dimension( * ) AP, real, dimension( * ) AFP, integer, dimension( *
) IPIV, 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)

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

✨ Purpose ✨

SSPSVX uses the diagonal pivoting factorization A = U*D*U**T or A = L*D*L**T to compute the solution to a real system of linear equations A * X = B, where A is an N-by-N symmetric 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. Like, not just the answer — the receipts about how much you can trust the answer. 🔥

🌈 Description 🌈

The following steps are performed:

  1. If FACT = 'N', the diagonal pivoting method is used to factor A as
   A = U * D * U**T,  if UPLO = 'U', or
   A = L * D * L**T,  if UPLO = 'L',

where U (or L) is a product of permutation and unit upper (lower) triangular matrices and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

  1. If some D(i,i)=0, so that D is exactly singular, 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 ≠ stop, bestie — it keeps going. 🫶)
  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.

🎀 Parameters 🎀

  • FACT
    • FACT is CHARACTER*1
    • Specifies whether or not the factored form of A has been supplied on entry.
    • = 'F': On entry, AFP and IPIV contain the factored form of A. AP, AFP and IPIV will not be modified.
    • = 'N': The matrix A will be 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 REAL array, dimension (N*(N+1)/2)
    • The upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
    • See below for further details.
  • AFP
    • AFP is REAL array, dimension (N*(N+1)/2)
    • If FACT = 'F', then AFP is an input argument and on entry contains the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as a packed triangular matrix in the same storage format as A.
    • If FACT = 'N', then AFP is an output argument and on exit contains the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as a packed triangular matrix in the same storage format as A.
  • IPIV
    • IPIV is INTEGER array, dimension (N)
    • If FACT = 'F', then IPIV is an input argument and on entry contains details of the interchanges and the block structure of D, as determined by SSPTRF.
    • If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal block.
    • If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
    • If FACT = 'N', then IPIV is an output argument and on exit contains details of the interchanges and the block structure of D, as determined by SSPTRF.
  • B
    • B is REAL array, dimension (LDB,NRHS)
    • The N-by-NRHS right hand side matrix B.
  • LDB
  • LDB — LDB is INTEGER ✨ the leading dimension of the array B. LDB >= max(1,N). 💅
  • X — X is REAL array, dimension (LDX,NRHS) 🦋. If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 💖
  • LDX — LDX is INTEGER 🌸 the leading dimension of the array X. LDX >= max(1,N). ✨
  • RCOND — RCOND is REAL 🫶 The estimate of the reciprocal condition number of the matrix A. 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 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. 👑
  • 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). 💀
  • WORK — WORK is REAL array, dimension (3*N) 🌈
  • IWORK — IWORK is INTEGER array, dimension (N) 💜
  • INFO — INFO is INTEGER 🔥
    • = 0: successful exit 🎉
    • < 0: if INFO = -i, the i-th argument had an illegal value 😭
    • > 0: if INFO = i, and i is
      • <= N: D(i,i) is exactly zero. The factorization has been completed but the factor D is exactly singular, so the solution and error bounds could not be computed. RCOND = 0 is returned. 💀
      • = N+1: D 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 👑

  • 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 = aji)
               a44

Packed storage of the upper triangle of A:

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

👑 Author 👑

Generated automatically by Doxygen for LAPACK from the source code.

💖 That's the whole realOTHERsolve(3) page, bestie — LAPACK Version 3.11.0, dated Sun Nov 27 2022. Slay, you made it. 💅✨

📦 where this page comes from