diff options
Diffstat (limited to 'SRC/cgesvx.f')
| -rw-r--r-- | SRC/cgesvx.f | 571 |
1 files changed, 344 insertions, 227 deletions
diff --git a/SRC/cgesvx.f b/SRC/cgesvx.f index d675baab..49c317f3 100644 --- a/SRC/cgesvx.f +++ b/SRC/cgesvx.f @@ -1,11 +1,353 @@ +*> \brief <b> CGESVX computes the solution to system of linear equations A * X = B for GE matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition +* ========== +* +* SUBROUTINE CGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, +* EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, +* WORK, RWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER EQUED, FACT, TRANS +* INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS +* REAL RCOND +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* REAL BERR( * ), C( * ), FERR( * ), R( * ), +* $ RWORK( * ) +* COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), +* $ WORK( * ), X( LDX, * ) +* .. +* +* Purpose +* ======= +* +*>\details \b Purpose: +*>\verbatim +*> +*> CGESVX uses the LU factorization to compute the solution to a complex +*> system of linear equations +*> A * X = B, +*> where A is an N-by-N matrix and X and B are N-by-NRHS matrices. +*> +*> Error bounds on the solution and a condition estimate are also +*> provided. +*> +*> Description +*> =========== +*> +*> The following steps are performed: +*> +*> 1. If FACT = 'E', real scaling factors are 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))**T *inv(diag(R))*X = diag(C)*B +*> TRANS = 'C': (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 TRANS='N') +*> or diag(C)*B (if TRANS = 'T' or 'C'). +*> +*> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the +*> matrix A (after equilibration if FACT = 'E') as +*> A = P * L * U, +*> where P is a permutation matrix, L is a unit lower triangular +*> matrix, and U is upper triangular. +*> +*> 3. 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 = N+1 is returned as a warning, but the routine still goes on +*> to solve for X and compute error bounds as described below. +*> +*> 4. The system of equations is solved for X using the factored form +*> of A. +*> +*> 5. Iterative refinement is applied to improve the computed solution +*> matrix and calculate error bounds and backward error estimates +*> for it. +*> +*> 6. If equilibration was used, the matrix X is premultiplied by +*> diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so +*> that it solves the original system before equilibration. +*> +*>\endverbatim +* +* Arguments +* ========= +* +*> \param[in] FACT +*> \verbatim +*> FACT is CHARACTER*1 +*> Specifies whether or not the factored form of the matrix A is +*> supplied on entry, and if not, whether the matrix A should be +*> equilibrated before it is factored. +*> = 'F': On entry, AF and IPIV contain the factored form of A. +*> If EQUED is not 'N', the matrix A has been +*> equilibrated with scaling factors given by R and C. +*> A, AF, and IPIV are not modified. +*> = 'N': The matrix A will be copied to AF and factored. +*> = 'E': The matrix A will be equilibrated if necessary, then +*> copied to AF and factored. +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> Specifies the form of the system of equations: +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of linear equations, i.e., the order of the +*> matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrices B and X. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is +*> not 'N', then A must have been equilibrated by the scaling +*> factors in R and/or C. A is not modified if FACT = 'F' or +*> 'N', or if FACT = 'E' and EQUED = 'N' on exit. +*> \endverbatim +*> \verbatim +*> On exit, if EQUED .ne. 'N', A is scaled as follows: +*> EQUED = 'R': A := diag(R) * A +*> EQUED = 'C': A := A * diag(C) +*> EQUED = 'B': A := diag(R) * A * diag(C). +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in,out] AF +*> \verbatim +*> AF is or output) COMPLEX array, dimension (LDAF,N) +*> If FACT = 'F', then AF is an input argument and on entry +*> contains the factors L and U from the factorization +*> A = P*L*U as computed by CGETRF. If EQUED .ne. 'N', then +*> AF is the factored form of the equilibrated matrix A. +*> \endverbatim +*> \verbatim +*> If FACT = 'N', then AF is an output argument and on exit +*> returns the factors L and U from the factorization A = P*L*U +*> of the original matrix A. +*> \endverbatim +*> \verbatim +*> If FACT = 'E', then AF is an output argument and on exit +*> returns the factors L and U from the factorization A = P*L*U +*> of the equilibrated matrix A (see the description of A for +*> the form of the equilibrated matrix). +*> \endverbatim +*> +*> \param[in] LDAF +*> \verbatim +*> LDAF is INTEGER +*> The leading dimension of the array AF. LDAF >= max(1,N). +*> \endverbatim +*> +*> \param[in,out] IPIV +*> \verbatim +*> IPIV is or output) INTEGER array, dimension (N) +*> If FACT = 'F', then IPIV is an input argument and on entry +*> contains the pivot indices from the factorization A = P*L*U +*> as computed by CGETRF; row i of the matrix was interchanged +*> with row IPIV(i). +*> \endverbatim +*> \verbatim +*> If FACT = 'N', then IPIV is an output argument and on exit +*> contains the pivot indices from the factorization A = P*L*U +*> of the original matrix A. +*> \endverbatim +*> \verbatim +*> If FACT = 'E', then IPIV is an output argument and on exit +*> contains the pivot indices from the factorization A = P*L*U +*> of the equilibrated matrix A. +*> \endverbatim +*> +*> \param[in,out] EQUED +*> \verbatim +*> EQUED is or output) CHARACTER*1 +*> Specifies the form of equilibration that was done. +*> = 'N': No equilibration (always true if FACT = 'N'). +*> = '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). +*> EQUED is an input argument if FACT = 'F'; otherwise, it is an +*> output argument. +*> \endverbatim +*> +*> \param[in,out] R +*> \verbatim +*> R is or output) REAL array, dimension (N) +*> The row scale factors for A. If EQUED = 'R' or 'B', A is +*> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R +*> is not accessed. R is an input argument if FACT = 'F'; +*> otherwise, R is an output argument. If FACT = 'F' and +*> EQUED = 'R' or 'B', each element of R must be positive. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is or output) REAL array, dimension (N) +*> The column scale factors for A. If EQUED = 'C' or 'B', A is +*> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C +*> is not accessed. C is an input argument if FACT = 'F'; +*> otherwise, C is an output argument. If FACT = 'F' and +*> EQUED = 'C' or 'B', each element of C must be positive. +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX array, dimension (LDB,NRHS) +*> On entry, the N-by-NRHS right hand side matrix B. +*> On exit, +*> if EQUED = 'N', B is not modified; +*> if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by +*> diag(R)*B; +*> if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is +*> overwritten by diag(C)*B. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX array, dimension (LDX,NRHS) +*> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X +*> to the original system of equations. Note that A and B are +*> modified on exit if EQUED .ne. 'N', and the solution to the +*> equilibrated system is inv(diag(C))*X if TRANS = 'N' and +*> EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' +*> and EQUED = 'R' or 'B'. +*> \endverbatim +*> +*> \param[in] LDX +*> \verbatim +*> LDX is INTEGER +*> The leading dimension of the array X. LDX >= max(1,N). +*> \endverbatim +*> +*> \param[out] RCOND +*> \verbatim +*> RCOND is REAL +*> 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 +*> indicated by a return code of INFO > 0. +*> \endverbatim +*> +*> \param[out] FERR +*> \verbatim +*> FERR is REAL array, dimension (NRHS) +*> 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. +*> \endverbatim +*> +*> \param[out] BERR +*> \verbatim +*> BERR is REAL array, dimension (NRHS) +*> The componentwise relative backward error of each solution +*> vector X(j) (i.e., the smallest relative change in +*> any element of A or B that makes X(j) an exact solution). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (2*N) +*> On exit, RWORK(1) contains the reciprocal pivot growth +*> factor norm(A)/norm(U). The "max absolute element" norm is +*> used. If RWORK(1) is much less than 1, then the stability +*> of the LU factorization of the (equilibrated) matrix A +*> could be poor. This also means that the solution X, condition +*> estimator RCOND, and forward error bound FERR could be +*> unreliable. If factorization fails with 0<INFO<=N, then +*> RWORK(1) contains the reciprocal pivot growth factor for the +*> leading INFO columns of A. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, and i is +*> <= N: 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. RCOND = 0 is returned. +*> = N+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 RCOND would suggest. +*> \endverbatim +*> +* +* Authors +* ======= +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexGEsolve +* +* ===================================================================== SUBROUTINE CGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, $ EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, $ WORK, RWORK, INFO ) * -* -- LAPACK driver routine (version 3.2) -- +* -- LAPACK solve routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. CHARACTER EQUED, FACT, TRANS @@ -20,231 +362,6 @@ $ WORK( * ), X( LDX, * ) * .. * -* Purpose -* ======= -* -* CGESVX uses the LU factorization to compute the solution to a complex -* system of linear equations -* A * X = B, -* where A is an N-by-N matrix and X and B are N-by-NRHS matrices. -* -* Error bounds on the solution and a condition estimate are also -* provided. -* -* Description -* =========== -* -* The following steps are performed: -* -* 1. If FACT = 'E', real scaling factors are 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))**T *inv(diag(R))*X = diag(C)*B -* TRANS = 'C': (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 TRANS='N') -* or diag(C)*B (if TRANS = 'T' or 'C'). -* -* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the -* matrix A (after equilibration if FACT = 'E') as -* A = P * L * U, -* where P is a permutation matrix, L is a unit lower triangular -* matrix, and U is upper triangular. -* -* 3. 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 = N+1 is returned as a warning, but the routine still goes on -* to solve for X and compute error bounds as described below. -* -* 4. The system of equations is solved for X using the factored form -* of A. -* -* 5. Iterative refinement is applied to improve the computed solution -* matrix and calculate error bounds and backward error estimates -* for it. -* -* 6. If equilibration was used, the matrix X is premultiplied by -* diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so -* that it solves the original system before equilibration. -* -* Arguments -* ========= -* -* FACT (input) CHARACTER*1 -* Specifies whether or not the factored form of the matrix A is -* supplied on entry, and if not, whether the matrix A should be -* equilibrated before it is factored. -* = 'F': On entry, AF and IPIV contain the factored form of A. -* If EQUED is not 'N', the matrix A has been -* equilibrated with scaling factors given by R and C. -* A, AF, and IPIV are not modified. -* = 'N': The matrix A will be copied to AF and factored. -* = 'E': The matrix A will be equilibrated if necessary, then -* copied to AF and factored. -* -* TRANS (input) CHARACTER*1 -* Specifies the form of the system of equations: -* = 'N': A * X = B (No transpose) -* = 'T': A**T * X = B (Transpose) -* = 'C': A**H * X = B (Conjugate transpose) -* -* N (input) INTEGER -* The number of linear equations, i.e., the order of the -* matrix A. N >= 0. -* -* NRHS (input) INTEGER -* The number of right hand sides, i.e., the number of columns -* of the matrices B and X. NRHS >= 0. -* -* A (input/output) COMPLEX array, dimension (LDA,N) -* On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is -* not 'N', then A must have been equilibrated by the scaling -* factors in R and/or C. A is not modified if FACT = 'F' or -* 'N', or if FACT = 'E' and EQUED = 'N' on exit. -* -* On exit, if EQUED .ne. 'N', A is scaled as follows: -* EQUED = 'R': A := diag(R) * A -* EQUED = 'C': A := A * diag(C) -* EQUED = 'B': A := diag(R) * A * diag(C). -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* AF (input or output) COMPLEX array, dimension (LDAF,N) -* If FACT = 'F', then AF is an input argument and on entry -* contains the factors L and U from the factorization -* A = P*L*U as computed by CGETRF. If EQUED .ne. 'N', then -* AF is the factored form of the equilibrated matrix A. -* -* If FACT = 'N', then AF is an output argument and on exit -* returns the factors L and U from the factorization A = P*L*U -* of the original matrix A. -* -* If FACT = 'E', then AF is an output argument and on exit -* returns the factors L and U from the factorization A = P*L*U -* of the equilibrated matrix A (see the description of A for -* the form of the equilibrated matrix). -* -* LDAF (input) INTEGER -* The leading dimension of the array AF. LDAF >= max(1,N). -* -* IPIV (input or output) INTEGER array, dimension (N) -* If FACT = 'F', then IPIV is an input argument and on entry -* contains the pivot indices from the factorization A = P*L*U -* as computed by CGETRF; row i of the matrix was interchanged -* with row IPIV(i). -* -* If FACT = 'N', then IPIV is an output argument and on exit -* contains the pivot indices from the factorization A = P*L*U -* of the original matrix A. -* -* If FACT = 'E', then IPIV is an output argument and on exit -* contains the pivot indices from the factorization A = P*L*U -* of the equilibrated matrix A. -* -* EQUED (input or output) CHARACTER*1 -* Specifies the form of equilibration that was done. -* = 'N': No equilibration (always true if FACT = 'N'). -* = '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). -* EQUED is an input argument if FACT = 'F'; otherwise, it is an -* output argument. -* -* R (input or output) REAL array, dimension (N) -* The row scale factors for A. If EQUED = 'R' or 'B', A is -* multiplied on the left by diag(R); if EQUED = 'N' or 'C', R -* is not accessed. R is an input argument if FACT = 'F'; -* otherwise, R is an output argument. If FACT = 'F' and -* EQUED = 'R' or 'B', each element of R must be positive. -* -* C (input or output) REAL array, dimension (N) -* The column scale factors for A. If EQUED = 'C' or 'B', A is -* multiplied on the right by diag(C); if EQUED = 'N' or 'R', C -* is not accessed. C is an input argument if FACT = 'F'; -* otherwise, C is an output argument. If FACT = 'F' and -* EQUED = 'C' or 'B', each element of C must be positive. -* -* B (input/output) COMPLEX array, dimension (LDB,NRHS) -* On entry, the N-by-NRHS right hand side matrix B. -* On exit, -* if EQUED = 'N', B is not modified; -* if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by -* diag(R)*B; -* if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is -* overwritten by diag(C)*B. -* -* LDB (input) INTEGER -* The leading dimension of the array B. LDB >= max(1,N). -* -* X (output) COMPLEX array, dimension (LDX,NRHS) -* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X -* to the original system of equations. Note that A and B are -* modified on exit if EQUED .ne. 'N', and the solution to the -* equilibrated system is inv(diag(C))*X if TRANS = 'N' and -* EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' -* and EQUED = 'R' or 'B'. -* -* LDX (input) INTEGER -* The leading dimension of the array X. LDX >= max(1,N). -* -* RCOND (output) REAL -* 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 -* indicated by a return code of INFO > 0. -* -* FERR (output) REAL array, dimension (NRHS) -* 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. -* -* BERR (output) REAL array, dimension (NRHS) -* The componentwise relative backward error of each solution -* vector X(j) (i.e., the smallest relative change in -* any element of A or B that makes X(j) an exact solution). -* -* WORK (workspace) COMPLEX array, dimension (2*N) -* -* RWORK (workspace/output) REAL array, dimension (2*N) -* On exit, RWORK(1) contains the reciprocal pivot growth -* factor norm(A)/norm(U). The "max absolute element" norm is -* used. If RWORK(1) is much less than 1, then the stability -* of the LU factorization of the (equilibrated) matrix A -* could be poor. This also means that the solution X, condition -* estimator RCOND, and forward error bound FERR could be -* unreliable. If factorization fails with 0<INFO<=N, then -* RWORK(1) contains the reciprocal pivot growth factor for the -* leading INFO columns of A. -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: if INFO = i, and i is -* <= N: 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. RCOND = 0 is returned. -* = N+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 RCOND would suggest. -* * ===================================================================== * * .. Parameters .. |