Annotation of rpl/lapack/lapack/zptrfs.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZPTRFS( UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
! 2: $ FERR, BERR, WORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER UPLO
! 11: INTEGER INFO, LDB, LDX, N, NRHS
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
! 15: $ RWORK( * )
! 16: COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
! 17: $ X( LDX, * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * ZPTRFS improves the computed solution to a system of linear
! 24: * equations when the coefficient matrix is Hermitian positive definite
! 25: * and tridiagonal, and provides error bounds and backward error
! 26: * estimates for the solution.
! 27: *
! 28: * Arguments
! 29: * =========
! 30: *
! 31: * UPLO (input) CHARACTER*1
! 32: * Specifies whether the superdiagonal or the subdiagonal of the
! 33: * tridiagonal matrix A is stored and the form of the
! 34: * factorization:
! 35: * = 'U': E is the superdiagonal of A, and A = U**H*D*U;
! 36: * = 'L': E is the subdiagonal of A, and A = L*D*L**H.
! 37: * (The two forms are equivalent if A is real.)
! 38: *
! 39: * N (input) INTEGER
! 40: * The order of the matrix A. N >= 0.
! 41: *
! 42: * NRHS (input) INTEGER
! 43: * The number of right hand sides, i.e., the number of columns
! 44: * of the matrix B. NRHS >= 0.
! 45: *
! 46: * D (input) DOUBLE PRECISION array, dimension (N)
! 47: * The n real diagonal elements of the tridiagonal matrix A.
! 48: *
! 49: * E (input) COMPLEX*16 array, dimension (N-1)
! 50: * The (n-1) off-diagonal elements of the tridiagonal matrix A
! 51: * (see UPLO).
! 52: *
! 53: * DF (input) DOUBLE PRECISION array, dimension (N)
! 54: * The n diagonal elements of the diagonal matrix D from
! 55: * the factorization computed by ZPTTRF.
! 56: *
! 57: * EF (input) COMPLEX*16 array, dimension (N-1)
! 58: * The (n-1) off-diagonal elements of the unit bidiagonal
! 59: * factor U or L from the factorization computed by ZPTTRF
! 60: * (see UPLO).
! 61: *
! 62: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
! 63: * The right hand side matrix B.
! 64: *
! 65: * LDB (input) INTEGER
! 66: * The leading dimension of the array B. LDB >= max(1,N).
! 67: *
! 68: * X (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
! 69: * On entry, the solution matrix X, as computed by ZPTTRS.
! 70: * On exit, the improved solution matrix X.
! 71: *
! 72: * LDX (input) INTEGER
! 73: * The leading dimension of the array X. LDX >= max(1,N).
! 74: *
! 75: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
! 76: * The forward error bound for each solution vector
! 77: * X(j) (the j-th column of the solution matrix X).
! 78: * If XTRUE is the true solution corresponding to X(j), FERR(j)
! 79: * is an estimated upper bound for the magnitude of the largest
! 80: * element in (X(j) - XTRUE) divided by the magnitude of the
! 81: * largest element in X(j).
! 82: *
! 83: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
! 84: * The componentwise relative backward error of each solution
! 85: * vector X(j) (i.e., the smallest relative change in
! 86: * any element of A or B that makes X(j) an exact solution).
! 87: *
! 88: * WORK (workspace) COMPLEX*16 array, dimension (N)
! 89: *
! 90: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
! 91: *
! 92: * INFO (output) INTEGER
! 93: * = 0: successful exit
! 94: * < 0: if INFO = -i, the i-th argument had an illegal value
! 95: *
! 96: * Internal Parameters
! 97: * ===================
! 98: *
! 99: * ITMAX is the maximum number of steps of iterative refinement.
! 100: *
! 101: * =====================================================================
! 102: *
! 103: * .. Parameters ..
! 104: INTEGER ITMAX
! 105: PARAMETER ( ITMAX = 5 )
! 106: DOUBLE PRECISION ZERO
! 107: PARAMETER ( ZERO = 0.0D+0 )
! 108: DOUBLE PRECISION ONE
! 109: PARAMETER ( ONE = 1.0D+0 )
! 110: DOUBLE PRECISION TWO
! 111: PARAMETER ( TWO = 2.0D+0 )
! 112: DOUBLE PRECISION THREE
! 113: PARAMETER ( THREE = 3.0D+0 )
! 114: * ..
! 115: * .. Local Scalars ..
! 116: LOGICAL UPPER
! 117: INTEGER COUNT, I, IX, J, NZ
! 118: DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
! 119: COMPLEX*16 BI, CX, DX, EX, ZDUM
! 120: * ..
! 121: * .. External Functions ..
! 122: LOGICAL LSAME
! 123: INTEGER IDAMAX
! 124: DOUBLE PRECISION DLAMCH
! 125: EXTERNAL LSAME, IDAMAX, DLAMCH
! 126: * ..
! 127: * .. External Subroutines ..
! 128: EXTERNAL XERBLA, ZAXPY, ZPTTRS
! 129: * ..
! 130: * .. Intrinsic Functions ..
! 131: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
! 132: * ..
! 133: * .. Statement Functions ..
! 134: DOUBLE PRECISION CABS1
! 135: * ..
! 136: * .. Statement Function definitions ..
! 137: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
! 138: * ..
! 139: * .. Executable Statements ..
! 140: *
! 141: * Test the input parameters.
! 142: *
! 143: INFO = 0
! 144: UPPER = LSAME( UPLO, 'U' )
! 145: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 146: INFO = -1
! 147: ELSE IF( N.LT.0 ) THEN
! 148: INFO = -2
! 149: ELSE IF( NRHS.LT.0 ) THEN
! 150: INFO = -3
! 151: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 152: INFO = -9
! 153: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
! 154: INFO = -11
! 155: END IF
! 156: IF( INFO.NE.0 ) THEN
! 157: CALL XERBLA( 'ZPTRFS', -INFO )
! 158: RETURN
! 159: END IF
! 160: *
! 161: * Quick return if possible
! 162: *
! 163: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
! 164: DO 10 J = 1, NRHS
! 165: FERR( J ) = ZERO
! 166: BERR( J ) = ZERO
! 167: 10 CONTINUE
! 168: RETURN
! 169: END IF
! 170: *
! 171: * NZ = maximum number of nonzero elements in each row of A, plus 1
! 172: *
! 173: NZ = 4
! 174: EPS = DLAMCH( 'Epsilon' )
! 175: SAFMIN = DLAMCH( 'Safe minimum' )
! 176: SAFE1 = NZ*SAFMIN
! 177: SAFE2 = SAFE1 / EPS
! 178: *
! 179: * Do for each right hand side
! 180: *
! 181: DO 100 J = 1, NRHS
! 182: *
! 183: COUNT = 1
! 184: LSTRES = THREE
! 185: 20 CONTINUE
! 186: *
! 187: * Loop until stopping criterion is satisfied.
! 188: *
! 189: * Compute residual R = B - A * X. Also compute
! 190: * abs(A)*abs(x) + abs(b) for use in the backward error bound.
! 191: *
! 192: IF( UPPER ) THEN
! 193: IF( N.EQ.1 ) THEN
! 194: BI = B( 1, J )
! 195: DX = D( 1 )*X( 1, J )
! 196: WORK( 1 ) = BI - DX
! 197: RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
! 198: ELSE
! 199: BI = B( 1, J )
! 200: DX = D( 1 )*X( 1, J )
! 201: EX = E( 1 )*X( 2, J )
! 202: WORK( 1 ) = BI - DX - EX
! 203: RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
! 204: $ CABS1( E( 1 ) )*CABS1( X( 2, J ) )
! 205: DO 30 I = 2, N - 1
! 206: BI = B( I, J )
! 207: CX = DCONJG( E( I-1 ) )*X( I-1, J )
! 208: DX = D( I )*X( I, J )
! 209: EX = E( I )*X( I+1, J )
! 210: WORK( I ) = BI - CX - DX - EX
! 211: RWORK( I ) = CABS1( BI ) +
! 212: $ CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
! 213: $ CABS1( DX ) + CABS1( E( I ) )*
! 214: $ CABS1( X( I+1, J ) )
! 215: 30 CONTINUE
! 216: BI = B( N, J )
! 217: CX = DCONJG( E( N-1 ) )*X( N-1, J )
! 218: DX = D( N )*X( N, J )
! 219: WORK( N ) = BI - CX - DX
! 220: RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
! 221: $ CABS1( X( N-1, J ) ) + CABS1( DX )
! 222: END IF
! 223: ELSE
! 224: IF( N.EQ.1 ) THEN
! 225: BI = B( 1, J )
! 226: DX = D( 1 )*X( 1, J )
! 227: WORK( 1 ) = BI - DX
! 228: RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
! 229: ELSE
! 230: BI = B( 1, J )
! 231: DX = D( 1 )*X( 1, J )
! 232: EX = DCONJG( E( 1 ) )*X( 2, J )
! 233: WORK( 1 ) = BI - DX - EX
! 234: RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
! 235: $ CABS1( E( 1 ) )*CABS1( X( 2, J ) )
! 236: DO 40 I = 2, N - 1
! 237: BI = B( I, J )
! 238: CX = E( I-1 )*X( I-1, J )
! 239: DX = D( I )*X( I, J )
! 240: EX = DCONJG( E( I ) )*X( I+1, J )
! 241: WORK( I ) = BI - CX - DX - EX
! 242: RWORK( I ) = CABS1( BI ) +
! 243: $ CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
! 244: $ CABS1( DX ) + CABS1( E( I ) )*
! 245: $ CABS1( X( I+1, J ) )
! 246: 40 CONTINUE
! 247: BI = B( N, J )
! 248: CX = E( N-1 )*X( N-1, J )
! 249: DX = D( N )*X( N, J )
! 250: WORK( N ) = BI - CX - DX
! 251: RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
! 252: $ CABS1( X( N-1, J ) ) + CABS1( DX )
! 253: END IF
! 254: END IF
! 255: *
! 256: * Compute componentwise relative backward error from formula
! 257: *
! 258: * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
! 259: *
! 260: * where abs(Z) is the componentwise absolute value of the matrix
! 261: * or vector Z. If the i-th component of the denominator is less
! 262: * than SAFE2, then SAFE1 is added to the i-th components of the
! 263: * numerator and denominator before dividing.
! 264: *
! 265: S = ZERO
! 266: DO 50 I = 1, N
! 267: IF( RWORK( I ).GT.SAFE2 ) THEN
! 268: S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
! 269: ELSE
! 270: S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
! 271: $ ( RWORK( I )+SAFE1 ) )
! 272: END IF
! 273: 50 CONTINUE
! 274: BERR( J ) = S
! 275: *
! 276: * Test stopping criterion. Continue iterating if
! 277: * 1) The residual BERR(J) is larger than machine epsilon, and
! 278: * 2) BERR(J) decreased by at least a factor of 2 during the
! 279: * last iteration, and
! 280: * 3) At most ITMAX iterations tried.
! 281: *
! 282: IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
! 283: $ COUNT.LE.ITMAX ) THEN
! 284: *
! 285: * Update solution and try again.
! 286: *
! 287: CALL ZPTTRS( UPLO, N, 1, DF, EF, WORK, N, INFO )
! 288: CALL ZAXPY( N, DCMPLX( ONE ), WORK, 1, X( 1, J ), 1 )
! 289: LSTRES = BERR( J )
! 290: COUNT = COUNT + 1
! 291: GO TO 20
! 292: END IF
! 293: *
! 294: * Bound error from formula
! 295: *
! 296: * norm(X - XTRUE) / norm(X) .le. FERR =
! 297: * norm( abs(inv(A))*
! 298: * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
! 299: *
! 300: * where
! 301: * norm(Z) is the magnitude of the largest component of Z
! 302: * inv(A) is the inverse of A
! 303: * abs(Z) is the componentwise absolute value of the matrix or
! 304: * vector Z
! 305: * NZ is the maximum number of nonzeros in any row of A, plus 1
! 306: * EPS is machine epsilon
! 307: *
! 308: * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
! 309: * is incremented by SAFE1 if the i-th component of
! 310: * abs(A)*abs(X) + abs(B) is less than SAFE2.
! 311: *
! 312: DO 60 I = 1, N
! 313: IF( RWORK( I ).GT.SAFE2 ) THEN
! 314: RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
! 315: ELSE
! 316: RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
! 317: $ SAFE1
! 318: END IF
! 319: 60 CONTINUE
! 320: IX = IDAMAX( N, RWORK, 1 )
! 321: FERR( J ) = RWORK( IX )
! 322: *
! 323: * Estimate the norm of inv(A).
! 324: *
! 325: * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
! 326: *
! 327: * m(i,j) = abs(A(i,j)), i = j,
! 328: * m(i,j) = -abs(A(i,j)), i .ne. j,
! 329: *
! 330: * and e = [ 1, 1, ..., 1 ]'. Note M(A) = M(L)*D*M(L)'.
! 331: *
! 332: * Solve M(L) * x = e.
! 333: *
! 334: RWORK( 1 ) = ONE
! 335: DO 70 I = 2, N
! 336: RWORK( I ) = ONE + RWORK( I-1 )*ABS( EF( I-1 ) )
! 337: 70 CONTINUE
! 338: *
! 339: * Solve D * M(L)' * x = b.
! 340: *
! 341: RWORK( N ) = RWORK( N ) / DF( N )
! 342: DO 80 I = N - 1, 1, -1
! 343: RWORK( I ) = RWORK( I ) / DF( I ) +
! 344: $ RWORK( I+1 )*ABS( EF( I ) )
! 345: 80 CONTINUE
! 346: *
! 347: * Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
! 348: *
! 349: IX = IDAMAX( N, RWORK, 1 )
! 350: FERR( J ) = FERR( J )*ABS( RWORK( IX ) )
! 351: *
! 352: * Normalize error.
! 353: *
! 354: LSTRES = ZERO
! 355: DO 90 I = 1, N
! 356: LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
! 357: 90 CONTINUE
! 358: IF( LSTRES.NE.ZERO )
! 359: $ FERR( J ) = FERR( J ) / LSTRES
! 360: *
! 361: 100 CONTINUE
! 362: *
! 363: RETURN
! 364: *
! 365: * End of ZPTRFS
! 366: *
! 367: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zptrfs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack