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