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)