|
|
◆ zgbsvx()
| void zgbsvx |
( |
char |
fact, |
|
|
char |
trans, |
|
|
int |
n, |
|
|
int |
kl, |
|
|
int |
ku, |
|
|
int |
nrhs, |
|
|
int |
ldab, |
|
|
doublecomplex |
ab[], |
|
|
int |
ldafb, |
|
|
doublecomplex |
afb[], |
|
|
int |
ipiv[], |
|
|
char * |
equed, |
|
|
double |
r[], |
|
|
double |
c[], |
|
|
int |
ldb, |
|
|
doublecomplex |
b[], |
|
|
int |
ldx, |
|
|
doublecomplex |
x[], |
|
|
double * |
rcond, |
|
|
double |
ferr[], |
|
|
double |
berr[], |
|
|
doublecomplex |
work[], |
|
|
double |
rwork[], |
|
|
int * |
info |
|
) |
| |
(Expert driver) Solution to system of linear equations AX = B for a complex band matrix
- Purpose
- This routine uses the LU factorization to computes the solution to a complex system of linear equations
A * X = B, A^T * X = B or A^H * X = B
where A is a band matrix of order n with kl sub-diagonals and ku super-diagonals, and X and B are n x nrhs matrices.
Error bounds on the solution and a condition estimate are also provided.
- Description
- The following steps are performed by this subroutine:
- If fact = 'E', real scaling factors are computed to equilibrate the system:
trans = 'N': diag(R)*A*diag(C)*inv(diag(C))*X = diag(R)*B
trans = 'T': (diag(R)*A*diag(C))^T *inv(diag(R))*X = diag(C)*B
trans = 'C': (diag(R)*A*diag(C))^H *inv(diag(R))*X = diag(C)*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(R)*A*diag(C) and B by diag(R)*B (if trans = 'N') or diag(C)*B (if trans = 'T' or 'C').
- If fact = 'N' or 'E', the LU decomposition is used to factor the matrix A (after equilibration if fact = 'E') as where L is a product of permutation and unit lower triangular matrices with kl sub-diagonals, and U is upper triangular with kl+ku super-diagonals.
- If i-th diagonal element of U = 0, so that U 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.
- 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(C) (if trans = 'N') or diag(R) (if trans = 'T' or 'C') so that it solves the original system before equilibration.
- Parameters
-
| [in] | fact | 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': afb[][] and ipiv[] contain the factored form of A. If equed is not 'N', the matrix A has been equilibrated with scaling factors given by r[] and c[]. ab[][], afb[][] and ipiv[] are not 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. |
| [in] | trans | Specifies the form of the system of equations:
= 'N': A * X = B. (no transpose)
= 'T': A^T * X = B. (transpose)
= 'C': A^H * X = B. (conjugate transpose) |
| [in] | n | Number of linear equations, i.e., order of the matrix A. (n >= 0) (If n = 0, returns without computation) |
| [in] | kl | The number of sub-diagonals within the band of A. (kl >= 0) |
| [in] | ku | The number of super-diagonals within the band of A. (ku >= 0) |
| [in] | nrhs | Number of right hand sides, i.e., number of columns of the matrices B and X. (nrhs >= 0) (If nrhs = 0, returns without computation) |
| [in] | ldab | Leading dimension of the two dimensional array ab[][]. (ldab >= kl + ku + 1) |
| [in,out] | ab[][] | Array ab[lab][ldab] (lab >= n)
[in] Matrix A in band matrix form, in rows 1 to kl+ku+1. If fact = 'F' and equed != 'N', then A must have been equilibrated by the scaling factors in r[] and/or c[].
[out] Not modified if fact = 'F' or 'N', or if fact = 'E' and equed = 'N' on exit.
If fact = 'E' and equed != 'N' on exit, A is scaled as follows:
equed = 'R': A := diag(R)*A,
equed = 'C': A := A*diag(C), or
equed = 'B': A := diag(R)*A*diag(C). |
| [in] | ldafb | Leading dimension of the two dimensional array afb[][]. (ldafb >= 2kl + ku + 1) |
| [in,out] | afb[][] | Array afb[lafb][ldafb] (lafb >= n)
[in] If fact = 'F', details of the LU factorization of the band matrix A, as computed by zgbtrf. U is stored as an upper triangular band matrix with kl+ku super-diagonals in rows 1 to kl+ku+1, and the multipliers used during the factorization are stored in rows kl+ku+2 to 2*kl+ku+1. If equed != 'N', then afb[][] is the factored form of the equilibrated matrix A.
[out] If fact = 'N', returns details of the LU factorization of A.
If fact = 'E', returns details of the LU factorization of the equilibrated matrix A (see the description of ab[][] for the form of the equilibrated matrix). |
| [out] | ipiv[] | Array ipiv[lipiv] (lipiv >= n)
[in] If fact = 'F', the pivot indices from the factorization A = L*U as computed by zgbtrf; row i of the matrix was interchanged with row ipiv[i-1].
[out] If fact = 'N', the pivot indices from the factorization A = L*U of the original matrix A.
If fact = 'E', the pivot indices from the factorization A = L*U of the equilibrated matrix A. |
| [in,out] | equed | Specifies the form of equilibration that was done.
= 'N': No equilibration.
= 'R': Row equilibration, i.e., A has been premultiplied by diag(R).
= 'C': Column equilibration, i.e., A has been postmultiplied by diag(C).
= 'B': Both row and column equilibration, i.e., A has been replaced by diag(R)*A*diag(C).
[in] If fact = 'F', the form of equilibration of the supplied matrix in afb[][].
[out] If fact != 'F', the form of equilibration that was done ('N', 'R', 'C' or 'B'). If fact = 'N', equed is always 'N'. |
| [in,out] | r[] | Array r[lr] (lr >= n)
The row scale factors diag(R) for A. If equed = 'R' or 'B', A is multiplied on the left by diag(R); If equed = 'N' or 'C', r[] is not accessed.
[in] If fact = 'F', the row scale factors for the supplied matrix in afb[][]. (Each element must be > 0)
[out] If fact != 'F', the resulted row scale factors for A. |
| [in,out] | c[] | Array c[lc] (lc >= n)
The column scale factors diag(C) for A. If equed = 'C' or 'B', A is multiplied on the right by diag(C); If equed = 'N' or 'R', c[] is not accessed.
[in] If fact = 'F', the column scale factors for the supplied matrix in afb[][]. (Each element must be > 0)
[out] If fact != 'F', the resulted column scale factors for A. |
| [in] | ldb | Leading dimension of the two dimensional array b[][]. (ldb >= max(1, n)) |
| [in,out] | b[][] | Array b[lb][ldb] (lb >= nrhs)
[in] n x nrhs right hand side matrix B.
[out] If equed = 'N', b[][] is not modified.
If trans = 'N' and equed = 'R' or 'B', b[][] is overwritten by diag(R)*B.
If trans = 'T' or 'C' and equed = 'C' or 'B', b[][] is overwritten by diag(C)*B. |
| [in] | ldx | Leading dimension of the two dimensional array x[]. (ldx >= max(1, n)) |
| [out] | x[][] | Array x[lx][ldx] (lx >= nrhs)
If info = 0 or info = n+1, the n x nrhs solution matrix X to the original system of equations. Note that A and B are modified on exit if equed != 'N', and the solution to the equilibrated system is inv(diag(C))*X if trans = 'N' and equed = 'C' or 'B', or inv(diag(R))*X if trans = 'T' or 'C' and equed = 'R' or 'B'. |
| [out] | rcond | 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. |
| [out] | ferr[] | Array ferr[lferr] (lferr >= 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-1] 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. |
| [out] | berr[] | Array berr[lberr] (lberr >= 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). |
| [out] | work[] | Array work[lwork] (lwork >= 2*n)
Work array. |
| [out] | rwork[] | Array rwork[lrwork] (lrwork >= n)
Work array.
On exit, rwork[0] contains the reciprocal pivot growth factor norm(A)/norm(U). The "max absolute element" norm is used. If rwork[0] is much less than 1, then the stability of the LU factorization of the (equilibrated) matrix A could be poor. This also means that the solution X, condition estimator rcond, and forward error bound ferr could be unreliable. If factorization fails with 0 < info <= n, then rwork[0] contains the reciprocal pivot growth factor for the leading info columns of A. |
| [out] | info | = 0: Successful exit
= -1: The argument fact had an illegal value (fact != 'F', 'N' nor 'E')
= -2: The argument trans had an illegal value (trans != 'N', 'T' nor 'C')
= -3: The argument n had an illegal value (n < 0)
= -4: The argument kl had an illegal value (kl < 0)
= -5: The argument ku had an illegal value (ku < 0)
= -6: The argument nrhs had an illegal value (nrhs < 0)
= -7: The argument ldab had an illegal value (ldab < kl+ku+1)
= -9: The argument ldafb had an illegal value (ldafb < 2kl+ku+1)
= -12: The argument equed had an illegal value (fact = 'F' and equed != 'N', 'R', 'C' nor 'B')
= -13: The argument r had an illegal value (r[i] <= 0 when fact = 'F' and equed = 'R' or 'B')
= -14: The argument c had an illegal value (c[i] <= 0 when fact = 'F' and equed = 'C' or 'B')
= -15: The argument ldb had an illegal value (ldb < max(1, n))
= -17: The argument ldx had an illegal value (ldx < max(1, n))
= i (0 < i <= n): The i-th diagonal element of the factor U is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution and error bounds could not be 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. |
- Reference
- LAPACK
|