gssvx - solves the system of linear equations A*X=B or A'*X=B, using the LU factorization from sgstrf(). Error bounds on the solution and a condition estimate are also pro- vided.
#include <sunperf.h> void sgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r, int *etree, char *equed, float *R, float *C, Super- Matrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix *B, SuperMatrix *X, float *recip_pivot_growth, float *rcond, float *ferr, float *berr, mem_usage_t *mem_usage, SuperLUStat_t *stat, int *info) void dgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r, int *etree, char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMa- trix *B, SuperMatrix *X, double *recip_pivot_growth, double *rcond, double *ferr, double *berr, mem_usage_t *mem_usage, SuperLUStat_t *stat, int *info) void cgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r, int *etree, char *equed, float *R, float *C, Super- Matrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix *B, SuperMatrix *X, float *recip_pivot_growth, float *rcond, float *ferr, float *berr, mem_usage_t *mem_usage, SuperLUStat_t *stat, int *info) void zgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r, int *etree, char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMa- trix *B, SuperMatrix *X, double *recip_pivot_growth, double *rcond, double *ferr, double *berr, mem_usage_t *mem_usage, SuperLUStat_t *stat, int *info) void sgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long *perm_c, long *perm_r, long *etree, char *equed, float *R, float *C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, float *recip_pivot_growth, float *rcond, float *ferr, float *berr, mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info) void dgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long *perm_c, long *perm_r, long *etree, char *equed, double *R, dou- ble *C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, double *recip_pivot_growth, double *rcond, double *ferr, double *berr, mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info) void cgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long *perm_c, long *perm_r, long *etree, char *equed, float *R, float *C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, float *recip_pivot_growth, float *rcond, float *ferr, float *berr, mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info) void zgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long *perm_c, long *perm_r, long *etree, char *equed, double *R, dou- ble *C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, double *recip_pivot_growth, double *rcond, double *ferr, double *berr, mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info)
Oracle Solaris Studio Performance Library gssvx(3P)
NAME
gssvx: cgssvx, dgssvx, sgssvx, zgssvx - solves the system of linear
equations A*X=B or A'*X=B, using the LU factorization from sgstrf().
Error bounds on the solution and a condition estimate are also pro-
vided.
SYNOPSIS
#include <sunperf.h>
void sgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c,
int *perm_r, int *etree, char *equed, float *R, float *C, Super-
Matrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix
*B, SuperMatrix *X, float *recip_pivot_growth, float *rcond,
float *ferr, float *berr, mem_usage_t *mem_usage, SuperLUStat_t
*stat, int *info)
void dgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c,
int *perm_r, int *etree, char *equed, double *R, double *C,
SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMa-
trix *B, SuperMatrix *X, double *recip_pivot_growth, double
*rcond, double *ferr, double *berr, mem_usage_t *mem_usage,
SuperLUStat_t *stat, int *info)
void cgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c,
int *perm_r, int *etree, char *equed, float *R, float *C, Super-
Matrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix
*B, SuperMatrix *X, float *recip_pivot_growth, float *rcond,
float *ferr, float *berr, mem_usage_t *mem_usage, SuperLUStat_t
*stat, int *info)
void zgssvx(superlu_options_t *options, SuperMatrix *A, int *perm_c,
int *perm_r, int *etree, char *equed, double *R, double *C,
SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMa-
trix *B, SuperMatrix *X, double *recip_pivot_growth, double
*rcond, double *ferr, double *berr, mem_usage_t *mem_usage,
SuperLUStat_t *stat, int *info)
void sgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long
*perm_c, long *perm_r, long *etree, char *equed, float *R, float
*C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long
lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, float
*recip_pivot_growth, float *rcond, float *ferr, float *berr,
mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info)
void dgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long
*perm_c, long *perm_r, long *etree, char *equed, double *R, dou-
ble *C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long
lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, double
*recip_pivot_growth, double *rcond, double *ferr, double *berr,
mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info)
void cgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long
*perm_c, long *perm_r, long *etree, char *equed, float *R, float
*C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long
lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, float
*recip_pivot_growth, float *rcond, float *ferr, float *berr,
mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info)
void zgssvx_64(superlu_options_t_64 *options, SuperMatrix_64 *A, long
*perm_c, long *perm_r, long *etree, char *equed, double *R, dou-
ble *C, SuperMatrix_64 *L, SuperMatrix_64 *U, void *work, long
lwork, SuperMatrix_64 *B, SuperMatrix_64 *X, double
*recip_pivot_growth, double *rcond, double *ferr, double *berr,
mem_usage_t_64 *mem_usage, SuperLUStat_t_64 *stat, long *info)
PURPOSE
gssvx solves the system of linear equations A*X=B or A'*X=B, using the
LU factorization from sgstrf(). Error bounds on the solution and a con-
dition estimate are also provided.
If A is stored column-wise (A->Stype = SLU_NC):
o If options->Equil = YES, scaling factors are computed to equili-
brate the system:
options->Trans = NOTRANS:
diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
options->Trans = TRANS:
(diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
options->Trans = CONJ:
(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
options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS
or CONJ).
o Permute columns of A, forming A*Pc, where Pc is a permutation
matrix that usually preserves sparsity. For more details of this
step, see man page sp_preorder.
o If options->Fact != FACTORED, the LU decomposition is used to fac-
tor the matrix A (after equilibration if options->Equil = YES) as
Pr*A*Pc = L*U, with Pr determined by partial pivoting.
o Compute the reciprocal pivot growth factor.
o 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 =
A->ncol+1 is returned as a warning, but the routine still goes on
to solve for X and computes error bounds as described below.
o The system of equations is solved for X using the factored form of
A.
o If options->IterRefine != NOREFINE, iterative refinement is applied
to improve the computed solution matrix and calculate error bounds
and backward error estimates for it.
o If equilibration was used, the matrix X is premultiplied by diag(C)
(if options->Trans = NOTRANS) or diag(R) (if options->Trans = TRANS
or CONJ) so that it solves the original system before equilibra-
tion.
If A is stored row-wise (A->Stype = SLU_NR), apply the above algorithm
to the transpose of A:
o If options->Equil = YES, scaling factors are computed to equili-
brate the system:
options->Trans = NOTRANS:
diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
options->Trans = TRANS:
(diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
options->Trans = CONJ:
(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
options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS
or CONJ).
o Permute columns of transpose(A) (rows of A), forming trans-
pose(A)*Pc, where Pc is a permutation matrix that usually preserves
sparsity. For more details of this step, see man page sp_preorder.
o If options->Fact != FACTORED, the LU decomposition is used to fac-
tor the transpose(A) (after equilibration if options->Fact = YES)
as Pr*transpose(A)*Pc = L*U with the permutation Pr determined by
partial pivoting.
o Compute the reciprocal pivot growth factor.
o If some U(i,i) = 0, so that U is exactly singular, then the routine
returns with info = i. Otherwise, the factored form of transpose(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 = A->nrow+1 is returned as a warning, but the routine still
goes on to solve for X and computes error bounds as described
below.
o The system of equations is solved for X using the factored form of
transpose(A).
o If options->IterRefine != NOREFINE, iterative refinement is applied
to improve the computed solution matrix and calculate error bounds
and backward error estimates for it.
o If equilibration was used, the matrix X is premultiplied by diag(C)
(if options->Trans = NOTRANS) or diag(R) (if options->Trans = TRANS
or CONJ) so that it solves the original system before equilibra-
tion.
ARGUMENTS
superlu_options_t *options (input)
The structure defines the input parameters to control how the LU
decomposition will be performed and how the system will be
solved.
SuperMatrix *A (input/output)
Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
of linear equations is A->nrow. Currently, the type of A can be:
Stype = SLU_NC; Dtype = SLU_S; Mtype = SLU_GE.
In the future, more general A may be handled.
On entry, If options->Fact = FACTORED and equed is not 'N', then
A must have been equilibrated by the scaling factors in R and/or
C.
On exit, A is not modified if options->Equil = NO, or if
options->Equil = YES but equed = 'N' on exit.
Otherwise, if options->Equil = YES and equed is not 'N', A is
scaled as follows:
o If A->Stype = SLU_NC:
equed = 'R': A := diag(R) * A
equed = 'C': A := A * diag(C)
equed = 'B': A := diag(R) * A * diag(C).
o If A->Stype = SLU_NR:
equed = 'R': transpose(A) := diag(R) * transpose(A)
equed = 'C': transpose(A) := transpose(A) * diag(C)
equed = 'B': transpose(A) := diag(R) * transpose(A) *
diag(C).
int *perm_c (input/output)
If A->Stype = SLU_NC, column permutation vector of size
A->ncol which defines the permutation matrix Pc; perm_c[i]
= j means column i of A is in position j in A*Pc.
On exit, perm_c may be overwritten by the product of the
input perm_c and a permutation that postorders the elimina-
tion tree of Pc'*A'*A*Pc; perm_c is not changed if the
elimination tree is already in postorder.
If A->Stype = SLU_NR, column permutation vector of size
A->nrow which describes permutation of columns of trans-
pose(A) (rows of A) as described above.
int *perm_r (input/output)
If A->Stype = SLU_NC, row permutation vector of size
A->nrow, 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.
If A->Stype = SLU_NR, permutation vector of size A->ncol,
which determines permutation of rows of transpose(A) (col-
umns of A) as described above.
If options->Fact = SamePattern_SameRowPerm, the pivoting
routine will try to use the input perm_r, unless a certain
threshold criterion is violated. In that case, perm_r is
overwritten by a new permutation determined by partial piv-
oting or diagonal threshold pivoting.
Otherwise it is an output argument.
int *etree (input/output) An array of size (A->ncol) that contains
the elimination tree of Pc'*A'*A*Pc.
If options->Fact != FACTORED and options->Fact != DOFACT,
etree is an input argument; otherwise it is an output argu-
ment.
Note: etree is a vector of parent pointers for a forest
whose vertices are the integers 0 to A->ncol-1;
etree[root]==A->ncol.
char *equed (input/output)
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).
If options->Fact = FACTORED, equed is an input argu-
ment; otherwise it is an output argument.
float *R (input/output)
The row scale factors for A or transpose(A); dimension
(A->nrow).
If equed = 'R' or 'B', A (A->Stype = SLU_NC) or transpose(A)
(A->Stype = SLU_NR) is multiplied on the left by diag(R).
If equed = 'N' or 'C', R is not accessed.
If options->Fact = FACTORED, R is an input argument; otherwise,
R is output.
If options->zFact = FACTORED and equed = 'R' or 'B', each ele-
ment of R must be positive.
float *C (input/output)
The column scale factors for A or transpose(A); dimension
(A->ncol).
If equed = 'C' or 'B', A (A->Stype = SLU_NC) or transpose(A)
(A->Stype = SLU_NR) is multiplied on the right by diag(C).
If equed = 'N' or 'R', C is not accessed.
If options->Fact = FACTORED, C is an input argument; otherwise,
C is output.
If options->Fact = FACTORED and equed = 'C' or 'B', each element
of C must be positive.
SuperMatrix *L (output)
The factor L from the factorization
Pr*A*Pc=L*U (if A->Stype = SLU_NC) or
Pr*transpose(A)*Pc=L*U (if A->Stype = SLU_NR). Uses com-
pressed row subscripts storage for supernodes, i.e.,
L has types: Stype = SLU_SC, Dtype = SLU_S, Mtype = SLU_TRLU.
SuperMatrix *U (output)
The factor U from the factorization
Pr*A*Pc=L*U (if A->Stype = SLU_NC) or
Pr*transpose(A)*Pc=L*U (if A->Stype = SLU_NR).
Uses column-wise storage scheme, i.e., U has types: Stype =
SLU_NC, Dtype = SLU_S, Mtype = SLU_TRU.
void *work (workspace/output)
User-supplied workspace of length lwork, should be large enough
to hold data structures for factors L and U.
On exit, if fact is not 'F', L and U point to this array.
int lwork (input)
Specifies the size of work array in bytes.
= 0: allocate space internally by system malloc
> 0: use user-supplied work array of length lwork in bytes,
returns error if space runs out.
= -1: the routine guesses the amount of space needed without
performing the factorization, and returns it in info; no other
side effects.
SuperMatrix *B (input/output)
On entry, the right hand side matrix B.
On exit, the solution matrix if info = 0.
B has types: Stype = SLU_DN, Dtype = SLU_S, Mtype = SLU_GE.
If B->ncol = 0, only LU decomposition is performed, the triangu-
lar solve is skipped.
On exit,
if equed = 'N', B is not modified
if A->Stype = SLU_NC:
if options->Trans = NOTRANS and equed = 'R' or 'B', B is over-
written by diag(R)*B;
if options->Trans = TRANS or CONJ and equed = 'C' of 'B', B is
overwritten by diag(C)*B;
if A->Stype = SLU_NR:
if options->Trans = NOTRANS and equed = 'C' or 'B', B is over-
written by diag(C)*B;
if options->Trans = TRANS or CONJ and equed = 'R' of 'B', B is
overwritten by diag(R)*B.
SuperMatrix *X (output)
X has types: Stype = SLU_DN, Dtype = SLU_C, Mtype = SLU_GE.
If info = 0 or info = A->ncol+1, X contains the solution matrix
to the original system of equations. Note that A and B are modi-
fied on exit if equed is not 'N', and the solution to the equi-
librated system is inv(diag(C))*X if options->Trans = NOTRANS
and equed = 'C' or 'B', or inv(diag(R))*X if options->Trans =
'T' or 'C' and equed = 'R' or 'B'.
float *recip_pivot_growth (output)
The reciprocal pivot growth factor max_j( norm(A_j)/norm(U_j) ).
The infinity-norm is used. If recip_pivot_growth is much less
than 1, the stability of the LU factorization could be poor.
float *rcond (output) 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 indi-
cated by a return code of info > 0.
float *ferr (output)
dimension (B->ncol)
* 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 options->IterRefine = NOREFINE, ferr = 1.0.
float *berr (output)
dimension (B->ncol)
The componentwise relative backward error of each solution vec-
tor X(j) (i.e., the smallest relative change in any element of A
or B that makes X(j) an exact solution).
If options->IterRefine = NOREFINE, berr = 1.0.
mem_usage_t *mem_usage (output)
Record the memory usage statistics, consisting of following
fields:
o for_lu (float)
o The amount of space used in bytes for LU data structures.
o total_needed (float)
o The amount of space needed in bytes to perform factorization.
o expansions (int)
o The number of memory expansions during the LU factorization.
SuperLUStat_t *stat (output)
Records the statistics on runtime and floating-point operation
count.
int *info (output)
= 0:
successful exit
< 0: if info = -i, the i-th argument had an illegal value
> 0: if info = i, and i is
<= A->ncol: U(i,i) 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.
= A->ncol+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 the condition number wouldsuggest.
> A->ncol+1: number of bytes allocated when memory alloca-
tion failure occurred, plus A->ncol.
SEE ALSO
SuperMatrix
set_default_options
StatInit
StatFree
gstrf
sp_preorder
http://crd.lbl.gov/~xiaoye/SuperLU/
James W. Demmel, Stanley C. Eisenstat, John R. Gilbert, Xiaoye S. Li
and Joseph W. H. Liu, "A supernodal approach to sparse partial pivot-
ing", SIAM J. Matrix Analysis and Applications, Vol. 20, Num. 3, 1999,
pp. 720-755.
7 Nov 2015 gssvx(3P)