![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZPOSVX computes the solution to system of linear equations A * X = B for PO 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 ZPOSVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zposvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zposvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zposvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 00022 * S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 00023 * RWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER EQUED, FACT, UPLO 00027 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00028 * DOUBLE PRECISION RCOND 00029 * .. 00030 * .. Array Arguments .. 00031 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * ) 00032 * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00033 * $ WORK( * ), X( LDX, * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to 00043 *> compute the solution to a complex system of linear equations 00044 *> A * X = B, 00045 *> where A is an N-by-N Hermitian positive definite matrix and X and B 00046 *> are N-by-NRHS matrices. 00047 *> 00048 *> Error bounds on the solution and a condition estimate are also 00049 *> provided. 00050 *> \endverbatim 00051 * 00052 *> \par Description: 00053 * ================= 00054 *> 00055 *> \verbatim 00056 *> 00057 *> The following steps are performed: 00058 *> 00059 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate 00060 *> the system: 00061 *> diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B 00062 *> Whether or not the system will be equilibrated depends on the 00063 *> scaling of the matrix A, but if equilibration is used, A is 00064 *> overwritten by diag(S)*A*diag(S) and B by diag(S)*B. 00065 *> 00066 *> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to 00067 *> factor the matrix A (after equilibration if FACT = 'E') as 00068 *> A = U**H* U, if UPLO = 'U', or 00069 *> A = L * L**H, if UPLO = 'L', 00070 *> where U is an upper triangular matrix and L is a lower triangular 00071 *> matrix. 00072 *> 00073 *> 3. If the leading i-by-i principal minor is not positive definite, 00074 *> then the routine returns with INFO = i. Otherwise, the factored 00075 *> form of A is used to estimate the condition number of the matrix 00076 *> A. If the reciprocal of the condition number is less than machine 00077 *> precision, INFO = N+1 is returned as a warning, but the routine 00078 *> still goes on to solve for X and compute error bounds as 00079 *> described below. 00080 *> 00081 *> 4. The system of equations is solved for X using the factored form 00082 *> of A. 00083 *> 00084 *> 5. Iterative refinement is applied to improve the computed solution 00085 *> matrix and calculate error bounds and backward error estimates 00086 *> for it. 00087 *> 00088 *> 6. If equilibration was used, the matrix X is premultiplied by 00089 *> diag(S) so that it solves the original system before 00090 *> equilibration. 00091 *> \endverbatim 00092 * 00093 * Arguments: 00094 * ========== 00095 * 00096 *> \param[in] FACT 00097 *> \verbatim 00098 *> FACT is CHARACTER*1 00099 *> Specifies whether or not the factored form of the matrix A is 00100 *> supplied on entry, and if not, whether the matrix A should be 00101 *> equilibrated before it is factored. 00102 *> = 'F': On entry, AF contains the factored form of A. 00103 *> If EQUED = 'Y', the matrix A has been equilibrated 00104 *> with scaling factors given by S. A and AF will not 00105 *> be modified. 00106 *> = 'N': The matrix A will be copied to AF and factored. 00107 *> = 'E': The matrix A will be equilibrated if necessary, then 00108 *> copied to AF and factored. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] UPLO 00112 *> \verbatim 00113 *> UPLO is CHARACTER*1 00114 *> = 'U': Upper triangle of A is stored; 00115 *> = 'L': Lower triangle of A is stored. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] N 00119 *> \verbatim 00120 *> N is INTEGER 00121 *> The number of linear equations, i.e., the order of the 00122 *> matrix A. N >= 0. 00123 *> \endverbatim 00124 *> 00125 *> \param[in] NRHS 00126 *> \verbatim 00127 *> NRHS is INTEGER 00128 *> The number of right hand sides, i.e., the number of columns 00129 *> of the matrices B and X. NRHS >= 0. 00130 *> \endverbatim 00131 *> 00132 *> \param[in,out] A 00133 *> \verbatim 00134 *> A is COMPLEX*16 array, dimension (LDA,N) 00135 *> On entry, the Hermitian matrix A, except if FACT = 'F' and 00136 *> EQUED = 'Y', then A must contain the equilibrated matrix 00137 *> diag(S)*A*diag(S). If UPLO = 'U', the leading 00138 *> N-by-N upper triangular part of A contains the upper 00139 *> triangular part of the matrix A, and the strictly lower 00140 *> triangular part of A is not referenced. If UPLO = 'L', the 00141 *> leading N-by-N lower triangular part of A contains the lower 00142 *> triangular part of the matrix A, and the strictly upper 00143 *> triangular part of A is not referenced. A is not modified if 00144 *> FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. 00145 *> 00146 *> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by 00147 *> diag(S)*A*diag(S). 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 triangular factor U or L from the Cholesky 00161 *> factorization A = U**H *U or A = L*L**H, in the same storage 00162 *> format as A. If EQUED .ne. 'N', then AF is the factored form 00163 *> of the equilibrated matrix diag(S)*A*diag(S). 00164 *> 00165 *> If FACT = 'N', then AF is an output argument and on exit 00166 *> returns the triangular factor U or L from the Cholesky 00167 *> factorization A = U**H *U or A = L*L**H of the original 00168 *> matrix A. 00169 *> 00170 *> If FACT = 'E', then AF is an output argument and on exit 00171 *> returns the triangular factor U or L from the Cholesky 00172 *> factorization A = U**H *U or A = L*L**H of the equilibrated 00173 *> matrix A (see the description of A for the form of the 00174 *> equilibrated matrix). 00175 *> \endverbatim 00176 *> 00177 *> \param[in] LDAF 00178 *> \verbatim 00179 *> LDAF is INTEGER 00180 *> The leading dimension of the array AF. LDAF >= max(1,N). 00181 *> \endverbatim 00182 *> 00183 *> \param[in,out] EQUED 00184 *> \verbatim 00185 *> EQUED is CHARACTER*1 00186 *> Specifies the form of equilibration that was done. 00187 *> = 'N': No equilibration (always true if FACT = 'N'). 00188 *> = 'Y': Equilibration was done, i.e., A has been replaced by 00189 *> diag(S) * A * diag(S). 00190 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an 00191 *> output argument. 00192 *> \endverbatim 00193 *> 00194 *> \param[in,out] S 00195 *> \verbatim 00196 *> S is DOUBLE PRECISION array, dimension (N) 00197 *> The scale factors for A; not accessed if EQUED = 'N'. S is 00198 *> an input argument if FACT = 'F'; otherwise, S is an output 00199 *> argument. If FACT = 'F' and EQUED = 'Y', each element of S 00200 *> must be positive. 00201 *> \endverbatim 00202 *> 00203 *> \param[in,out] B 00204 *> \verbatim 00205 *> B is COMPLEX*16 array, dimension (LDB,NRHS) 00206 *> On entry, the N-by-NRHS righthand side matrix B. 00207 *> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', 00208 *> B is overwritten by diag(S) * B. 00209 *> \endverbatim 00210 *> 00211 *> \param[in] LDB 00212 *> \verbatim 00213 *> LDB is INTEGER 00214 *> The leading dimension of the array B. LDB >= max(1,N). 00215 *> \endverbatim 00216 *> 00217 *> \param[out] X 00218 *> \verbatim 00219 *> X is COMPLEX*16 array, dimension (LDX,NRHS) 00220 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to 00221 *> the original system of equations. Note that if EQUED = 'Y', 00222 *> A and B are modified on exit, and the solution to the 00223 *> equilibrated system is inv(diag(S))*X. 00224 *> \endverbatim 00225 *> 00226 *> \param[in] LDX 00227 *> \verbatim 00228 *> LDX is INTEGER 00229 *> The leading dimension of the array X. LDX >= max(1,N). 00230 *> \endverbatim 00231 *> 00232 *> \param[out] RCOND 00233 *> \verbatim 00234 *> RCOND is DOUBLE PRECISION 00235 *> The estimate of the reciprocal condition number of the matrix 00236 *> A after equilibration (if done). If RCOND is less than the 00237 *> machine precision (in particular, if RCOND = 0), the matrix 00238 *> is singular to working precision. This condition is 00239 *> indicated by a return code of INFO > 0. 00240 *> \endverbatim 00241 *> 00242 *> \param[out] FERR 00243 *> \verbatim 00244 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00245 *> The estimated forward error bound for each solution vector 00246 *> X(j) (the j-th column of the solution matrix X). 00247 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00248 *> is an estimated upper bound for the magnitude of the largest 00249 *> element in (X(j) - XTRUE) divided by the magnitude of the 00250 *> largest element in X(j). The estimate is as reliable as 00251 *> the estimate for RCOND, and is almost always a slight 00252 *> overestimate of the true error. 00253 *> \endverbatim 00254 *> 00255 *> \param[out] BERR 00256 *> \verbatim 00257 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00258 *> The componentwise relative backward error of each solution 00259 *> vector X(j) (i.e., the smallest relative change in 00260 *> any element of A or B that makes X(j) an exact solution). 00261 *> \endverbatim 00262 *> 00263 *> \param[out] WORK 00264 *> \verbatim 00265 *> WORK is COMPLEX*16 array, dimension (2*N) 00266 *> \endverbatim 00267 *> 00268 *> \param[out] RWORK 00269 *> \verbatim 00270 *> RWORK is DOUBLE PRECISION array, dimension (N) 00271 *> \endverbatim 00272 *> 00273 *> \param[out] INFO 00274 *> \verbatim 00275 *> INFO is INTEGER 00276 *> = 0: successful exit 00277 *> < 0: if INFO = -i, the i-th argument had an illegal value 00278 *> > 0: if INFO = i, and i is 00279 *> <= N: the leading minor of order i of A is 00280 *> not positive definite, so the factorization 00281 *> could not be completed, and the solution has not 00282 *> been computed. RCOND = 0 is returned. 00283 *> = N+1: U is nonsingular, but RCOND is less than machine 00284 *> precision, meaning that the matrix is singular 00285 *> to working precision. Nevertheless, the 00286 *> solution and error bounds are computed because 00287 *> there are a number of situations where the 00288 *> computed solution can be more accurate than the 00289 *> value of RCOND would suggest. 00290 *> \endverbatim 00291 * 00292 * Authors: 00293 * ======== 00294 * 00295 *> \author Univ. of Tennessee 00296 *> \author Univ. of California Berkeley 00297 *> \author Univ. of Colorado Denver 00298 *> \author NAG Ltd. 00299 * 00300 *> \date April 2012 00301 * 00302 *> \ingroup complex16POsolve 00303 * 00304 * ===================================================================== 00305 SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 00306 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 00307 $ RWORK, INFO ) 00308 * 00309 * -- LAPACK driver routine (version 3.4.1) -- 00310 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00311 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00312 * April 2012 00313 * 00314 * .. Scalar Arguments .. 00315 CHARACTER EQUED, FACT, UPLO 00316 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 00317 DOUBLE PRECISION RCOND 00318 * .. 00319 * .. Array Arguments .. 00320 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * ) 00321 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00322 $ WORK( * ), X( LDX, * ) 00323 * .. 00324 * 00325 * ===================================================================== 00326 * 00327 * .. Parameters .. 00328 DOUBLE PRECISION ZERO, ONE 00329 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00330 * .. 00331 * .. Local Scalars .. 00332 LOGICAL EQUIL, NOFACT, RCEQU 00333 INTEGER I, INFEQU, J 00334 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 00335 * .. 00336 * .. External Functions .. 00337 LOGICAL LSAME 00338 DOUBLE PRECISION DLAMCH, ZLANHE 00339 EXTERNAL LSAME, DLAMCH, ZLANHE 00340 * .. 00341 * .. External Subroutines .. 00342 EXTERNAL XERBLA, ZLACPY, ZLAQHE, ZPOCON, ZPOEQU, ZPORFS, 00343 $ ZPOTRF, ZPOTRS 00344 * .. 00345 * .. Intrinsic Functions .. 00346 INTRINSIC MAX, MIN 00347 * .. 00348 * .. Executable Statements .. 00349 * 00350 INFO = 0 00351 NOFACT = LSAME( FACT, 'N' ) 00352 EQUIL = LSAME( FACT, 'E' ) 00353 IF( NOFACT .OR. EQUIL ) THEN 00354 EQUED = 'N' 00355 RCEQU = .FALSE. 00356 ELSE 00357 RCEQU = LSAME( EQUED, 'Y' ) 00358 SMLNUM = DLAMCH( 'Safe minimum' ) 00359 BIGNUM = ONE / SMLNUM 00360 END IF 00361 * 00362 * Test the input parameters. 00363 * 00364 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 00365 $ THEN 00366 INFO = -1 00367 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 00368 $ THEN 00369 INFO = -2 00370 ELSE IF( N.LT.0 ) THEN 00371 INFO = -3 00372 ELSE IF( NRHS.LT.0 ) THEN 00373 INFO = -4 00374 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00375 INFO = -6 00376 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00377 INFO = -8 00378 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 00379 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 00380 INFO = -9 00381 ELSE 00382 IF( RCEQU ) THEN 00383 SMIN = BIGNUM 00384 SMAX = ZERO 00385 DO 10 J = 1, N 00386 SMIN = MIN( SMIN, S( J ) ) 00387 SMAX = MAX( SMAX, S( J ) ) 00388 10 CONTINUE 00389 IF( SMIN.LE.ZERO ) THEN 00390 INFO = -10 00391 ELSE IF( N.GT.0 ) THEN 00392 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 00393 ELSE 00394 SCOND = ONE 00395 END IF 00396 END IF 00397 IF( INFO.EQ.0 ) THEN 00398 IF( LDB.LT.MAX( 1, N ) ) THEN 00399 INFO = -12 00400 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00401 INFO = -14 00402 END IF 00403 END IF 00404 END IF 00405 * 00406 IF( INFO.NE.0 ) THEN 00407 CALL XERBLA( 'ZPOSVX', -INFO ) 00408 RETURN 00409 END IF 00410 * 00411 IF( EQUIL ) THEN 00412 * 00413 * Compute row and column scalings to equilibrate the matrix A. 00414 * 00415 CALL ZPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU ) 00416 IF( INFEQU.EQ.0 ) THEN 00417 * 00418 * Equilibrate the matrix. 00419 * 00420 CALL ZLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 00421 RCEQU = LSAME( EQUED, 'Y' ) 00422 END IF 00423 END IF 00424 * 00425 * Scale the right hand side. 00426 * 00427 IF( RCEQU ) THEN 00428 DO 30 J = 1, NRHS 00429 DO 20 I = 1, N 00430 B( I, J ) = S( I )*B( I, J ) 00431 20 CONTINUE 00432 30 CONTINUE 00433 END IF 00434 * 00435 IF( NOFACT .OR. EQUIL ) THEN 00436 * 00437 * Compute the Cholesky factorization A = U**H *U or A = L*L**H. 00438 * 00439 CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 00440 CALL ZPOTRF( UPLO, N, AF, LDAF, INFO ) 00441 * 00442 * Return if INFO is non-zero. 00443 * 00444 IF( INFO.GT.0 )THEN 00445 RCOND = ZERO 00446 RETURN 00447 END IF 00448 END IF 00449 * 00450 * Compute the norm of the matrix A. 00451 * 00452 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 00453 * 00454 * Compute the reciprocal of the condition number of A. 00455 * 00456 CALL ZPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO ) 00457 * 00458 * Compute the solution matrix X. 00459 * 00460 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00461 CALL ZPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 00462 * 00463 * Use iterative refinement to improve the computed solution and 00464 * compute error bounds and backward error estimates for it. 00465 * 00466 CALL ZPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, 00467 $ FERR, BERR, WORK, RWORK, INFO ) 00468 * 00469 * Transform the solution matrix X to a solution of the original 00470 * system. 00471 * 00472 IF( RCEQU ) THEN 00473 DO 50 J = 1, NRHS 00474 DO 40 I = 1, N 00475 X( I, J ) = S( I )*X( I, J ) 00476 40 CONTINUE 00477 50 CONTINUE 00478 DO 60 J = 1, NRHS 00479 FERR( J ) = FERR( J ) / SCOND 00480 60 CONTINUE 00481 END IF 00482 * 00483 * Set INFO = N+1 if the matrix is singular to working precision. 00484 * 00485 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00486 $ INFO = N + 1 00487 * 00488 RETURN 00489 * 00490 * End of ZPOSVX 00491 * 00492 END