![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CPBSVX computes the solution to system of linear equations A * X = B for OTHER 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 CPBSVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpbsvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpbsvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpbsvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, 00022 * EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, 00023 * WORK, RWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER EQUED, FACT, UPLO 00027 * INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS 00028 * REAL RCOND 00029 * .. 00030 * .. Array Arguments .. 00031 * REAL BERR( * ), FERR( * ), RWORK( * ), S( * ) 00032 * COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 00033 * $ WORK( * ), X( LDX, * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> CPBSVX 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 band matrix and X 00046 *> and B 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 band matrix, and L is a lower 00071 *> triangular band 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, AFB contains the factored form of A. 00103 *> If EQUED = 'Y', the matrix A has been equilibrated 00104 *> with scaling factors given by S. AB and AFB will not 00105 *> be modified. 00106 *> = 'N': The matrix A will be copied to AFB and factored. 00107 *> = 'E': The matrix A will be equilibrated if necessary, then 00108 *> copied to AFB 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] KD 00126 *> \verbatim 00127 *> KD is INTEGER 00128 *> The number of superdiagonals of the matrix A if UPLO = 'U', 00129 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] NRHS 00133 *> \verbatim 00134 *> NRHS is INTEGER 00135 *> The number of right-hand sides, i.e., the number of columns 00136 *> of the matrices B and X. NRHS >= 0. 00137 *> \endverbatim 00138 *> 00139 *> \param[in,out] AB 00140 *> \verbatim 00141 *> AB is COMPLEX array, dimension (LDAB,N) 00142 *> On entry, the upper or lower triangle of the Hermitian band 00143 *> matrix A, stored in the first KD+1 rows of the array, except 00144 *> if FACT = 'F' and EQUED = 'Y', then A must contain the 00145 *> equilibrated matrix diag(S)*A*diag(S). The j-th column of A 00146 *> is stored in the j-th column of the array AB as follows: 00147 *> if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j; 00148 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD). 00149 *> See below for further details. 00150 *> 00151 *> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by 00152 *> diag(S)*A*diag(S). 00153 *> \endverbatim 00154 *> 00155 *> \param[in] LDAB 00156 *> \verbatim 00157 *> LDAB is INTEGER 00158 *> The leading dimension of the array A. LDAB >= KD+1. 00159 *> \endverbatim 00160 *> 00161 *> \param[in,out] AFB 00162 *> \verbatim 00163 *> AFB is COMPLEX array, dimension (LDAFB,N) 00164 *> If FACT = 'F', then AFB is an input argument and on entry 00165 *> contains the triangular factor U or L from the Cholesky 00166 *> factorization A = U**H*U or A = L*L**H of the band matrix 00167 *> A, in the same storage format as A (see AB). If EQUED = 'Y', 00168 *> then AFB is the factored form of the equilibrated matrix A. 00169 *> 00170 *> If FACT = 'N', then AFB 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. 00173 *> 00174 *> If FACT = 'E', then AFB is an output argument and on exit 00175 *> returns the triangular factor U or L from the Cholesky 00176 *> factorization A = U**H*U or A = L*L**H of the equilibrated 00177 *> matrix A (see the description of A for the form of the 00178 *> equilibrated matrix). 00179 *> \endverbatim 00180 *> 00181 *> \param[in] LDAFB 00182 *> \verbatim 00183 *> LDAFB is INTEGER 00184 *> The leading dimension of the array AFB. LDAFB >= KD+1. 00185 *> \endverbatim 00186 *> 00187 *> \param[in,out] EQUED 00188 *> \verbatim 00189 *> EQUED is CHARACTER*1 00190 *> Specifies the form of equilibration that was done. 00191 *> = 'N': No equilibration (always true if FACT = 'N'). 00192 *> = 'Y': Equilibration was done, i.e., A has been replaced by 00193 *> diag(S) * A * diag(S). 00194 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an 00195 *> output argument. 00196 *> \endverbatim 00197 *> 00198 *> \param[in,out] S 00199 *> \verbatim 00200 *> S is REAL array, dimension (N) 00201 *> The scale factors for A; not accessed if EQUED = 'N'. S is 00202 *> an input argument if FACT = 'F'; otherwise, S is an output 00203 *> argument. If FACT = 'F' and EQUED = 'Y', each element of S 00204 *> must be positive. 00205 *> \endverbatim 00206 *> 00207 *> \param[in,out] B 00208 *> \verbatim 00209 *> B is COMPLEX array, dimension (LDB,NRHS) 00210 *> On entry, the N-by-NRHS right hand side matrix B. 00211 *> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', 00212 *> B is overwritten by diag(S) * B. 00213 *> \endverbatim 00214 *> 00215 *> \param[in] LDB 00216 *> \verbatim 00217 *> LDB is INTEGER 00218 *> The leading dimension of the array B. LDB >= max(1,N). 00219 *> \endverbatim 00220 *> 00221 *> \param[out] X 00222 *> \verbatim 00223 *> X is COMPLEX array, dimension (LDX,NRHS) 00224 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to 00225 *> the original system of equations. Note that if EQUED = 'Y', 00226 *> A and B are modified on exit, and the solution to the 00227 *> equilibrated system is inv(diag(S))*X. 00228 *> \endverbatim 00229 *> 00230 *> \param[in] LDX 00231 *> \verbatim 00232 *> LDX is INTEGER 00233 *> The leading dimension of the array X. LDX >= max(1,N). 00234 *> \endverbatim 00235 *> 00236 *> \param[out] RCOND 00237 *> \verbatim 00238 *> RCOND is REAL 00239 *> The estimate of the reciprocal condition number of the matrix 00240 *> A after equilibration (if done). If RCOND is less than the 00241 *> machine precision (in particular, if RCOND = 0), the matrix 00242 *> is singular to working precision. This condition is 00243 *> indicated by a return code of INFO > 0. 00244 *> \endverbatim 00245 *> 00246 *> \param[out] FERR 00247 *> \verbatim 00248 *> FERR is REAL array, dimension (NRHS) 00249 *> The estimated forward error bound for each solution vector 00250 *> X(j) (the j-th column of the solution matrix X). 00251 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00252 *> is an estimated upper bound for the magnitude of the largest 00253 *> element in (X(j) - XTRUE) divided by the magnitude of the 00254 *> largest element in X(j). The estimate is as reliable as 00255 *> the estimate for RCOND, and is almost always a slight 00256 *> overestimate of the true error. 00257 *> \endverbatim 00258 *> 00259 *> \param[out] BERR 00260 *> \verbatim 00261 *> BERR is REAL array, dimension (NRHS) 00262 *> The componentwise relative backward error of each solution 00263 *> vector X(j) (i.e., the smallest relative change in 00264 *> any element of A or B that makes X(j) an exact solution). 00265 *> \endverbatim 00266 *> 00267 *> \param[out] WORK 00268 *> \verbatim 00269 *> WORK is COMPLEX array, dimension (2*N) 00270 *> \endverbatim 00271 *> 00272 *> \param[out] RWORK 00273 *> \verbatim 00274 *> RWORK is REAL array, dimension (N) 00275 *> \endverbatim 00276 *> 00277 *> \param[out] INFO 00278 *> \verbatim 00279 *> INFO is INTEGER 00280 *> = 0: successful exit 00281 *> < 0: if INFO = -i, the i-th argument had an illegal value 00282 *> > 0: if INFO = i, and i is 00283 *> <= N: the leading minor of order i of A is 00284 *> not positive definite, so the factorization 00285 *> could not be completed, and the solution has not 00286 *> been computed. RCOND = 0 is returned. 00287 *> = N+1: U is nonsingular, but RCOND is less than machine 00288 *> precision, meaning that the matrix is singular 00289 *> to working precision. Nevertheless, the 00290 *> solution and error bounds are computed because 00291 *> there are a number of situations where the 00292 *> computed solution can be more accurate than the 00293 *> value of RCOND would suggest. 00294 *> \endverbatim 00295 * 00296 * Authors: 00297 * ======== 00298 * 00299 *> \author Univ. of Tennessee 00300 *> \author Univ. of California Berkeley 00301 *> \author Univ. of Colorado Denver 00302 *> \author NAG Ltd. 00303 * 00304 *> \date April 2012 00305 * 00306 *> \ingroup complexOTHERsolve 00307 * 00308 *> \par Further Details: 00309 * ===================== 00310 *> 00311 *> \verbatim 00312 *> 00313 *> The band storage scheme is illustrated by the following example, when 00314 *> N = 6, KD = 2, and UPLO = 'U': 00315 *> 00316 *> Two-dimensional storage of the Hermitian matrix A: 00317 *> 00318 *> a11 a12 a13 00319 *> a22 a23 a24 00320 *> a33 a34 a35 00321 *> a44 a45 a46 00322 *> a55 a56 00323 *> (aij=conjg(aji)) a66 00324 *> 00325 *> Band storage of the upper triangle of A: 00326 *> 00327 *> * * a13 a24 a35 a46 00328 *> * a12 a23 a34 a45 a56 00329 *> a11 a22 a33 a44 a55 a66 00330 *> 00331 *> Similarly, if UPLO = 'L' the format of A is as follows: 00332 *> 00333 *> a11 a22 a33 a44 a55 a66 00334 *> a21 a32 a43 a54 a65 * 00335 *> a31 a42 a53 a64 * * 00336 *> 00337 *> Array elements marked * are not used by the routine. 00338 *> \endverbatim 00339 *> 00340 * ===================================================================== 00341 SUBROUTINE CPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, 00342 $ EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, 00343 $ WORK, RWORK, INFO ) 00344 * 00345 * -- LAPACK driver routine (version 3.4.1) -- 00346 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00347 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00348 * April 2012 00349 * 00350 * .. Scalar Arguments .. 00351 CHARACTER EQUED, FACT, UPLO 00352 INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS 00353 REAL RCOND 00354 * .. 00355 * .. Array Arguments .. 00356 REAL BERR( * ), FERR( * ), RWORK( * ), S( * ) 00357 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 00358 $ WORK( * ), X( LDX, * ) 00359 * .. 00360 * 00361 * ===================================================================== 00362 * 00363 * .. Parameters .. 00364 REAL ZERO, ONE 00365 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00366 * .. 00367 * .. Local Scalars .. 00368 LOGICAL EQUIL, NOFACT, RCEQU, UPPER 00369 INTEGER I, INFEQU, J, J1, J2 00370 REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 00371 * .. 00372 * .. External Functions .. 00373 LOGICAL LSAME 00374 REAL CLANHB, SLAMCH 00375 EXTERNAL LSAME, CLANHB, SLAMCH 00376 * .. 00377 * .. External Subroutines .. 00378 EXTERNAL CCOPY, CLACPY, CLAQHB, CPBCON, CPBEQU, CPBRFS, 00379 $ CPBTRF, CPBTRS, XERBLA 00380 * .. 00381 * .. Intrinsic Functions .. 00382 INTRINSIC MAX, MIN 00383 * .. 00384 * .. Executable Statements .. 00385 * 00386 INFO = 0 00387 NOFACT = LSAME( FACT, 'N' ) 00388 EQUIL = LSAME( FACT, 'E' ) 00389 UPPER = LSAME( UPLO, 'U' ) 00390 IF( NOFACT .OR. EQUIL ) THEN 00391 EQUED = 'N' 00392 RCEQU = .FALSE. 00393 ELSE 00394 RCEQU = LSAME( EQUED, 'Y' ) 00395 SMLNUM = SLAMCH( 'Safe minimum' ) 00396 BIGNUM = ONE / SMLNUM 00397 END IF 00398 * 00399 * Test the input parameters. 00400 * 00401 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 00402 $ THEN 00403 INFO = -1 00404 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00405 INFO = -2 00406 ELSE IF( N.LT.0 ) THEN 00407 INFO = -3 00408 ELSE IF( KD.LT.0 ) THEN 00409 INFO = -4 00410 ELSE IF( NRHS.LT.0 ) THEN 00411 INFO = -5 00412 ELSE IF( LDAB.LT.KD+1 ) THEN 00413 INFO = -7 00414 ELSE IF( LDAFB.LT.KD+1 ) THEN 00415 INFO = -9 00416 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 00417 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 00418 INFO = -10 00419 ELSE 00420 IF( RCEQU ) THEN 00421 SMIN = BIGNUM 00422 SMAX = ZERO 00423 DO 10 J = 1, N 00424 SMIN = MIN( SMIN, S( J ) ) 00425 SMAX = MAX( SMAX, S( J ) ) 00426 10 CONTINUE 00427 IF( SMIN.LE.ZERO ) THEN 00428 INFO = -11 00429 ELSE IF( N.GT.0 ) THEN 00430 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 00431 ELSE 00432 SCOND = ONE 00433 END IF 00434 END IF 00435 IF( INFO.EQ.0 ) THEN 00436 IF( LDB.LT.MAX( 1, N ) ) THEN 00437 INFO = -13 00438 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00439 INFO = -15 00440 END IF 00441 END IF 00442 END IF 00443 * 00444 IF( INFO.NE.0 ) THEN 00445 CALL XERBLA( 'CPBSVX', -INFO ) 00446 RETURN 00447 END IF 00448 * 00449 IF( EQUIL ) THEN 00450 * 00451 * Compute row and column scalings to equilibrate the matrix A. 00452 * 00453 CALL CPBEQU( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, INFEQU ) 00454 IF( INFEQU.EQ.0 ) THEN 00455 * 00456 * Equilibrate the matrix. 00457 * 00458 CALL CLAQHB( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED ) 00459 RCEQU = LSAME( EQUED, 'Y' ) 00460 END IF 00461 END IF 00462 * 00463 * Scale the right-hand side. 00464 * 00465 IF( RCEQU ) THEN 00466 DO 30 J = 1, NRHS 00467 DO 20 I = 1, N 00468 B( I, J ) = S( I )*B( I, J ) 00469 20 CONTINUE 00470 30 CONTINUE 00471 END IF 00472 * 00473 IF( NOFACT .OR. EQUIL ) THEN 00474 * 00475 * Compute the Cholesky factorization A = U**H *U or A = L*L**H. 00476 * 00477 IF( UPPER ) THEN 00478 DO 40 J = 1, N 00479 J1 = MAX( J-KD, 1 ) 00480 CALL CCOPY( J-J1+1, AB( KD+1-J+J1, J ), 1, 00481 $ AFB( KD+1-J+J1, J ), 1 ) 00482 40 CONTINUE 00483 ELSE 00484 DO 50 J = 1, N 00485 J2 = MIN( J+KD, N ) 00486 CALL CCOPY( J2-J+1, AB( 1, J ), 1, AFB( 1, J ), 1 ) 00487 50 CONTINUE 00488 END IF 00489 * 00490 CALL CPBTRF( UPLO, N, KD, AFB, LDAFB, INFO ) 00491 * 00492 * Return if INFO is non-zero. 00493 * 00494 IF( INFO.GT.0 )THEN 00495 RCOND = ZERO 00496 RETURN 00497 END IF 00498 END IF 00499 * 00500 * Compute the norm of the matrix A. 00501 * 00502 ANORM = CLANHB( '1', UPLO, N, KD, AB, LDAB, RWORK ) 00503 * 00504 * Compute the reciprocal of the condition number of A. 00505 * 00506 CALL CPBCON( UPLO, N, KD, AFB, LDAFB, ANORM, RCOND, WORK, RWORK, 00507 $ INFO ) 00508 * 00509 * Compute the solution matrix X. 00510 * 00511 CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00512 CALL CPBTRS( UPLO, N, KD, NRHS, AFB, LDAFB, X, LDX, INFO ) 00513 * 00514 * Use iterative refinement to improve the computed solution and 00515 * compute error bounds and backward error estimates for it. 00516 * 00517 CALL CPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B, LDB, X, 00518 $ LDX, FERR, BERR, WORK, RWORK, INFO ) 00519 * 00520 * Transform the solution matrix X to a solution of the original 00521 * system. 00522 * 00523 IF( RCEQU ) THEN 00524 DO 70 J = 1, NRHS 00525 DO 60 I = 1, N 00526 X( I, J ) = S( I )*X( I, J ) 00527 60 CONTINUE 00528 70 CONTINUE 00529 DO 80 J = 1, NRHS 00530 FERR( J ) = FERR( J ) / SCOND 00531 80 CONTINUE 00532 END IF 00533 * 00534 * Set INFO = N+1 if the matrix is singular to working precision. 00535 * 00536 IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 00537 $ INFO = N + 1 00538 * 00539 RETURN 00540 * 00541 * End of CPBSVX 00542 * 00543 END