?posvx
?posvx
Uses the Cholesky factorization to compute the solution to the system of linear equations with a symmetric or Hermitian positive-definite coefficient matrix A, and provides error bounds on the solution.
Syntax
lapack_int LAPACKE_sposvx
(
int
matrix_layout
,
char
fact
,
char
uplo
,
lapack_int
n
,
lapack_int
nrhs
,
float*
a
,
lapack_int
lda
,
float*
af
,
lapack_int
ldaf
,
char*
equed
,
float*
s
,
float*
b
,
lapack_int
ldb
,
float*
x
,
lapack_int
ldx
,
float*
rcond
,
float*
ferr
,
float*
berr
);
lapack_int LAPACKE_dposvx
(
int
matrix_layout
,
char
fact
,
char
uplo
,
lapack_int
n
,
lapack_int
nrhs
,
double*
a
,
lapack_int
lda
,
double*
af
,
lapack_int
ldaf
,
char*
equed
,
double*
s
,
double*
b
,
lapack_int
ldb
,
double*
x
,
lapack_int
ldx
,
double*
rcond
,
double*
ferr
,
double*
berr
);
lapack_int LAPACKE_cposvx
(
int
matrix_layout
,
char
fact
,
char
uplo
,
lapack_int
n
,
lapack_int
nrhs
,
lapack_complex_float*
a
,
lapack_int
lda
,
lapack_complex_float*
af
,
lapack_int
ldaf
,
char*
equed
,
float*
s
,
lapack_complex_float*
b
,
lapack_int
ldb
,
lapack_complex_float*
x
,
lapack_int
ldx
,
float*
rcond
,
float*
ferr
,
float*
berr
);
lapack_int LAPACKE_zposvx
(
int
matrix_layout
,
char
fact
,
char
uplo
,
lapack_int
n
,
lapack_int
nrhs
,
lapack_complex_double*
a
,
lapack_int
lda
,
lapack_complex_double*
af
,
lapack_int
ldaf
,
char*
equed
,
double*
s
,
lapack_complex_double*
b
,
lapack_int
ldb
,
lapack_complex_double*
x
,
lapack_int
ldx
,
double*
rcond
,
double*
ferr
,
double*
berr
);
Include Files
- mkl.h
Description
The routine uses the (real flavors) / (complex flavors) or (real flavors) / (complex flavors) to compute the solution to a real or complex system of linear equations , where
Cholesky
factorization A
=U
T
*U
A
=U
H
*U
A
=L*L
T
A
=L*L
H
A*X
= B
A
is a n
-by-n
real symmetric/Hermitian positive definite matrix, the columns of matrix B
are individual right-hand sides, and the columns of X
are the corresponding solutions.Error bounds on the solution and a condition estimate are also provided.
The routine
?posvx
performs the following steps:- If, real scaling factorsfact='E'sare computed to equilibrate the system:diag(s)*A*diag(s)*inv(diag(s))*X = diag(s)*B.Whether or not the system will be equilibrated depends on the scaling of the matrixA, but if equilibration is used,Ais overwritten byanddiag(s)*A*diag(s)Bby.diag(s)*B
- Iforfact='N''E', the Cholesky decomposition is used to factor the matrixA(after equilibration if) asfact='E'(real),A=UT*U(complex), ifA=UH*U,uplo='U'or(real),A=L*LT(complex), ifA=L*LH,uplo='L'whereUis an upper triangular matrix andLis a lower triangular matrix.
- If the leadingi-by-iprincipal minor is not positive-definite, then the routine returns with. Otherwise, the factored form ofinfo=iAis used to estimate the condition number of the matrixA. If the reciprocal of the condition number is less than machine precision,is returned as a warning, but the routine still goes on to solve forinfo=n+ 1Xand compute error bounds as described below.
- The system of equations is solved forXusing the factored form ofA.
- Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
- If equilibration was used, the matrixXis premultiplied byso that it solves the original system before equilibration.diag(s)
Input Parameters
- matrix_layout
- Specifies whether matrix storage layout is row major (LAPACK_ROW_MAJOR) or column major (LAPACK_COL_MAJOR).
- fact
- Must be'F','N', or'E'.Specifies whether or not the factored form of the matrixAis supplied on entry, and if not, whether the matrixAshould be equilibrated before it is factored.If: on entry,fact='F'afcontains the factored form ofA. If, the matrixequed='Y'Ahas been equilibrated with scaling factors given bys.aandafwill not be modified.If, the matrixfact='N'Awill be copied toafand factored.If, the matrixfact='E'Awill be equilibrated if necessary, then copied toafand factored.
- uplo
- Must be'U'or'L'.Indicates whether the upper or lower triangular part ofAis stored:If, the upper triangle ofuplo='U'Ais stored.If, the lower triangle ofuplo='L'Ais stored.
- n
- The order of matrixA;n≥0.
- nrhs
- The number of right-hand sides, the number of columns in.B;nrhs≥0
- a,af,b
- Arrays:a(size max(1,,lda*n))af(size max(1,,ldaf*n))b, size max(.ldb*nrhs) for column major layout and max(ldb*n) for row major layout,The arrayacontains the matrixAas specified byuplo. Ifandfact='F'equed='Y', thenAmust have been equilibrated by the scaling factors ins, andamust contain the equilibrated matrix.diag(s)*A*diag(s)The arrayafis an input argument if. It contains the triangular factorfact='F'UorLfrom the Cholesky factorization ofAin the same storage format asA. Ifequedis not'N', thenafis the factored form of the equilibrated matrix.diag(s)*A*diag(s)The arraybcontains the matrixBwhose columns are the right-hand sides for the systems of equations.
- lda
- The leading dimension ofa;.lda≥max(1,n)
- The leading dimension ofaf;.ldaf≥max(1,n)
- ldb
- The leading dimension ofb;.ldb≥max(1,n) for column major layout andldb≥nrhsfor row major layout
- equed
- Must be'N'or'Y'.equedis an input argument if. It specifies the form of equilibration that was done:fact='F'if, no equilibration was done (always true ifequed='N');fact='N'if, equilibration was done, that is,equed='Y'Ahas been replaced by.diag(s)*A*diag(s)
- s
- Array, size (n). The arrayscontains the scale factors forA. This array is an input argument ifonly; otherwise it is an output argument.fact='F'If,equed='N'sis not accessed.Ifandfact='F'equed='Y', each element ofsmust be positive.
- ldx
- The leading dimension of the output arrayx;.ldx≥max(1,n) for column major layout andldx≥nrhsfor row major layout
Output Parameters
- x
- Array, sizemax(1,.ldx*nrhs) for column major layout and max(1,ldx*n) for row major layoutIforinfo= 0, the arrayinfo=n+1xcontains the solution matrixXto theoriginalsystem of equations. Note that if,equed='Y'AandBare modified on exit, and the solution to the equilibrated system isinv(.diag(s))*X
- a
- Arrayais not modified on exit iforfact='F''N', or ifandfact='E'equed='N'.Ifandfact='E'equed='Y',Ais overwritten by.diag(s)*A*diag(s)
- af
- Iforfact='N''E', thenafis an output argument and on exit returns the triangular factorUorLfrom the Cholesky factorizationorA=UT*U(real routines),A=L*LTorA=UH*U(complex routines) of the original matrixA=L*LHA(if), or of the equilibrated matrixfact='N'A(if). See the description offact='E'afor the form of the equilibrated matrix.
- b
- Overwritten by, ifdiag(s)*B; not changed ifequed='Y'equed='N'.
- s
- This array is an output argument if. See the description offact≠'F'sinInput Argumentssection.
- rcond
- An estimate of the reciprocal condition number of the matrixAafter equilibration (if done). Ifrcondis less than the machine precision (in particular, ifrcond=0), the matrix is singular to working precision. This condition is indicated by a return code ofinfo>0.
- ferr
- Array, size at leastmax(1,. Contains the estimated forward error bound for each solution vectornrhs)(thexjj-th column of the solution matrixX). Ifxtrueis the true solution corresponding to,xjis an estimated upper bound for the magnitude of the largest element in (ferr[j-1]) -xjxtrue) divided by the magnitude of the largest element in. The estimate is as reliable as the estimate forxjrcond, and is almost always a slight overestimate of the true error.
- berr
- Array, size at leastmax(1,. Contains the component-wise relative backward error for each solution vectornrhs), that is, the smallest relative change in any element ofxjAorBthat makesan exact solution.xj
- equed
- If, thenfact≠'F'equedis an output argument. It specifies the form of equilibration that was done (see the description ofequedinInput Argumentssection).
Return Values
This function returns a value
info
.If , the execution is successful.
info
= 0If , parameter
info
= -i
i
had an illegal value. If , and , the leading minor of order
info
= i
i
≤
n
i
(and therefore the matrix A
itself) is not positive-definite, so the factorization could not be completed, and the solution and error bounds could not be computed; rcond
=0 is returned. If , and
info
= i
i
= n
+ 1, then 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.