Uses the Cholesky factorization to compute the solution to the system of linear equations with a symmetric (Hermitian) positive definite packed coefficient matrix A, and provides error bounds on the solution.
The routine uses the Cholesky factorization
(real flavors) /
(complex flavors) or
(real flavors) /
(complex flavors) to compute the solution to a real or complex system of linear equations
nsymmetric or Hermitian positive-definite matrix stored in packed format, the columns of matrix
Bare individual right-hand sides, and the columns of
Xare the corresponding solutions.
Error bounds on the solution and a condition estimate are also provided.
?ppsvxperforms 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).Bbydiag(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)
- Specifies whether matrix storage layout is row major (LAPACK_ROW_MAJOR) or column major (LAPACK_COL_MAJOR).
- 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'afpcontains the factored form ofA. If, the matrixequed='Y'Ahas been equilibrated with scaling factors given bys.apandafpwill not be modified.If, the matrixfact='N'Awill be copied toafpand factored.If, the matrixfact='E'Awill be equilibrated if necessary, then copied toafpand factored.
- 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.
- The order of matrixA;n≥0.
- The number of right-hand sides; the number of columns in.B;nrhs≥0
- Arrays:(size max(1,,n*(n+1)/2)afp(size max(1,,n*(n+1)/2)bof size max(1,.ldb*nrhs) for column major layout and max(1,ldb*n) for row major layoutThe arrayafpis an input argument ifand contains the triangular factorfact='F'UorLfrom the Cholesky factorization ofAin the same storage format asA. Ifequedis not'N', thenafpis the factored form of the equilibrated matrixA.The arraybcontains the matrixBwhose columns are the right-hand sides for the systems of equations.
- The leading dimension ofb;.ldb≥max(1,n) for column major layout andldb≥nrhsfor row major layout
- 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)
- 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.
- The leading dimension of the output arrayx;.ldx≥max(1,n) for column major layout andldx≥nrhsfor row major layout
- 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
- Arrayapis not modified on exit iforfact='F''N', or ifandfact='E'equed='N'.Ifandfact='E'equed='Y',apis overwritten bydiag(s)*A*diag(s).
- Iforfact='N''E', thenafpis 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'apfor the form of the equilibrated matrix.
- Overwritten by, ifdiag(s)*B; not changed ifequed='Y'equed='N'.
- This array is an output argument if. See the description offact≠'F'sinInput Argumentssection.
- 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.
- 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 inferr[j-1](divided by the magnitude of the largest element inx-jxtrue). The estimate is as reliable as the estimate forxjrcond, and is almost always a slight overestimate of the true error.
- 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
- If, thenfact≠'F'equedis an output argument. It specifies the form of equilibration that was done (see the description ofequedinInput Argumentssection).
This function returns a value
, the execution is successful.
ihad an illegal value.
, the leading minor of order
i(and therefore the matrix
Aitself) 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.
n+ 1, then
Uis nonsingular, but
rcondis 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