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