Function

NumCosmoMathlapack_dsysvxx

Declaration [src] 

gint
ncm_lapack_dsysvxx (
  gchar fact,
  gchar uplo,
  gint n,
  gint nrhs,
  gdouble* a,
  gint lda,
  gdouble* af,
  gint ldaf,
  gint* ipiv,
  gchar* equed,
  gdouble* s,
  gdouble* b,
  gint ldb,
  gdouble* x,
  gint ldx,
  gdouble* rcond,
  gdouble* rpvgrw,
  gdouble* berr,
  const gint n_err_bnds,
  gdouble* err_bnds_norm,
  gdouble* err_bnds_comp,
  const gint nparams,
  gdouble* params,
  gdouble* work,
  gint* iwork
)

Description [src]

Purpose

DSYSVXX uses the diagonal pivoting factorization to compute the solution to a double precision system of linear equations A * X = B, where A is an N-by-N symmetric matrix and X and B are N-by-NRHS matrices.

If requested, both norm-wise and maximum component-wise error bounds are returned. DSYSVXX will return a solution with a tiny guaranteed error (O(eps) where eps is the working machine precision) unless the matrix is very ill-conditioned, in which case a warning is returned. Relevant condition numbers also are calculated and returned.

DSYSVXX accepts user-provided factorizations and equilibration factors; see the definitions of the FACT and EQUED options. Solving with refinement and using a factorization from a previous DSYSVXX call will also produce a solution with either O(eps) errors or warnings, but we cannot make that claim for general user-provided factorizations and equilibration factors if they differ from what DSYSVXX would itself produce.

Description

The following steps are performed:

  1. If FACT = ‘E’, double precision scaling factors are computed to equilibrate the system:
  2. diag(S)Adiag(S) inv(diag(S))X = diag(S)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(S)Adiag(S) and B by diag(S)B.
  3. If FACT = ‘N’ or ‘E’, the LU decomposition is used to factor the matrix A (after equilibration if FACT = ‘E’) as
  4. A = U * D * U**T, if UPLO = ‘U’, or
  5. A = L * D * L**T, if UPLO = ‘L’, where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
  6. If some D(i,i)=0, so that D 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 (see argument RCOND). If the reciprocal of the condition number is less than machine precision, the routine still goes on to solve for X and compute error bounds as described below.
  7. The system of equations is solved for X using the factored form of A.
  8. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero), the routine will use iterative refinement to try to get a small error and error bounds. Refinement calculates the residual to at least twice the working precision.
  9. If equilibration was used, the matrix X is premultiplied by diag(R) so that it solves the original system before equilibration.

Some optional parameters are bundled in the PARAMS array. These settings determine how refinement is performed, but often the defaults are acceptable. If the defaults are acceptable, users can pass NPARAMS = 0 which prevents the source code from accessing the PARAMS argument.

Parameters

fact

Type: gchar

FACT is CHARACTER*1.

uplo

Type: gchar

UPLO is CHARACTER*1.

n

Type: gint

N is INTEGER.

nrhs

Type: gint

NRHS is INTEGER.

a

Type: gdouble*

A is DOUBLE PRECISION array, dimension (LDA,N).

The data is owned by the caller of the function.
lda

Type: gint

LDA is INTEGER.

af

Type: gdouble*

AF is DOUBLE PRECISION array, dimension (LDAF,N).

The data is owned by the caller of the function.
ldaf

Type: gint

LDA is INTEGER.

ipiv

Type: gint*

IPIV is INTEGER array, dimension (N).

The data is owned by the caller of the function.
equed

Type: gchar*

EQUED is CHARACTER*1.

The data is owned by the caller of the function.
The value is a NUL terminated UTF-8 string.
s

Type: gdouble*

S is DOUBLE PRECISION array, dimension (N).

The data is owned by the caller of the function.
b

Type: gdouble*

B is DOUBLE PRECISION array, dimension (LDB,NRHS).

The data is owned by the caller of the function.
ldb

Type: gint

LDB is INTEGER.

x

Type: gdouble*

X is DOUBLE PRECISION array, dimension (LDX,NRHS).

The data is owned by the caller of the function.
ldx

Type: gint

LDX is INTEGER.

rcond

Type: gdouble*

RCOND is DOUBLE PRECISION.

The data is owned by the caller of the function.
rpvgrw

Type: gdouble*

RPVGRW is DOUBLE PRECISION.

The data is owned by the caller of the function.
berr

Type: gdouble*

BERR is DOUBLE PRECISION array, dimension (NRHS).

The data is owned by the caller of the function.
n_err_bnds

Type: const gint

N_ERR_BNDS is INTEGER.

err_bnds_norm

Type: gdouble*

ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS).

The data is owned by the caller of the function.
err_bnds_comp

Type: gdouble*

ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS).

The data is owned by the caller of the function.
nparams

Type: const gint

NPARAMS is INTEGER.

params

Type: gdouble*

PARAMS is DOUBLE PRECISION array, dimension (NPARAMS).

The data is owned by the caller of the function.
work

Type: gdouble*

WORK is DOUBLE PRECISION array, dimension (4*N).

The data is owned by the caller of the function.
iwork

Type: gint*

IWORK is INTEGER array, dimension (N).

The data is owned by the caller of the function.

Return value

Type: gint

INFO is INTEGER - = 0: Successful exit. The solution to every right-hand side is guaranteed. - < 0: If INFO = -i, the i-th argument had an illegal value - > 0 and <= N: U(INFO,INFO) 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+J: The solution corresponding to the Jth right-hand side is not guaranteed. The solutions corresponding to other right- hand sides K with K > J may not be guaranteed as well, but only the first such right-hand side is reported. If a small componentwise error is not requested (PARAMS(3) = 0.0) then the Jth right-hand side is the first with a normwise error bound that is not guaranteed (the smallest J such that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) the Jth right-hand side is the first with either a normwise or componentwise error bound that is not guaranteed (the smallest J such that either ERR_BNDS_NORM(J,1) = 0.0 or ERR_BNDS_COMP(J,1) = 0.0). See the definition of ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information about all of the right-hand sides check ERR_BNDS_NORM or ERR_BNDS_COMP.