![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SPTRFS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SPTRFS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sptrfs.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sptrfs.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sptrfs.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, 00022 * BERR, WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDB, LDX, N, NRHS 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL B( LDB, * ), BERR( * ), D( * ), DF( * ), 00029 * $ E( * ), EF( * ), FERR( * ), WORK( * ), 00030 * $ X( LDX, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SPTRFS improves the computed solution to a system of linear 00040 *> equations when the coefficient matrix is symmetric positive definite 00041 *> and tridiagonal, and provides error bounds and backward error 00042 *> estimates for the solution. 00043 *> \endverbatim 00044 * 00045 * Arguments: 00046 * ========== 00047 * 00048 *> \param[in] N 00049 *> \verbatim 00050 *> N is INTEGER 00051 *> The order of the matrix A. N >= 0. 00052 *> \endverbatim 00053 *> 00054 *> \param[in] NRHS 00055 *> \verbatim 00056 *> NRHS is INTEGER 00057 *> The number of right hand sides, i.e., the number of columns 00058 *> of the matrix B. NRHS >= 0. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] D 00062 *> \verbatim 00063 *> D is REAL array, dimension (N) 00064 *> The n diagonal elements of the tridiagonal matrix A. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] E 00068 *> \verbatim 00069 *> E is REAL array, dimension (N-1) 00070 *> The (n-1) subdiagonal elements of the tridiagonal matrix A. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] DF 00074 *> \verbatim 00075 *> DF is REAL array, dimension (N) 00076 *> The n diagonal elements of the diagonal matrix D from the 00077 *> factorization computed by SPTTRF. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] EF 00081 *> \verbatim 00082 *> EF is REAL array, dimension (N-1) 00083 *> The (n-1) subdiagonal elements of the unit bidiagonal factor 00084 *> L from the factorization computed by SPTTRF. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] B 00088 *> \verbatim 00089 *> B is REAL array, dimension (LDB,NRHS) 00090 *> The right hand side matrix B. 00091 *> \endverbatim 00092 *> 00093 *> \param[in] LDB 00094 *> \verbatim 00095 *> LDB is INTEGER 00096 *> The leading dimension of the array B. LDB >= max(1,N). 00097 *> \endverbatim 00098 *> 00099 *> \param[in,out] X 00100 *> \verbatim 00101 *> X is REAL array, dimension (LDX,NRHS) 00102 *> On entry, the solution matrix X, as computed by SPTTRS. 00103 *> On exit, the improved solution matrix X. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] LDX 00107 *> \verbatim 00108 *> LDX is INTEGER 00109 *> The leading dimension of the array X. LDX >= max(1,N). 00110 *> \endverbatim 00111 *> 00112 *> \param[out] FERR 00113 *> \verbatim 00114 *> FERR is REAL array, dimension (NRHS) 00115 *> The forward error bound for each solution vector 00116 *> X(j) (the j-th column of the solution matrix X). 00117 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00118 *> is an estimated upper bound for the magnitude of the largest 00119 *> element in (X(j) - XTRUE) divided by the magnitude of the 00120 *> largest element in X(j). 00121 *> \endverbatim 00122 *> 00123 *> \param[out] BERR 00124 *> \verbatim 00125 *> BERR is REAL array, dimension (NRHS) 00126 *> The componentwise relative backward error of each solution 00127 *> vector X(j) (i.e., the smallest relative change in 00128 *> any element of A or B that makes X(j) an exact solution). 00129 *> \endverbatim 00130 *> 00131 *> \param[out] WORK 00132 *> \verbatim 00133 *> WORK is REAL array, dimension (2*N) 00134 *> \endverbatim 00135 *> 00136 *> \param[out] INFO 00137 *> \verbatim 00138 *> INFO is INTEGER 00139 *> = 0: successful exit 00140 *> < 0: if INFO = -i, the i-th argument had an illegal value 00141 *> \endverbatim 00142 * 00143 *> \par Internal Parameters: 00144 * ========================= 00145 *> 00146 *> \verbatim 00147 *> ITMAX is the maximum number of steps of iterative refinement. 00148 *> \endverbatim 00149 * 00150 * Authors: 00151 * ======== 00152 * 00153 *> \author Univ. of Tennessee 00154 *> \author Univ. of California Berkeley 00155 *> \author Univ. of Colorado Denver 00156 *> \author NAG Ltd. 00157 * 00158 *> \date November 2011 00159 * 00160 *> \ingroup realOTHERcomputational 00161 * 00162 * ===================================================================== 00163 SUBROUTINE SPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, 00164 $ BERR, WORK, INFO ) 00165 * 00166 * -- LAPACK computational routine (version 3.4.0) -- 00167 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00169 * November 2011 00170 * 00171 * .. Scalar Arguments .. 00172 INTEGER INFO, LDB, LDX, N, NRHS 00173 * .. 00174 * .. Array Arguments .. 00175 REAL B( LDB, * ), BERR( * ), D( * ), DF( * ), 00176 $ E( * ), EF( * ), FERR( * ), WORK( * ), 00177 $ X( LDX, * ) 00178 * .. 00179 * 00180 * ===================================================================== 00181 * 00182 * .. Parameters .. 00183 INTEGER ITMAX 00184 PARAMETER ( ITMAX = 5 ) 00185 REAL ZERO 00186 PARAMETER ( ZERO = 0.0E+0 ) 00187 REAL ONE 00188 PARAMETER ( ONE = 1.0E+0 ) 00189 REAL TWO 00190 PARAMETER ( TWO = 2.0E+0 ) 00191 REAL THREE 00192 PARAMETER ( THREE = 3.0E+0 ) 00193 * .. 00194 * .. Local Scalars .. 00195 INTEGER COUNT, I, IX, J, NZ 00196 REAL BI, CX, DX, EPS, EX, LSTRES, S, SAFE1, SAFE2, 00197 $ SAFMIN 00198 * .. 00199 * .. External Subroutines .. 00200 EXTERNAL SAXPY, SPTTRS, XERBLA 00201 * .. 00202 * .. Intrinsic Functions .. 00203 INTRINSIC ABS, MAX 00204 * .. 00205 * .. External Functions .. 00206 INTEGER ISAMAX 00207 REAL SLAMCH 00208 EXTERNAL ISAMAX, SLAMCH 00209 * .. 00210 * .. Executable Statements .. 00211 * 00212 * Test the input parameters. 00213 * 00214 INFO = 0 00215 IF( N.LT.0 ) THEN 00216 INFO = -1 00217 ELSE IF( NRHS.LT.0 ) THEN 00218 INFO = -2 00219 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00220 INFO = -8 00221 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00222 INFO = -10 00223 END IF 00224 IF( INFO.NE.0 ) THEN 00225 CALL XERBLA( 'SPTRFS', -INFO ) 00226 RETURN 00227 END IF 00228 * 00229 * Quick return if possible 00230 * 00231 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 00232 DO 10 J = 1, NRHS 00233 FERR( J ) = ZERO 00234 BERR( J ) = ZERO 00235 10 CONTINUE 00236 RETURN 00237 END IF 00238 * 00239 * NZ = maximum number of nonzero elements in each row of A, plus 1 00240 * 00241 NZ = 4 00242 EPS = SLAMCH( 'Epsilon' ) 00243 SAFMIN = SLAMCH( 'Safe minimum' ) 00244 SAFE1 = NZ*SAFMIN 00245 SAFE2 = SAFE1 / EPS 00246 * 00247 * Do for each right hand side 00248 * 00249 DO 90 J = 1, NRHS 00250 * 00251 COUNT = 1 00252 LSTRES = THREE 00253 20 CONTINUE 00254 * 00255 * Loop until stopping criterion is satisfied. 00256 * 00257 * Compute residual R = B - A * X. Also compute 00258 * abs(A)*abs(x) + abs(b) for use in the backward error bound. 00259 * 00260 IF( N.EQ.1 ) THEN 00261 BI = B( 1, J ) 00262 DX = D( 1 )*X( 1, J ) 00263 WORK( N+1 ) = BI - DX 00264 WORK( 1 ) = ABS( BI ) + ABS( DX ) 00265 ELSE 00266 BI = B( 1, J ) 00267 DX = D( 1 )*X( 1, J ) 00268 EX = E( 1 )*X( 2, J ) 00269 WORK( N+1 ) = BI - DX - EX 00270 WORK( 1 ) = ABS( BI ) + ABS( DX ) + ABS( EX ) 00271 DO 30 I = 2, N - 1 00272 BI = B( I, J ) 00273 CX = E( I-1 )*X( I-1, J ) 00274 DX = D( I )*X( I, J ) 00275 EX = E( I )*X( I+1, J ) 00276 WORK( N+I ) = BI - CX - DX - EX 00277 WORK( I ) = ABS( BI ) + ABS( CX ) + ABS( DX ) + ABS( EX ) 00278 30 CONTINUE 00279 BI = B( N, J ) 00280 CX = E( N-1 )*X( N-1, J ) 00281 DX = D( N )*X( N, J ) 00282 WORK( N+N ) = BI - CX - DX 00283 WORK( N ) = ABS( BI ) + ABS( CX ) + ABS( DX ) 00284 END IF 00285 * 00286 * Compute componentwise relative backward error from formula 00287 * 00288 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) ) 00289 * 00290 * where abs(Z) is the componentwise absolute value of the matrix 00291 * or vector Z. If the i-th component of the denominator is less 00292 * than SAFE2, then SAFE1 is added to the i-th components of the 00293 * numerator and denominator before dividing. 00294 * 00295 S = ZERO 00296 DO 40 I = 1, N 00297 IF( WORK( I ).GT.SAFE2 ) THEN 00298 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) ) 00299 ELSE 00300 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) / 00301 $ ( WORK( I )+SAFE1 ) ) 00302 END IF 00303 40 CONTINUE 00304 BERR( J ) = S 00305 * 00306 * Test stopping criterion. Continue iterating if 00307 * 1) The residual BERR(J) is larger than machine epsilon, and 00308 * 2) BERR(J) decreased by at least a factor of 2 during the 00309 * last iteration, and 00310 * 3) At most ITMAX iterations tried. 00311 * 00312 IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND. 00313 $ COUNT.LE.ITMAX ) THEN 00314 * 00315 * Update solution and try again. 00316 * 00317 CALL SPTTRS( N, 1, DF, EF, WORK( N+1 ), N, INFO ) 00318 CALL SAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 ) 00319 LSTRES = BERR( J ) 00320 COUNT = COUNT + 1 00321 GO TO 20 00322 END IF 00323 * 00324 * Bound error from formula 00325 * 00326 * norm(X - XTRUE) / norm(X) .le. FERR = 00327 * norm( abs(inv(A))* 00328 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X) 00329 * 00330 * where 00331 * norm(Z) is the magnitude of the largest component of Z 00332 * inv(A) is the inverse of A 00333 * abs(Z) is the componentwise absolute value of the matrix or 00334 * vector Z 00335 * NZ is the maximum number of nonzeros in any row of A, plus 1 00336 * EPS is machine epsilon 00337 * 00338 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B)) 00339 * is incremented by SAFE1 if the i-th component of 00340 * abs(A)*abs(X) + abs(B) is less than SAFE2. 00341 * 00342 DO 50 I = 1, N 00343 IF( WORK( I ).GT.SAFE2 ) THEN 00344 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) 00345 ELSE 00346 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1 00347 END IF 00348 50 CONTINUE 00349 IX = ISAMAX( N, WORK, 1 ) 00350 FERR( J ) = WORK( IX ) 00351 * 00352 * Estimate the norm of inv(A). 00353 * 00354 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by 00355 * 00356 * m(i,j) = abs(A(i,j)), i = j, 00357 * m(i,j) = -abs(A(i,j)), i .ne. j, 00358 * 00359 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T. 00360 * 00361 * Solve M(L) * x = e. 00362 * 00363 WORK( 1 ) = ONE 00364 DO 60 I = 2, N 00365 WORK( I ) = ONE + WORK( I-1 )*ABS( EF( I-1 ) ) 00366 60 CONTINUE 00367 * 00368 * Solve D * M(L)**T * x = b. 00369 * 00370 WORK( N ) = WORK( N ) / DF( N ) 00371 DO 70 I = N - 1, 1, -1 00372 WORK( I ) = WORK( I ) / DF( I ) + WORK( I+1 )*ABS( EF( I ) ) 00373 70 CONTINUE 00374 * 00375 * Compute norm(inv(A)) = max(x(i)), 1<=i<=n. 00376 * 00377 IX = ISAMAX( N, WORK, 1 ) 00378 FERR( J ) = FERR( J )*ABS( WORK( IX ) ) 00379 * 00380 * Normalize error. 00381 * 00382 LSTRES = ZERO 00383 DO 80 I = 1, N 00384 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) ) 00385 80 CONTINUE 00386 IF( LSTRES.NE.ZERO ) 00387 $ FERR( J ) = FERR( J ) / LSTRES 00388 * 00389 90 CONTINUE 00390 * 00391 RETURN 00392 * 00393 * End of SPTRFS 00394 * 00395 END