Computes the solution to the system of linear equations with a square coefficient matrix A and multiple right-hand sides, and provides error bounds on the solution.
The routine uses the
LUfactorization to compute the solution to a real or complex system of linear equations
nmatrix, 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.
?gesvxperforms the following steps:
- If, real scaling factorsfact='E'randcare computed to equilibrate the system::trans='N'diag(r)*A*diag(c)*inv(diag(c))*X=diag(r)*B:trans='T'(diag(r)*A*diag(c))*inv(Tdiag(r))*X=diag(c)*B:trans='C'(diag(r)*A*diag(c))*inv(Hdiag(r))*X=diag(c)*BWhether or not the system will be equilibrated depends on the scaling of the matrixA, but if equilibration is used,Ais overwritten byanddiag(r)*A*diag(c)Bby(ifdiag(r)*Bortrans='N')(ifdiag(c)*Bortrans='T').'C'
- Iforfact='N', the'E'LUdecomposition is used to factor the matrixA(after equilibration if) asfact='E', whereA=P*L*UPis a permutation matrix,Lis a unit lower triangular matrix, andUis upper triangular.
- If some= 0, so thatUi,iUis exactly singular, 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 by(ifdiag(c)) ortrans='N'(ifdiag(r)ortrans='T''C') so that it solves the original system before equilibration.
- 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'afandipivcontain the factored form ofA. Ifequedis not'N', the matrixAhas been equilibrated with scaling factors given byrandc.a,af, andipivare not modified.If, the matrixfact='N'Awill be copied toafand factored.If, the matrixfact='E'Awill be equilibrated if necessary, then copied toafand factored.
- Must be'N','T', or'C'.Specifies the form of the system of equations:If, the system has the formtrans='N'(No transpose).A*X=BIf, the system has the formtrans='T'AT*X=B(Transpose).If, the system has the formtrans='C'AH*X=B(Transpose for real flavors, conjugate transpose for complex flavors).
- The number of linear equations; the order of the matrixA;n≥0.
- The number of right hand sides; the number of columns of the matricesBandX;nrhs≥0.
- The arraya(size max(1,contains the matrixlda*n))A. Ifandfact='F'equedis not'N', thenAmust have been equilibrated by the scaling factors inrand/orc.
- The arrayafaf(size max(1,is an input argument ifldaf*n)). It contains the factored form of the matrixfact='F'A, that is, the factorsLandUfrom the factorizationas computed byA=P*L*U?getrf. Ifequedis not'N', thenafis the factored form of the equilibrated matrixA.
- The arraybbof size max(1,contains the matrixldb*nrhs) for column major layout and max(1,ldb*n) for row major layoutBwhose columns are the right-hand sides for the systems of equations.
- The leading dimension ofa;.lda≥max(1,n)
- The leading dimension ofaf;.ldaf≥max(1,n)
- The leading dimension ofb;.ldb≥max(1,n) for column major layout andldb≥nrhsfor row major layout
- Array, size at leastmax(1,. The arrayn)ipivis an input argument if. It contains the pivot indices from the factorizationfact='F'as computed byA=P*L*U?getrf; rowiof the matrix was interchanged with row.ipiv[i-1]
- Must be'N','R','C', or'B'.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, row equilibration was done, that is,equed='R'Ahas been premultiplied bydiag(r).If, column equilibration was done, that is,equed='C'Ahas been postmultiplied bydiag(c).If, both row and column equilibration was done, that is,equed='B'Ahas been replaced by.diag(r)*A*diag(c)
- Arrays:,r(sizen). The arrayc(sizen)rcontains the row scale factors forA, and the arrayccontains the column scale factors forA. These arrays are input arguments ifonly; otherwise they are output arguments.fact='F'Iforequed='R''B',Ais multiplied on the left bydiag(r); iforequed='N''C',ris not accessed.Ifandfact='F'orequed='R''B', each element ofrmust be positive.Iforequed='C''B',Ais multiplied on the right bydiag(c); iforequed='N''R',cis not accessed.Ifandfact='F'orequed='C''B', each element ofcmust 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 thatAandBare modified on exit if, and the solution to the equilibrated system is:equed≠'N', ifdiag(C)-1*Xandtrans='N'orequed='C''B';, ifdiag(R)-1*Xortrans='T''C'andorequed='R''B'. The second dimension ofxmust be at leastmax(1,.nrhs)
- Arrayais not modified on exit iforfact='F''N', or ifandfact='E'equed='N'. If,equed≠'N'Ais scaled on exit as follows:equed='R':A=diag(R)*Aequed='C':A=A*diag(c)equed='B':.A=diag(R)*A*diag(c)
- Iforfact='N''E', thenafis an output argument and on exit returns the factorsLandUfrom the factorizationof the original matrixA=PLUA(if) or of the equilibrated matrixfact='N'A(if). See the description offact='E'afor the form of the equilibrated matrix.
- Overwritten byifdiag(r)*Bandtrans='N'orequed='R';'B'overwritten byifdiag(c)*Bortrans='T'and'C'orequed='C';'B'not changed ifequed='N'.
- These arrays are output arguments if. See the description offact≠'F'r,cinInput 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
- Iforfact='N''E', thenipivis an output argument and on exit contains the pivot indices from the factorizationof the original matrixA=P*L*UA(if) or of the equilibrated matrixfact='N'A(if).fact='E'
- If, thenfact≠'F'equedis an output argument. It specifies the form of equilibration that was done (see the description ofequedinInput Argumentssection).
- On exit,contains the reciprocal pivot growth factor:rpivotIfis much less than 1, then the stability of therpivotLUfactorization of the (equilibrated) matrixAcould be poor. This also means that the solutionx, condition estimatorrcond, and forward error boundferrcould be unreliable. If factorization fails with0 <, theninfo≤ncontains the reciprocal pivot growth factor for the leadingrpivotinfocolumns ofA.
This function returns a value
, the execution is successful.
ihad an illegal value.
i) is exactly zero. The factorization has been completed, but the factor
Uis exactly singular, so the solution and error bounds could not be computed;
rcond= 0 is returned.
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 computed solution can be more accurate than the value of