![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZGESVX computes the solution to system of linear equations A * X = B for GE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZGESVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, 00022 * EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, 00023 * WORK, RWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER EQUED, FACT, TRANS 00027 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00028 * DOUBLE PRECISION RCOND 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IPIV( * ) 00032 * DOUBLE PRECISION BERR( * ), C( * ), FERR( * ), R( * ), 00033 * $ RWORK( * ) 00034 * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00035 * $ WORK( * ), X( LDX, * ) 00036 * .. 00037 * 00038 * 00039 *> \par Purpose: 00040 * ============= 00041 *> 00042 *> \verbatim 00043 *> 00044 *> ZGESVX uses the LU factorization to compute the solution to a complex 00045 *> system of linear equations 00046 *> A * X = B, 00047 *> where A is an N-by-N matrix and X and B are N-by-NRHS matrices. 00048 *> 00049 *> Error bounds on the solution and a condition estimate are also 00050 *> provided. 00051 *> \endverbatim 00052 * 00053 *> \par Description: 00054 * ================= 00055 *> 00056 *> \verbatim 00057 *> 00058 *> The following steps are performed: 00059 *> 00060 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate 00061 *> the system: 00062 *> TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B 00063 *> TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B 00064 *> TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B 00065 *> Whether or not the system will be equilibrated depends on the 00066 *> scaling of the matrix A, but if equilibration is used, A is 00067 *> overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') 00068 *> or diag(C)*B (if TRANS = 'T' or 'C'). 00069 *> 00070 *> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the 00071 *> matrix A (after equilibration if FACT = 'E') as 00072 *> A = P * L * U, 00073 *> where P is a permutation matrix, L is a unit lower triangular 00074 *> matrix, and U is upper triangular. 00075 *> 00076 *> 3. If some U(i,i)=0, so that U is exactly singular, then the routine 00077 *> returns with INFO = i. Otherwise, the factored form of A is used 00078 *> to estimate the condition number of the matrix A. If the 00079 *> reciprocal of the condition number is less than machine precision, 00080 *> INFO = N+1 is returned as a warning, but the routine still goes on 00081 *> to solve for X and compute error bounds as described below. 00082 *> 00083 *> 4. The system of equations is solved for X using the factored form 00084 *> of A. 00085 *> 00086 *> 5. Iterative refinement is applied to improve the computed solution 00087 *> matrix and calculate error bounds and backward error estimates 00088 *> for it. 00089 *> 00090 *> 6. If equilibration was used, the matrix X is premultiplied by 00091 *> diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so 00092 *> that it solves the original system before equilibration. 00093 *> \endverbatim 00094 * 00095 * Arguments: 00096 * ========== 00097 * 00098 *> \param[in] FACT 00099 *> \verbatim 00100 *> FACT is CHARACTER*1 00101 *> Specifies whether or not the factored form of the matrix A is 00102 *> supplied on entry, and if not, whether the matrix A should be 00103 *> equilibrated before it is factored. 00104 *> = 'F': On entry, AF and IPIV contain the factored form of A. 00105 *> If EQUED is not 'N', the matrix A has been 00106 *> equilibrated with scaling factors given by R and C. 00107 *> A, AF, and IPIV are not modified. 00108 *> = 'N': The matrix A will be copied to AF and factored. 00109 *> = 'E': The matrix A will be equilibrated if necessary, then 00110 *> copied to AF and factored. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] TRANS 00114 *> \verbatim 00115 *> TRANS is CHARACTER*1 00116 *> Specifies the form of the system of equations: 00117 *> = 'N': A * X = B (No transpose) 00118 *> = 'T': A**T * X = B (Transpose) 00119 *> = 'C': A**H * X = B (Conjugate transpose) 00120 *> \endverbatim 00121 *> 00122 *> \param[in] N 00123 *> \verbatim 00124 *> N is INTEGER 00125 *> The number of linear equations, i.e., the order of the 00126 *> matrix A. N >= 0. 00127 *> \endverbatim 00128 *> 00129 *> \param[in] NRHS 00130 *> \verbatim 00131 *> NRHS is INTEGER 00132 *> The number of right hand sides, i.e., the number of columns 00133 *> of the matrices B and X. NRHS >= 0. 00134 *> \endverbatim 00135 *> 00136 *> \param[in,out] A 00137 *> \verbatim 00138 *> A is COMPLEX*16 array, dimension (LDA,N) 00139 *> On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is 00140 *> not 'N', then A must have been equilibrated by the scaling 00141 *> factors in R and/or C. A is not modified if FACT = 'F' or 00142 *> 'N', or if FACT = 'E' and EQUED = 'N' on exit. 00143 *> 00144 *> On exit, if EQUED .ne. 'N', A is scaled as follows: 00145 *> EQUED = 'R': A := diag(R) * A 00146 *> EQUED = 'C': A := A * diag(C) 00147 *> EQUED = 'B': A := diag(R) * A * diag(C). 00148 *> \endverbatim 00149 *> 00150 *> \param[in] LDA 00151 *> \verbatim 00152 *> LDA is INTEGER 00153 *> The leading dimension of the array A. LDA >= max(1,N). 00154 *> \endverbatim 00155 *> 00156 *> \param[in,out] AF 00157 *> \verbatim 00158 *> AF is COMPLEX*16 array, dimension (LDAF,N) 00159 *> If FACT = 'F', then AF is an input argument and on entry 00160 *> contains the factors L and U from the factorization 00161 *> A = P*L*U as computed by ZGETRF. If EQUED .ne. 'N', then 00162 *> AF is the factored form of the equilibrated matrix A. 00163 *> 00164 *> If FACT = 'N', then AF is an output argument and on exit 00165 *> returns the factors L and U from the factorization A = P*L*U 00166 *> of the original matrix A. 00167 *> 00168 *> If FACT = 'E', then AF is an output argument and on exit 00169 *> returns the factors L and U from the factorization A = P*L*U 00170 *> of the equilibrated matrix A (see the description of A for 00171 *> the form of the equilibrated matrix). 00172 *> \endverbatim 00173 *> 00174 *> \param[in] LDAF 00175 *> \verbatim 00176 *> LDAF is INTEGER 00177 *> The leading dimension of the array AF. LDAF >= max(1,N). 00178 *> \endverbatim 00179 *> 00180 *> \param[in,out] IPIV 00181 *> \verbatim 00182 *> IPIV is INTEGER array, dimension (N) 00183 *> If FACT = 'F', then IPIV is an input argument and on entry 00184 *> contains the pivot indices from the factorization A = P*L*U 00185 *> as computed by ZGETRF; row i of the matrix was interchanged 00186 *> with row IPIV(i). 00187 *> 00188 *> If FACT = 'N', then IPIV is an output argument and on exit 00189 *> contains the pivot indices from the factorization A = P*L*U 00190 *> of the original matrix A. 00191 *> 00192 *> If FACT = 'E', then IPIV is an output argument and on exit 00193 *> contains the pivot indices from the factorization A = P*L*U 00194 *> of the equilibrated matrix A. 00195 *> \endverbatim 00196 *> 00197 *> \param[in,out] EQUED 00198 *> \verbatim 00199 *> EQUED is CHARACTER*1 00200 *> Specifies the form of equilibration that was done. 00201 *> = 'N': No equilibration (always true if FACT = 'N'). 00202 *> = 'R': Row equilibration, i.e., A has been premultiplied by 00203 *> diag(R). 00204 *> = 'C': Column equilibration, i.e., A has been postmultiplied 00205 *> by diag(C). 00206 *> = 'B': Both row and column equilibration, i.e., A has been 00207 *> replaced by diag(R) * A * diag(C). 00208 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an 00209 *> output argument. 00210 *> \endverbatim 00211 *> 00212 *> \param[in,out] R 00213 *> \verbatim 00214 *> R is DOUBLE PRECISION array, dimension (N) 00215 *> The row scale factors for A. If EQUED = 'R' or 'B', A is 00216 *> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R 00217 *> is not accessed. R is an input argument if FACT = 'F'; 00218 *> otherwise, R is an output argument. If FACT = 'F' and 00219 *> EQUED = 'R' or 'B', each element of R must be positive. 00220 *> \endverbatim 00221 *> 00222 *> \param[in,out] C 00223 *> \verbatim 00224 *> C is DOUBLE PRECISION array, dimension (N) 00225 *> The column scale factors for A. If EQUED = 'C' or 'B', A is 00226 *> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C 00227 *> is not accessed. C is an input argument if FACT = 'F'; 00228 *> otherwise, C is an output argument. If FACT = 'F' and 00229 *> EQUED = 'C' or 'B', each element of C must be positive. 00230 *> \endverbatim 00231 *> 00232 *> \param[in,out] B 00233 *> \verbatim 00234 *> B is COMPLEX*16 array, dimension (LDB,NRHS) 00235 *> On entry, the N-by-NRHS right hand side matrix B. 00236 *> On exit, 00237 *> if EQUED = 'N', B is not modified; 00238 *> if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by 00239 *> diag(R)*B; 00240 *> if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is 00241 *> overwritten by diag(C)*B. 00242 *> \endverbatim 00243 *> 00244 *> \param[in] LDB 00245 *> \verbatim 00246 *> LDB is INTEGER 00247 *> The leading dimension of the array B. LDB >= max(1,N). 00248 *> \endverbatim 00249 *> 00250 *> \param[out] X 00251 *> \verbatim 00252 *> X is COMPLEX*16 array, dimension (LDX,NRHS) 00253 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X 00254 *> to the original system of equations. Note that A and B are 00255 *> modified on exit if EQUED .ne. 'N', and the solution to the 00256 *> equilibrated system is inv(diag(C))*X if TRANS = 'N' and 00257 *> EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' 00258 *> and EQUED = 'R' or 'B'. 00259 *> \endverbatim 00260 *> 00261 *> \param[in] LDX 00262 *> \verbatim 00263 *> LDX is INTEGER 00264 *> The leading dimension of the array X. LDX >= max(1,N). 00265 *> \endverbatim 00266 *> 00267 *> \param[out] RCOND 00268 *> \verbatim 00269 *> RCOND is DOUBLE PRECISION 00270 *> The estimate of the reciprocal condition number of the matrix 00271 *> A after equilibration (if done). If RCOND is less than the 00272 *> machine precision (in particular, if RCOND = 0), the matrix 00273 *> is singular to working precision. This condition is 00274 *> indicated by a return code of INFO > 0. 00275 *> \endverbatim 00276 *> 00277 *> \param[out] FERR 00278 *> \verbatim 00279 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00280 *> The estimated forward error bound for each solution vector 00281 *> X(j) (the j-th column of the solution matrix X). 00282 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00283 *> is an estimated upper bound for the magnitude of the largest 00284 *> element in (X(j) - XTRUE) divided by the magnitude of the 00285 *> largest element in X(j). The estimate is as reliable as 00286 *> the estimate for RCOND, and is almost always a slight 00287 *> overestimate of the true error. 00288 *> \endverbatim 00289 *> 00290 *> \param[out] BERR 00291 *> \verbatim 00292 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00293 *> The componentwise relative backward error of each solution 00294 *> vector X(j) (i.e., the smallest relative change in 00295 *> any element of A or B that makes X(j) an exact solution). 00296 *> \endverbatim 00297 *> 00298 *> \param[out] WORK 00299 *> \verbatim 00300 *> WORK is COMPLEX*16 array, dimension (2*N) 00301 *> \endverbatim 00302 *> 00303 *> \param[out] RWORK 00304 *> \verbatim 00305 *> RWORK is DOUBLE PRECISION array, dimension (2*N) 00306 *> On exit, RWORK(1) contains the reciprocal pivot growth 00307 *> factor norm(A)/norm(U). The "max absolute element" norm is 00308 *> used. If RWORK(1) is much less than 1, then the stability 00309 *> of the LU factorization of the (equilibrated) matrix A 00310 *> could be poor. This also means that the solution X, condition 00311 *> estimator RCOND, and forward error bound FERR could be 00312 *> unreliable. If factorization fails with 0<INFO<=N, then 00313 *> RWORK(1) contains the reciprocal pivot growth factor for the 00314 *> leading INFO columns of A. 00315 *> \endverbatim 00316 *> 00317 *> \param[out] INFO 00318 *> \verbatim 00319 *> INFO is INTEGER 00320 *> = 0: successful exit 00321 *> < 0: if INFO = -i, the i-th argument had an illegal value 00322 *> > 0: if INFO = i, and i is 00323 *> <= N: U(i,i) is exactly zero. The factorization has 00324 *> been completed, but the factor U is exactly 00325 *> singular, so the solution and error bounds 00326 *> could not be computed. RCOND = 0 is returned. 00327 *> = N+1: U is nonsingular, but RCOND is less than machine 00328 *> precision, meaning that the matrix is singular 00329 *> to working precision. Nevertheless, the 00330 *> solution and error bounds are computed because 00331 *> there are a number of situations where the 00332 *> computed solution can be more accurate than the 00333 *> value of RCOND would suggest. 00334 *> \endverbatim 00335 * 00336 * Authors: 00337 * ======== 00338 * 00339 *> \author Univ. of Tennessee 00340 *> \author Univ. of California Berkeley 00341 *> \author Univ. of Colorado Denver 00342 *> \author NAG Ltd. 00343 * 00344 *> \date April 2012 00345 * 00346 *> \ingroup complex16GEsolve 00347 * 00348 * ===================================================================== 00349 SUBROUTINE ZGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, 00350 $ EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, 00351 $ WORK, RWORK, INFO ) 00352 * 00353 * -- LAPACK driver routine (version 3.4.1) -- 00354 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00355 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00356 * April 2012 00357 * 00358 * .. Scalar Arguments .. 00359 CHARACTER EQUED, FACT, TRANS 00360 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00361 DOUBLE PRECISION RCOND 00362 * .. 00363 * .. Array Arguments .. 00364 INTEGER IPIV( * ) 00365 DOUBLE PRECISION BERR( * ), C( * ), FERR( * ), R( * ), 00366 $ RWORK( * ) 00367 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00368 $ WORK( * ), X( LDX, * ) 00369 * .. 00370 * 00371 * ===================================================================== 00372 * 00373 * .. Parameters .. 00374 DOUBLE PRECISION ZERO, ONE 00375 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00376 * .. 00377 * .. Local Scalars .. 00378 LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU 00379 CHARACTER NORM 00380 INTEGER I, INFEQU, J 00381 DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN, 00382 $ ROWCND, RPVGRW, SMLNUM 00383 * .. 00384 * .. External Functions .. 00385 LOGICAL LSAME 00386 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANTR 00387 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANTR 00388 * .. 00389 * .. External Subroutines .. 00390 EXTERNAL XERBLA, ZGECON, ZGEEQU, ZGERFS, ZGETRF, ZGETRS, 00391 $ ZLACPY, ZLAQGE 00392 * .. 00393 * .. Intrinsic Functions .. 00394 INTRINSIC MAX, MIN 00395 * .. 00396 * .. Executable Statements .. 00397 * 00398 INFO = 0 00399 NOFACT = LSAME( FACT, 'N' ) 00400 EQUIL = LSAME( FACT, 'E' ) 00401 NOTRAN = LSAME( TRANS, 'N' ) 00402 IF( NOFACT .OR. EQUIL ) THEN 00403 EQUED = 'N' 00404 ROWEQU = .FALSE. 00405 COLEQU = .FALSE. 00406 ELSE 00407 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' ) 00408 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' ) 00409 SMLNUM = DLAMCH( 'Safe minimum' ) 00410 BIGNUM = ONE / SMLNUM 00411 END IF 00412 * 00413 * Test the input parameters. 00414 * 00415 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 00416 $ THEN 00417 INFO = -1 00418 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 00419 $ LSAME( TRANS, 'C' ) ) THEN 00420 INFO = -2 00421 ELSE IF( N.LT.0 ) THEN 00422 INFO = -3 00423 ELSE IF( NRHS.LT.0 ) THEN 00424 INFO = -4 00425 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00426 INFO = -6 00427 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00428 INFO = -8 00429 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 00430 $ ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 00431 INFO = -10 00432 ELSE 00433 IF( ROWEQU ) THEN 00434 RCMIN = BIGNUM 00435 RCMAX = ZERO 00436 DO 10 J = 1, N 00437 RCMIN = MIN( RCMIN, R( J ) ) 00438 RCMAX = MAX( RCMAX, R( J ) ) 00439 10 CONTINUE 00440 IF( RCMIN.LE.ZERO ) THEN 00441 INFO = -11 00442 ELSE IF( N.GT.0 ) THEN 00443 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00444 ELSE 00445 ROWCND = ONE 00446 END IF 00447 END IF 00448 IF( COLEQU .AND. INFO.EQ.0 ) THEN 00449 RCMIN = BIGNUM 00450 RCMAX = ZERO 00451 DO 20 J = 1, N 00452 RCMIN = MIN( RCMIN, C( J ) ) 00453 RCMAX = MAX( RCMAX, C( J ) ) 00454 20 CONTINUE 00455 IF( RCMIN.LE.ZERO ) THEN 00456 INFO = -12 00457 ELSE IF( N.GT.0 ) THEN 00458 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00459 ELSE 00460 COLCND = ONE 00461 END IF 00462 END IF 00463 IF( INFO.EQ.0 ) THEN 00464 IF( LDB.LT.MAX( 1, N ) ) THEN 00465 INFO = -14 00466 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00467 INFO = -16 00468 END IF 00469 END IF 00470 END IF 00471 * 00472 IF( INFO.NE.0 ) THEN 00473 CALL XERBLA( 'ZGESVX', -INFO ) 00474 RETURN 00475 END IF 00476 * 00477 IF( EQUIL ) THEN 00478 * 00479 * Compute row and column scalings to equilibrate the matrix A. 00480 * 00481 CALL ZGEEQU( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFEQU ) 00482 IF( INFEQU.EQ.0 ) THEN 00483 * 00484 * Equilibrate the matrix. 00485 * 00486 CALL ZLAQGE( N, N, A, LDA, R, C, ROWCND, COLCND, AMAX, 00487 $ EQUED ) 00488 ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' ) 00489 COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' ) 00490 END IF 00491 END IF 00492 * 00493 * Scale the right hand side. 00494 * 00495 IF( NOTRAN ) THEN 00496 IF( ROWEQU ) THEN 00497 DO 40 J = 1, NRHS 00498 DO 30 I = 1, N 00499 B( I, J ) = R( I )*B( I, J ) 00500 30 CONTINUE 00501 40 CONTINUE 00502 END IF 00503 ELSE IF( COLEQU ) THEN 00504 DO 60 J = 1, NRHS 00505 DO 50 I = 1, N 00506 B( I, J ) = C( I )*B( I, J ) 00507 50 CONTINUE 00508 60 CONTINUE 00509 END IF 00510 * 00511 IF( NOFACT .OR. EQUIL ) THEN 00512 * 00513 * Compute the LU factorization of A. 00514 * 00515 CALL ZLACPY( 'Full', N, N, A, LDA, AF, LDAF ) 00516 CALL ZGETRF( N, N, AF, LDAF, IPIV, INFO ) 00517 * 00518 * Return if INFO is non-zero. 00519 * 00520 IF( INFO.GT.0 ) THEN 00521 * 00522 * Compute the reciprocal pivot growth factor of the 00523 * leading rank-deficient INFO columns of A. 00524 * 00525 RPVGRW = ZLANTR( 'M', 'U', 'N', INFO, INFO, AF, LDAF, 00526 $ RWORK ) 00527 IF( RPVGRW.EQ.ZERO ) THEN 00528 RPVGRW = ONE 00529 ELSE 00530 RPVGRW = ZLANGE( 'M', N, INFO, A, LDA, RWORK ) / 00531 $ RPVGRW 00532 END IF 00533 RWORK( 1 ) = RPVGRW 00534 RCOND = ZERO 00535 RETURN 00536 END IF 00537 END IF 00538 * 00539 * Compute the norm of the matrix A and the 00540 * reciprocal pivot growth factor RPVGRW. 00541 * 00542 IF( NOTRAN ) THEN 00543 NORM = '1' 00544 ELSE 00545 NORM = 'I' 00546 END IF 00547 ANORM = ZLANGE( NORM, N, N, A, LDA, RWORK ) 00548 RPVGRW = ZLANTR( 'M', 'U', 'N', N, N, AF, LDAF, RWORK ) 00549 IF( RPVGRW.EQ.ZERO ) THEN 00550 RPVGRW = ONE 00551 ELSE 00552 RPVGRW = ZLANGE( 'M', N, N, A, LDA, RWORK ) / RPVGRW 00553 END IF 00554 * 00555 * Compute the reciprocal of the condition number of A. 00556 * 00557 CALL ZGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO ) 00558 * 00559 * Compute the solution matrix X. 00560 * 00561 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00562 CALL ZGETRS( TRANS, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 00563 * 00564 * Use iterative refinement to improve the computed solution and 00565 * compute error bounds and backward error estimates for it. 00566 * 00567 CALL ZGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 00568 $ LDX, FERR, BERR, WORK, RWORK, INFO ) 00569 * 00570 * Transform the solution matrix X to a solution of the original 00571 * system. 00572 * 00573 IF( NOTRAN ) THEN 00574 IF( COLEQU ) THEN 00575 DO 80 J = 1, NRHS 00576 DO 70 I = 1, N 00577 X( I, J ) = C( I )*X( I, J ) 00578 70 CONTINUE 00579 80 CONTINUE 00580 DO 90 J = 1, NRHS 00581 FERR( J ) = FERR( J ) / COLCND 00582 90 CONTINUE 00583 END IF 00584 ELSE IF( ROWEQU ) THEN 00585 DO 110 J = 1, NRHS 00586 DO 100 I = 1, N 00587 X( I, J ) = R( I )*X( I, J ) 00588 100 CONTINUE 00589 110 CONTINUE 00590 DO 120 J = 1, NRHS 00591 FERR( J ) = FERR( J ) / ROWCND 00592 120 CONTINUE 00593 END IF 00594 * 00595 * Set INFO = N+1 if the matrix is singular to working precision. 00596 * 00597 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00598 $ INFO = N + 1 00599 * 00600 RWORK( 1 ) = RPVGRW 00601 RETURN 00602 * 00603 * End of ZGESVX 00604 * 00605 END