|
|
◆ zgssvx()
| void zgssvx |
( |
char |
equil, |
|
|
char |
col_perm, |
|
|
char |
trans, |
|
|
char |
refine, |
|
|
double |
thresh, |
|
|
char |
sym_mode, |
|
|
int |
n, |
|
|
doublecomplex |
val[], |
|
|
const int |
ptr[], |
|
|
const int |
ind[], |
|
|
int |
format, |
|
|
int |
base, |
|
|
int |
perm_c[], |
|
|
int |
perm_r[], |
|
|
int |
etree[], |
|
|
char * |
equed, |
|
|
double |
r[], |
|
|
double |
c[], |
|
|
int |
nrhs, |
|
|
int |
ldb, |
|
|
doublecomplex |
b[], |
|
|
int |
ldx, |
|
|
doublecomplex |
x[], |
|
|
double * |
rpg, |
|
|
double * |
rcond, |
|
|
double |
ferr[], |
|
|
double |
berr[], |
|
|
int * |
info |
|
) |
| |
Solves the system of linear equations A*X = B, A^T*X = B or A^H*X = B (complex sparse matrix, direct method (SuperLU)) (expert driver)
- Purpose
- This routine solves the system of linear equations A*X = B, A^T*X = B or A^H*X = B, using the LU factorization, where A is a general N x N complex sparse matrix. Error bounds on the solution and a condition estimate are also provided. It performs the following steps.
- If A is stored column-wise (CSC):
1.1. If equil = 'E', scaling factors are computed to equilibrate the system as follows.
If trans = 'N': diag(R)*A*diag(C) * inv(diag(C))*X = diag(R)*B
If trans = 'T' or 'C': (diag(R)*A*diag(C))^T * 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').
1.2. Permute the columns of A, forming A*Pc, where Pc is a permutation matrix that usually preserves sparsity.
1.3. The LU decomposition is used to factor the matrix A (after equilibration if equil = 'E') as Pr*A*Pc = L*U, with Pr determined by partial pivoting.
1.4. Compute the reciprocal pivot growth factor.
1.5. If some U(i, i) = 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 computes error bounds as described below.
1.6. The system of equations is solved for X using the factored form of A.
1.7. If refine != 'N', iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
1.8. 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.
- If A is stored row-wise (CSR), apply the above algorithm to the transpose of A:
2.1. If equil = 'E', scaling factors are computed to equilibrate the system as follows.
If trans = 'N': diag(R)*A*diag(C) * inv(diag(C))*X = diag(R)*B
If trans = 'T' or 'C': (diag(R)*A*diag(C))^T * 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^T is overwritten by diag(R)*A^T*diag(C) and B by diag(R)*B (if trans = 'N') or diag(C)*B (if trans = 'T' or 'C').
2.2. Permute columns of A^T (rows of A), forming A^T*Pc, where Pc is a permutation matrix that usually preserves sparsity.
2.3. The LU decomposition is used to factor the A^T (after equilibration if equil = 'E') as Pr*A^T*Pc = L*U with the permutation Pr determined by partial pivoting.
2.4. Compute the reciprocal pivot growth factor.
2.5. If some U(i,i) = 0, so that U is exactly singular, then the routine returns with info = i. Otherwise, the factored form of A^T 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 computes error bounds as described below.
2.6. The system of equations is solved for X using the factored form of A^T.
2.7. If refine != 'N', iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
2.8. 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] | equil | = 'E': Equilibrates the system (scale A's row and columns to have unit norm).
= 'N': Do not equilibrates the system. |
| [in] | col_perm | Specifies the type of column ordering to reduce fill-in.
= 'A': Approximate minimum degree (AMD) ordering.
= 'N': Natural ordering (No ordering) (Pc = I).
= 'M': Multiple minimum degree (MMD) ordering on A^T*A.
= 'P': Multiple minimum degree (MMD) ordering on A^T + A.
= 'U': User specified ordering. |
| [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] | refine | Specifies whether to perform iterative refinement.
= 'N': No iterative refinement.
= 'S': Perform iterative refinement in single precision.
= 'D': Perform iterative refinement in double precision.
= 'X': Perform iterative refinement in extra precision. |
| [in] | thresh | Specifies the threshold used for a diagonal entry to be an acceptable pivot. (0 <= thresh <= 1) |
| [in] | sym_mode | Specifies whether to use symmetric mode. Symmetric mode gives preference to diagonal pivots, and uses an (A^T + A)-based column permutation algorithm.
= 'N': Do not use symmetric mode.
= 'S': Use symmetric mode. |
| [in] | n | Dimension of the matrix. (n >= 0) (If n = 0, returns without computation) |
| [in] | val[] | Array val[lval] (lval >= nnz)
Values of nonzero elements of the matrix A in CSC/CSR format (where nnz is the number of nonzero elements). |
| [in] | ptr[] | Array ptr[lptr] (lptr >= n + 1)
Column/row pointers of the matrix A in CSC/CSR format. |
| [in] | ind[] | Array ind[lind] (lind >= nnz)
Row/column indices of the matrix A in CSC/CSR format (where nnz is the number of nonzero elements). |
| [in] | format | Sparse matrix format.
= 0: CSR format.
= 1: CSC format. |
| [in] | base | Indexing of ptr[] and ind[].
= 0: Zero-based (C style) indexing: Starting index is 0.
= 1: One-based (Fortran style) indexing: Starting index is 1. |
| [in,out] | perm_c[] | Array perm_c[lperm_c] (lperm_c >= n)
[in] Input permutation vector if option col_perm = 'U'. Otherwise, an output argument.
[out] In the case of CSC format, column permutation vector which defines the permutation matrix Pc. perm_c[i] = j means column i of A is in position j in A*Pc. Similarly, in the case of CSR, column permutation vector which describes permutation of columns of A^T (rows of A).
Input permutation vector may be overwritten by the product of the input perm_c[] and a permutation that postorders the elimination tree of Pc^T*A^T*A*Pc. perm_c[] is not changed if the elimination tree is already in postorder. |
| [out] | perm_r[] | Array perm_r[lperm_r] (lperm_r >= n)
In the case of CSC format, row permutation vector which defines the permutation matrix Pr, and is determined by partial pivoting. perm_r[i] = j means row i of A is in position j in Pr*A. In the case of CSR format, permutation vector which determines permutation of rows of A^T (columns of A) as described above. |
| [out] | etree[] | Array etree[letree] (letree >= n)
Elimination tree of Pc^T*A^T*A*Pc.
Note: etree[] is a vector of parent pointers for a forest whose vertices are the integers 0 to n - 1; Etree(root) = n. |
| [in,out] | equed | Specifies the form of equilibration that was done.
= 'N': No equilibration.
= 'R': Row equilibration, i.e., A was premultiplied by diag(R).
= 'C': Column equilibration, i.e., A was postmultiplied by diag(C).
= 'B': Both row and column equilibration, i.e., A was replaced by diag(R)*A*diag(C). |
| [out] | r[] | Array r[lr] (lr >= n)
The row scale factors for A or A^T.
If equed = 'R' or 'B', A (if CSC) or A^T (if CSR) is multiplied on the left by diag(R).
If equed = 'N' or 'C', r[] is not accessed. |
| [out] | c[] | Array c[lc] (lc >= n)
The column scale factors for A or A^T.
If equed = 'C' or 'B', A (if CSC) or A^T (if CSR) is multiplied on the right by diag(C).
If equed = 'N' or 'R', c[] is not accessed. |
| [in] | nrhs | Number of right hand sides, i.e., number of columns of the matrix B. (nrhs >= 0) (If nrhs = 0, returns without computation) |
| [in] | ldb | Leading dimension of the two dimensional array b[][]. (ldb >= max(0, n)) |
| [in,out] | b[][] | Array b[lb][ldb] (lb >= nrhs)
[in] The right hand side matrix B.
[out] If equed = 'N' on exit, b[][] is not modified.
Otherwise, b[][] is modified as follows.
If A is in CSC format:
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.
If A is in CSR format:
If trans = 'N' and equed = 'C' or 'B', b[][] is overwritten by diag(C)*B.
If trans = 'T' or 'C' and equed = 'R' or 'B', b[][] is overwritten by diag(R)*B. |
| [in] | ldx | Leading dimension of the two dimensional array x[][]. (ldx >= max(0, n)) |
| [out] | x[][] | Array x[lx][ldx] (lx >= nrhs)
Solution matrix X.
If info = 0 or info = n + 1, X contains the solution matrix to the original system of equations (before equilibration).
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] | rpg | The reciprocal pivot growth factor max_j(norm(A_j)/norm(U_j)). The infinity norm is used.
If rpg is much less than 1, the stability of the LU factorization could be poor. |
| [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 - 1) (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] 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. If refine = 'N', ferr[j] = 1. |
| [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). If refine = 'N', berr[j] = 1. |
| [out] | info | = 0: Successful exit.
< 0: The (-info)-th argument is invalid.
= -10000: Unrecoverble error occured in SuperLU routine.
= i > 0 and <= n: U(i, i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.
= 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.
= i > n + 1: Number of bytes allocated when memory allocation failure occurred. i - n - 1 is the number of bytes allocated |
- Reference
- SuperLU
|
1.9.7