pbsvx(3)
okay so this is pbsvx(3) from the Library Functions Manual โ yes, the same pbsvx(3) up in the header โ and she is here to factor and solve like an expert. ๐
โจ
โจ NAME โจ
pbsvx โ pbsvx: factor and solve, expert
๐ SYNOPSIS ๐
๐ซถ Functions ๐ซถ
cpbsvx โ CPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine cpbsvx (fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb, equed,
s, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
dpbsvx โ DPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine dpbsvx (fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb, equed,
s, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
spbsvx โ SPBSVX 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)
zpbsvx โ ZPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zpbsvx (fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb, equed,
s, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
๐ฆ Detailed Description ๐ฆ
No extra detail here in this part, bestie โ the page just moves straight into the function documentation. ๐
๐ Function Documentation ๐
๐งฎ cpbsvx โ complex expert mode
subroutine cpbsvx (character fact, character uplo, integer n, integer kd,
integer nrhs, complex, dimension( ldab, * ) ab, integer ldab, complex,
dimension( ldafb, * ) afb, integer ldafb, 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)
CPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
๐ฏ Purpose
okay so CPBSVX 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 band matrix and X and B are N-by-NRHS matrices. Error bounds on the solution and a condition estimate are also provided. ๐
๐ฆ Description
The following steps are performed:
- 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.
- If
FACT = 'N'or'E', the Cholesky decomposition is used to factor the matrix A (after equilibration ifFACT = 'E') as
A = U**H * U, if UPLO = 'U', or
A = L * L**H, if UPLO = 'L',
where U is an upper triangular band matrix, and L is a lower triangular band matrix.
- 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+1is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. ๐ warning behavior, bestie โ not a full stop. - The system of equations is solved for X using the factored form of A.
- Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
- If equilibration was used, the matrix X is premultiplied by
diag(S)so that it solves the original system before equilibration.
๐ 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,AFBcontains the factored form of A. IfEQUED = 'Y', the matrix A has been equilibrated with scaling factors given byS.ABandAFBwill not be modified.= 'N': The matrix A will be copied toAFBand factored.= 'E': The matrix A will be equilibrated if necessary, then copied toAFBand 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.KDโINTEGER. The number of superdiagonals of the matrix A ifUPLO = 'U', or the number of subdiagonals ifUPLO = 'L'.KD >= 0.NRHSโINTEGER. The number of right-hand sides, i.e., the number of columns of the matrices B and X.NRHS >= 0.ABโCOMPLEX array, dimension (LDAB,N). On entry, the upper or lower triangle of the Hermitian band matrix A, stored in the firstKD+1rows of the array, except ifFACT = 'F'andEQUED = 'Y', then A must contain the equilibrated matrixdiag(S)*A*diag(S). The j-th column of A is stored in the j-th column of the arrayABas follows: ifUPLO = 'U',AB(KD+1+i-j,j) = A(i,j)formax(1,j-KD)<=i<=j; ifUPLO = 'L',AB(1+i-j,j) = A(i,j)forj<=i<=min(N,j+KD). See below for further details. On exit, ifFACT = 'E'andEQUED = 'Y', A is overwritten bydiag(S)*A*diag(S).LDABโINTEGER. The leading dimension of the array A.LDAB >= KD+1.AFBโCOMPLEX array, dimension (LDAFB,N).- If
FACT = 'F', thenAFBis an input argument and on entry contains the triangular factor U or L from the Cholesky factorizationA = U**H*UorA = L*L**Hof the band matrix A, in the same storage format as A (seeAB). IfEQUED = 'Y', thenAFBis the factored form of the equilibrated matrix A. - If
FACT = 'N', thenAFBis an output argument and on exit returns the triangular factor U or L from the Cholesky factorizationA = U**H*UorA = L*L**H. - If
FACT = 'E', thenAFBis an output argument and on exit returns the triangular factor U or L from the Cholesky factorizationA = U**H*UorA = L*L**Hof the equilibrated matrix A (see the description of A for the form of the equilibrated matrix).
- If
LDAFBโINTEGER. The leading dimension of the array AFB.LDAFB >= KD+1.EQUEDโCHARACTER*1. Specifies the form of equilibration that was done.= 'N': No equilibration (always true ifFACT = 'N').= 'Y': Equilibration was done, i.e., A has been replaced bydiag(S) * A * diag(S).EQUEDis an input argument ifFACT = 'F'; otherwise, it is an output argument.
SโREAL array, dimension (N). The scale factors for A; not accessed ifEQUED = 'N'.Sis an input argument ifFACT = 'F'; otherwise,Sis an output argument. IfFACT = 'F'andEQUED = 'Y', each element ofSmust be positive.BโCOMPLEX array, dimension (LDB,NRHS). On entry, the N-by-NRHS right hand side matrix B. On exit, ifEQUED = 'N', B is not modified; ifEQUED = 'Y', B is overwritten bydiag(S) * B.LDBโINTEGER. The leading dimension of the array B.LDB >= max(1,N).XโCOMPLEX array, dimension (LDX,NRHS). IfINFO = 0orINFO = N+1, the N-by-NRHS solution matrix X to the original system of equations. Note that ifEQUED = 'Y', A and B are modified on exit, and the solution to the equilibrated system isinv(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). IfRCONDis less than the machine precision (in particular, ifRCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code ofINFO > 0.FERRโREAL array, dimension (NRHS). The estimated forward error bound for each solution vectorX(j)(the j-th column of the solution matrix X). IfXTRUEis the true solution corresponding toX(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 inX(j). The estimate is as reliable as the estimate forRCOND, and is almost always a slight overestimate of the true error.BERRโREAL array, dimension (NRHS). The componentwise relative backward error of each solution vectorX(j)(i.e., the smallest relative change in any element of A or B that makesX(j)an exact solution).WORKโCOMPLEX array, dimension (2*N).RWORKโREAL array, dimension (N).INFOโ (the original page cuts off right here mid-parameter, bestie โINFO's description continues in the next part ๐ )
INFO is INTEGER โ okay bestie, here's the return-code tea you actually need:
= 0: successful exit< 0: ifINFO = -i, thei-th argument had an illegal value> 0: ifINFO = i, andiis<= N: the leading principal minor of orderiofAis not positive, so the factorization could not be completed, and the solution has not been computed.RCOND = 0is returned.= N+1:Uis nonsingular, butRCONDis 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 ofRCONDwould 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 Hermitian 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.
Definition at line 339 of file cpbsvx.f.
๐ง subroutine dpbsvx ๐ง
subroutine dpbsvx (character fact, character uplo, integer n, integer kd,
integer nrhs, double precision, dimension( ldab, * ) ab, integer ldab,
double precision, dimension( ldafb, * ) afb, integer ldafb, 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)
DPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
โญ Purpose
DPBSVX 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. It's giving โsolve it and also hand you the receiptsโ energy, fr fr.
๐ซ Description
The following steps are performed:
- 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 matrixA, but if equilibration is used,Ais overwritten bydiag(S)*A*diag(S)andBbydiag(S)*B. - If
FACT = 'N'or'E', the Cholesky decomposition is used to factor the matrixA(after equilibration ifFACT = 'E') asA = U**T * U, ifUPLO = 'U', orA = L * L**T, ifUPLO = 'L', whereUis an upper triangular band matrix, andLis a lower triangular band matrix. - If the leading principal minor of order
iis not positive, then the routine returns withINFO = i. Otherwise, the factored form ofAis used to estimate the condition number of the matrixA. If the reciprocal of the condition number is less than machine precision,INFO = N+1is returned as a warning, but the routine still goes on to solve forXand compute error bounds as described below. - The system of equations is solved for
Xusing the factored form ofA. - Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
- If equilibration was used, the matrix
Xis premultiplied bydiag(S)so that it solves the original system before equilibration.
๐ Parameters
FACTโFACTisCHARACTER*1. Specifies whether or not the factored form of the matrixAis supplied on entry, and if not, whether the matrixAshould be equilibrated before it is factored.= 'F': On entry,AFBcontains the factored form ofA. IfEQUED = 'Y', the matrixAhas been equilibrated with scaling factors given byS.ABandAFBwill not be modified.= 'N': The matrixAwill be copied toAFBand factored.= 'E': The matrixAwill be equilibrated if necessary, then copied toAFBand factored.
UPLOโUPLOisCHARACTER*1.= 'U': Upper triangle ofAis stored;= 'L': Lower triangle ofAis stored.
NโNisINTEGER. The number of linear equations, i.e., the order of the matrixA.N >= 0.KDโKDisINTEGER. The number of superdiagonals of the matrixAifUPLO = 'U', or the number of subdiagonals ifUPLO = 'L'.KD >= 0.NRHSโNRHSisINTEGER. The number of right-hand sides, i.e., the number of columns of the matricesBandX.NRHS >= 0.ABโABisDOUBLE PRECISION array, dimension (LDAB,N). On entry, the upper or lower triangle of the symmetric band matrixA, stored in the firstKD+1rows of the array, except ifFACT = 'F'andEQUED = 'Y', thenAmust contain the equilibrated matrixdiag(S)*A*diag(S). Thej-th column ofAis stored in thej-th column of the arrayABas follows: ifUPLO = 'U',AB(KD+1+i-j,j) = A(i,j)formax(1,j-KD)<=i<=j; ifUPLO = 'L',AB(1+i-j,j) = A(i,j)forj<=i<=min(N,j+KD). See below for further details.- On exit, if
FACT = 'E'andEQUED = 'Y',Ais overwritten bydiag(S)*A*diag(S).
- On exit, if
LDABโLDABisINTEGER. The leading dimension of the arrayA.LDAB >= KD+1.AFBโAFBisDOUBLE PRECISION array, dimension (LDAFB,N). IfFACT = 'F', thenAFBis an input argument and on entry contains the triangular factorUorLfrom the Cholesky factorizationA = U**T*UorA = L*L**Tof the band matrixA, in the same storage format asA(seeAB). IfEQUED = 'Y', thenAFBis the factored form of the equilibrated matrixA.- If
FACT = 'N', thenAFBis an output argument and on exit returns the triangular factorUorLfrom the Cholesky factorizationA = U**T*UorA = L*L**T. - If
FACT = 'E', thenAFBis an output argument and on exit returns the triangular factorUorLfrom the Cholesky factorizationA = U**T*UorA = L*L**Tof the equilibrated matrixA(see the description ofAfor the form of the equilibrated matrix).
- If
LDAFBโLDAFBisINTEGER. The leading dimension of the arrayAFB.LDAFB >= KD+1.EQUEDโEQUEDisCHARACTER*1. Specifies the form of equilibration that was done.= 'N': No equilibration (always true ifFACT = 'N').= 'Y': Equilibration was done, i.e.,Ahas been replaced bydiag(S) * A * diag(S).EQUEDis an input argument ifFACT = 'F'; otherwise, it is an output argument.
SโSisDOUBLE PRECISION array, dimension (N). The scale factors forA; not accessed ifEQUED = 'N'.Sis an input argument ifFACT = 'F'; otherwise,Sis an output argument. IfFACT = 'F'andEQUED = 'Y', each element ofSmust be positive.BโBisDOUBLE PRECISION array, dimension (LDB,NRHS). On entry, theN-by-NRHSright hand side matrixB. On exit, ifEQUED = 'N',Bis not modified; ifEQUED = 'Y',Bis overwritten bydiag(S) * B.LDBโLDBisINTEGER. The leading dimension of the arrayB.LDB >= max(1,N).XโXisDOUBLE PRECISION array, dimension (LDX,NRHS). IfINFO = 0orINFO = N+1, theN-by-NRHSsolution matrixXto the original system of equations. Note that ifEQUED = 'Y',AandBare modified on exit, and the solution to the equilibrated system isinv(diag(S))*X.LDXโLDXisINTEGER. The leading dimension of the arrayX.LDX >= max(1,N).RCONDโRCONDisDOUBLE PRECISION. The estimate of the reciprocal condition number of the matrixAafter equilibration (if done). IfRCONDis less than the machine precision (in particular, ifRCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code ofINFO > 0.
FERR
FERRโFERR is DOUBLE PRECISION array, dimension (NRHS). This is the estimated forward error bound for each solution vectorX(j)(the j-th column of the solution matrix X). IfXTRUEis the true solution corresponding toX(j), thenFERR(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 inX(j). The estimate is as reliable as the estimate forRCOND, and is almost always a slight overestimate of the true error. ๐BERRโBERR is DOUBLE PRECISION array, dimension (NRHS). This is the componentwise relative backward error of each solution vectorX(j)(i.e., the smallest relative change in any element of A or B that makesX(j)an exact solution). โจWORKโWORK is DOUBLE PRECISION array, dimension (3*N). ๐IWORKโIWORK is INTEGER array, dimension (N). ๐ฆINFOโINFO is INTEGER.= 0: successful exit ๐ธ< 0: ifINFO = -i, the i-th argument had an illegal value> 0: ifINFO = 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 = 0is returned.= N+1: U is nonsingular, butRCONDis 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 ofRCONDwould 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.
Definition at line 340 of file dpbsvx.f.
๐ 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 โจ
Okay so 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. Obsessed, honestly.
๐ธ Description ๐ธ
The following steps are performed:
- 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.
- If
FACT = 'N'or'E', the Cholesky decomposition is used to factor the matrix A (after equilibration ifFACT = '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.
- 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+1is returned as a warning, but the routine still goes on to solve for X and compute error bounds as described below. - The system of equations is solved for X using the factored form of A.
- Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
- 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,AFBcontains the factored form of A. IfEQUED = 'Y', the matrix A has been equilibrated with scaling factors given byS.ABandAFBwill not be modified.= 'N': The matrix A will be copied toAFBand factored.= 'E': The matrix A will be equilibrated if necessary, then copied toAFBand 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.KDโKD is INTEGER. The number of superdiagonals of the matrix A ifUPLO = 'U', or the number of subdiagonals ifUPLO = '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 firstKD+1rows of the array, except ifFACT = 'F'andEQUED = 'Y', then A must contain the equilibrated matrixdiag(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 factorizationA = U**T*UorA = L*L**Tof the band matrix A, in the same storage format as A (seeAB). IfEQUED = '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 factorizationA = U**T*UorA = 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 factorizationA = U**T*UorA = L*L**Tof the equilibrated matrix A (see the description of A for the form of the equilibrated matrix).
- If
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 ifFACT = 'N').= 'Y': Equilibration was done, i.e., A has been replaced bydiag(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 ifEQUED = 'N'.Sis an input argument ifFACT = 'F'; otherwise,Sis an output argument. IfFACT = 'F'andEQUED = '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, ifEQUED = 'N', B is not modified; ifEQUED = 'Y', B is overwritten bydiag(S) * B.LDBโLDB is INTEGER. The leading dimension of the array B.LDB >= max(1,N).XXโX is REAL array, dimension (LDX,NRHS)โจ so like, ifINFO = 0orINFO = N+1, this is the N-by-NRHS solution matrix X to the original system of equations. Bestie, read this part twice: ifEQUED = 'Y', thenAandBare modified on exit, and the solution to the equilibrated system isinv(diag(S))*X. That's the whole vibe check right there ๐LDXโLDX is INTEGER, the leading dimension of the arrayX. Constraint, no cap:LDX >= max(1,N).RCONDโRCOND is REAL, the estimate of the reciprocal condition number of the matrixAafter equilibration (if equilibration happened, that is). IfRCONDis less than the machine precision โ and especially ifRCOND = 0โ the matrix is singular to working precision ๐ญ and that whole situation gets flagged by a return code ofINFO > 0.FERRโFERR is REAL array, dimension (NRHS), the estimated forward error bound for each solution vectorX(j)(that's the j-th column of the solution matrixX). Okay so: ifXTRUEis the true solution corresponding toX(j), thenFERR(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 inX(j). This estimate is exactly as reliable as the estimate forRCOND, and it's almost always a slight overestimate of the true error. Slight. Overestimate. Every time. ๐ซถBERRโBERR is REAL array, dimension (NRHS), the componentwise relative backward error of each solution vectorX(j)(i.e., the smallest relative change in any element ofAorBthat makesX(j)an exact solution). Iconic little measurement honestly.WORKโWORK is REAL array, dimension (3*N). That's it, that's the workspace ๐IWORKโIWORK is INTEGER array, dimension (N). Integer workspace, she's doing her thing.INFOโINFO is INTEGER, and here's the breakdown bestie, memorize it:= 0: successful exit ๐ slay< 0: ifINFO = -i, the i-th argument had an illegal value ๐ (you passed something bad, not my fault)> 0: ifINFO = i, and i is:<= N: the leading principal minor of order i ofAis not positive, so the factorization could not be completed, and the solution has not been computed.RCOND = 0is returned. That's a wrap, no solution for you ๐= N+1:Uis nonsingular, butRCONDis 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 ofRCONDwould suggest. Not me being generous with the answers โจ
๐ 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'. Get a load of this, it's lowkey gorgeous:
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. Decor. Background characters ๐
Definition at line 340 of file spbsvx.f.
๐ subroutine zpbsvx ๐
Like, real talk, this is the complex*16 glow-up of everything above. Same energy, double the precision, main character energy fr fr ๐ฅ
subroutine zpbsvx (character fact, character uplo, integer n, integer kd,
integer nrhs, complex*16, dimension( ldab, * ) ab, integer ldab,
complex*16, dimension( ldafb, * ) afb, integer ldafb, 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)
๐ ZPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices ๐
โจ Purpose โจ
So ZPBSVX 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 band matrix and X and B are N-by-NRHS matrices. ๐
Error bounds on the solution and a condition estimate are also provided. Yes, she does the extra credit too, no cap.
๐ Description ๐
The following steps are performed, in this exact order, bestie:
- 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. โ ๏ธ overwritten, as in, gone, as in, that's the trade.
- If
FACT = 'N'or'E', the Cholesky decomposition is used to factor the matrixA(after equilibration ifFACT = 'E') as
A = U**H * U, if UPLO = 'U', or
A = L * L**H, if UPLO = 'L',
where U is an upper triangular band matrix, and L is a lower triangular band matrix. ๐
- If the leading principal minor of order i is not positive, then the routine returns with
INFO = i. Otherwise, the factored form ofAis used to estimate the condition number of the matrixA. If the reciprocal of the condition number is less than machine precision,INFO = N+1is returned as a warning, but the routine still goes on to solve forXand compute error bounds as described below. She's not gonna leave you hanging just because of a little singularity scare ๐ซถ - The system of equations is solved for
Xusing the factored form ofA. - Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it. Glow-up phase โจ
- If equilibration was used, the matrix
Xis premultiplied bydiag(S)so that it solves the original system before equilibration. Full circle, periodt.
๐ Parameters ๐
FACTโFACT is CHARACTER*1. Specifies whether or not the factored form of the matrixAis supplied on entry, and if not, whether the matrixAshould be equilibrated before it is factored.= 'F': On entry,AFBcontains the factored form ofA. IfEQUED = 'Y', the matrixAhas been equilibrated with scaling factors given byS.ABandAFBwill not be modified. Respect the input, bestie ๐= 'N': The matrixAwill be copied toAFBand factored.= 'E': The matrixAwill be equilibrated if necessary, then copied toAFBand factored.
UPLOโUPLO is CHARACTER*1, aka which half of the triangle you're serving:= 'U': Upper triangle ofAis stored;= 'L': Lower triangle ofAis stored.
NโN is INTEGER, the number of linear equations, i.e., the order of the matrixA. Constraint:N >= 0.KDโKD is INTEGER, the number of superdiagonals of the matrixAifUPLO = 'U', or the number of subdiagonals ifUPLO = 'L'. Constraint:KD >= 0.NRHSโNRHS is INTEGER, the number of right-hand sides, i.e., the number of columns of the matricesBandX. Constraint:NRHS >= 0.ABโAB is COMPLEX*16 array, dimension (LDAB,N). Okay so on entry, it's the upper or lower triangle of the Hermitian band matrixA, stored in the firstKD+1rows of the array โ except ifFACT = 'F'andEQUED = 'Y', thenAmust contain the equilibrated matrixdiag(S)*A*diag(S). The j-th column ofAis stored in the j-th column of the arrayABas 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). Yeah. Your input's gone. That's the whole deal.
LDABโLDAB is INTEGER, the leading dimension of the arrayA. Constraint:LDAB >= KD+1.AFBโAFB is COMPLEX*16 array, dimension (LDAFB,N). Depends entirely onFACT, so pay attention:- If
FACT = 'F', thenAFBis an input argument and on entry contains the triangular factorUorLfrom the Cholesky factorizationA = U**H *UorA = L*L**Hof the band matrixA, in the same storage format asA(seeAB). IfEQUED = 'Y', thenAFBis the factored form of the equilibrated matrixA. - If
FACT = 'N', thenAFBis an output argument and on exit returns the triangular factorUorLfrom the Cholesky factorizationA = U**H *UorA = L*L**H. - If
FACT = 'E', thenAFBis an output argument and on exit returns the triangular factorUorLfrom the Cholesky factorizationA = U**H *UorA = L*L**Hof the equilibrated matrixA(see the description ofAfor the form of the equilibrated matrix).
- If
LDAFBโLDAFB is INTEGER, the leading dimension of the arrayAFB. Constraint:LDAFB >= KD+1.EQUEDEQUEDโ ok so this oneโs aCHARACTER*1, and it specifies the form of equilibration that was done. vibe check:= 'N': no equilibration (always true ifFACT = 'N').= 'Y': equilibration was done, i.e.,Ahas been replaced bydiag(S) * A * diag(S).
EQUED is an input argument if FACT = 'F'; otherwise, it is an output argument. ๐
SโSisDOUBLE PRECISIONarray, dimension(N). The scale factors forA; not accessed ifEQUED = 'N'.Sis an input argument ifFACT = 'F'; otherwise,Sis an output argument. IfFACT = 'F'andEQUED = 'Y', each element ofSmust be positive. no cap, positivity is required.BโBisCOMPLEX*16array, dimension(LDB,NRHS). On entry, theN-by-NRHSright hand side matrixB. On exit, ifEQUED = 'N',Bis not modified; ifEQUED = 'Y',Bis overwritten bydiag(S) * B. so watch out, bestie โ it might get rewritten.LDBโLDBisINTEGER. The leading dimension of the arrayB.LDB >= max(1,N). โจXโXisCOMPLEX*16array, dimension(LDX,NRHS). IfINFO = 0orINFO = N+1, theN-by-NRHSsolution matrixXto the original system of equations. Note that ifEQUED = 'Y',AandBare modified on exit, and the solution to the equilibrated system isinv(diag(S))*X. iconic little transformation, honestly.LDXโLDXisINTEGER. The leading dimension of the arrayX.LDX >= max(1,N). ๐ซถRCONDโRCONDisDOUBLE PRECISION. The estimate of the reciprocal condition number of the matrixAafter equilibration (if done). IfRCONDis less than the machine precision (in particular, ifRCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code ofINFO > 0. translating: tinyRCOND= red flag city.FERRโFERRisDOUBLE PRECISIONarray, dimension(NRHS). The estimated forward error bound for each solution vectorX(j)(thej-th column of the solution matrixX). IfXTRUEis the true solution corresponding toX(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 inX(j). The estimate is as reliable as the estimate forRCOND, and is almost always a slight overestimate of the true error. itโs giving โclose enough, and honest about it.โBERRโBERRisDOUBLE PRECISIONarray, dimension(NRHS). The componentwise relative backward error of each solution vectorX(j)(i.e., the smallest relative change in any element ofAorBthat makesX(j)an exact solution). ๐WORKโWORKisCOMPLEX*16array, dimension(2*N). scratch space bestie, donโt touch.RWORKโRWORKisDOUBLE PRECISIONarray, dimension(N). ๐ธINFOโINFOisINTEGER, aka the final vibe report.= 0: successful exit.< 0: ifINFO = -i, thei-th argument had an illegal value.> 0: ifINFO = i, andiis:<= N: the leading principal minor of orderiofAis not positive, so the factorization could not be completed, and the solution has not been computed.RCOND = 0is returned.= N+1:Uis nonsingular, butRCONDis 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 ofRCONDwould suggest. plot twist: it still tries, and sometimes it slays.
๐ Further Details ๐
okay so the band storage scheme is illustrated by the following example, when N = 6, KD = 2, and UPLO = 'U':
Two-dimensional storage of the Hermitian 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. leave them alone, theyโre just vibing.
๐ Author ๐
- Univ. of Tennessee
- Univ. of California Berkeley
- Univ. of Colorado Denver
- NAG Ltd.
absolute legends, fr fr.
Definition at line 339 of file zpbsvx.f.
๐ Author (Doxygenโs version) ๐
Generated automatically by Doxygen for LAPACK from the source code.
okay bestie, thatโs the whole pbsvx(3) page โ source: LAPACK, version 3.12.0, and the pbsvx(3) footer tag. you ate, no crumbs ๐
โจ๐
๐ฆ where this page comes from
- arch rolling: lapack-doc 3.12.1-2 ยท custom
- ๐ girlified by deepseek-ai/DeepSeek-V4.1-Flash-1017 on : a language model's rewrite of the original, under its terms