gssvx(3P) Sun 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)
PURPOSEgssvx 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 elimination 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 transpose(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) (columns 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 pivoting 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 argument.
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 argument;
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.
COPYRIGHT
Copyright (c) 2003, The Regents of the University of California,
through Lawrence Berkeley National Laboratory (subject to receipt of
any required approvals from U.S. Dept. of Energy)
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.
6 Mar 2009 gssvx(3P)
