001:       SUBROUTINE ZPTRFS( UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
002:      $                   FERR, BERR, WORK, RWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO
010:       INTEGER            INFO, LDB, LDX, N, NRHS
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   BERR( * ), D( * ), DF( * ), FERR( * ),
014:      $                   RWORK( * )
015:       COMPLEX*16         B( LDB, * ), E( * ), EF( * ), WORK( * ),
016:      $                   X( LDX, * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  ZPTRFS improves the computed solution to a system of linear
023: *  equations when the coefficient matrix is Hermitian positive definite
024: *  and tridiagonal, and provides error bounds and backward error
025: *  estimates for the solution.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  UPLO    (input) CHARACTER*1
031: *          Specifies whether the superdiagonal or the subdiagonal of the
032: *          tridiagonal matrix A is stored and the form of the
033: *          factorization:
034: *          = 'U':  E is the superdiagonal of A, and A = U**H*D*U;
035: *          = 'L':  E is the subdiagonal of A, and A = L*D*L**H.
036: *          (The two forms are equivalent if A is real.)
037: *
038: *  N       (input) INTEGER
039: *          The order of the matrix A.  N >= 0.
040: *
041: *  NRHS    (input) INTEGER
042: *          The number of right hand sides, i.e., the number of columns
043: *          of the matrix B.  NRHS >= 0.
044: *
045: *  D       (input) DOUBLE PRECISION array, dimension (N)
046: *          The n real diagonal elements of the tridiagonal matrix A.
047: *
048: *  E       (input) COMPLEX*16 array, dimension (N-1)
049: *          The (n-1) off-diagonal elements of the tridiagonal matrix A
050: *          (see UPLO).
051: *
052: *  DF      (input) DOUBLE PRECISION array, dimension (N)
053: *          The n diagonal elements of the diagonal matrix D from
054: *          the factorization computed by ZPTTRF.
055: *
056: *  EF      (input) COMPLEX*16 array, dimension (N-1)
057: *          The (n-1) off-diagonal elements of the unit bidiagonal
058: *          factor U or L from the factorization computed by ZPTTRF
059: *          (see UPLO).
060: *
061: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
062: *          The right hand side matrix B.
063: *
064: *  LDB     (input) INTEGER
065: *          The leading dimension of the array B.  LDB >= max(1,N).
066: *
067: *  X       (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
068: *          On entry, the solution matrix X, as computed by ZPTTRS.
069: *          On exit, the improved solution matrix X.
070: *
071: *  LDX     (input) INTEGER
072: *          The leading dimension of the array X.  LDX >= max(1,N).
073: *
074: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
075: *          The forward error bound for each solution vector
076: *          X(j) (the j-th column of the solution matrix X).
077: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
078: *          is an estimated upper bound for the magnitude of the largest
079: *          element in (X(j) - XTRUE) divided by the magnitude of the
080: *          largest element in X(j).
081: *
082: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
083: *          The componentwise relative backward error of each solution
084: *          vector X(j) (i.e., the smallest relative change in
085: *          any element of A or B that makes X(j) an exact solution).
086: *
087: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
088: *
089: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
090: *
091: *  INFO    (output) INTEGER
092: *          = 0:  successful exit
093: *          < 0:  if INFO = -i, the i-th argument had an illegal value
094: *
095: *  Internal Parameters
096: *  ===================
097: *
098: *  ITMAX is the maximum number of steps of iterative refinement.
099: *
100: *  =====================================================================
101: *
102: *     .. Parameters ..
103:       INTEGER            ITMAX
104:       PARAMETER          ( ITMAX = 5 )
105:       DOUBLE PRECISION   ZERO
106:       PARAMETER          ( ZERO = 0.0D+0 )
107:       DOUBLE PRECISION   ONE
108:       PARAMETER          ( ONE = 1.0D+0 )
109:       DOUBLE PRECISION   TWO
110:       PARAMETER          ( TWO = 2.0D+0 )
111:       DOUBLE PRECISION   THREE
112:       PARAMETER          ( THREE = 3.0D+0 )
113: *     ..
114: *     .. Local Scalars ..
115:       LOGICAL            UPPER
116:       INTEGER            COUNT, I, IX, J, NZ
117:       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
118:       COMPLEX*16         BI, CX, DX, EX, ZDUM
119: *     ..
120: *     .. External Functions ..
121:       LOGICAL            LSAME
122:       INTEGER            IDAMAX
123:       DOUBLE PRECISION   DLAMCH
124:       EXTERNAL           LSAME, IDAMAX, DLAMCH
125: *     ..
126: *     .. External Subroutines ..
127:       EXTERNAL           XERBLA, ZAXPY, ZPTTRS
128: *     ..
129: *     .. Intrinsic Functions ..
130:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
131: *     ..
132: *     .. Statement Functions ..
133:       DOUBLE PRECISION   CABS1
134: *     ..
135: *     .. Statement Function definitions ..
136:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
137: *     ..
138: *     .. Executable Statements ..
139: *
140: *     Test the input parameters.
141: *
142:       INFO = 0
143:       UPPER = LSAME( UPLO, 'U' )
144:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
145:          INFO = -1
146:       ELSE IF( N.LT.0 ) THEN
147:          INFO = -2
148:       ELSE IF( NRHS.LT.0 ) THEN
149:          INFO = -3
150:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
151:          INFO = -9
152:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
153:          INFO = -11
154:       END IF
155:       IF( INFO.NE.0 ) THEN
156:          CALL XERBLA( 'ZPTRFS', -INFO )
157:          RETURN
158:       END IF
159: *
160: *     Quick return if possible
161: *
162:       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
163:          DO 10 J = 1, NRHS
164:             FERR( J ) = ZERO
165:             BERR( J ) = ZERO
166:    10    CONTINUE
167:          RETURN
168:       END IF
169: *
170: *     NZ = maximum number of nonzero elements in each row of A, plus 1
171: *
172:       NZ = 4
173:       EPS = DLAMCH( 'Epsilon' )
174:       SAFMIN = DLAMCH( 'Safe minimum' )
175:       SAFE1 = NZ*SAFMIN
176:       SAFE2 = SAFE1 / EPS
177: *
178: *     Do for each right hand side
179: *
180:       DO 100 J = 1, NRHS
181: *
182:          COUNT = 1
183:          LSTRES = THREE
184:    20    CONTINUE
185: *
186: *        Loop until stopping criterion is satisfied.
187: *
188: *        Compute residual R = B - A * X.  Also compute
189: *        abs(A)*abs(x) + abs(b) for use in the backward error bound.
190: *
191:          IF( UPPER ) THEN
192:             IF( N.EQ.1 ) THEN
193:                BI = B( 1, J )
194:                DX = D( 1 )*X( 1, J )
195:                WORK( 1 ) = BI - DX
196:                RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
197:             ELSE
198:                BI = B( 1, J )
199:                DX = D( 1 )*X( 1, J )
200:                EX = E( 1 )*X( 2, J )
201:                WORK( 1 ) = BI - DX - EX
202:                RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
203:      $                      CABS1( E( 1 ) )*CABS1( X( 2, J ) )
204:                DO 30 I = 2, N - 1
205:                   BI = B( I, J )
206:                   CX = DCONJG( E( I-1 ) )*X( I-1, J )
207:                   DX = D( I )*X( I, J )
208:                   EX = E( I )*X( I+1, J )
209:                   WORK( I ) = BI - CX - DX - EX
210:                   RWORK( I ) = CABS1( BI ) +
211:      $                         CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
212:      $                         CABS1( DX ) + CABS1( E( I ) )*
213:      $                         CABS1( X( I+1, J ) )
214:    30          CONTINUE
215:                BI = B( N, J )
216:                CX = DCONJG( E( N-1 ) )*X( N-1, J )
217:                DX = D( N )*X( N, J )
218:                WORK( N ) = BI - CX - DX
219:                RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
220:      $                      CABS1( X( N-1, J ) ) + CABS1( DX )
221:             END IF
222:          ELSE
223:             IF( N.EQ.1 ) THEN
224:                BI = B( 1, J )
225:                DX = D( 1 )*X( 1, J )
226:                WORK( 1 ) = BI - DX
227:                RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
228:             ELSE
229:                BI = B( 1, J )
230:                DX = D( 1 )*X( 1, J )
231:                EX = DCONJG( E( 1 ) )*X( 2, J )
232:                WORK( 1 ) = BI - DX - EX
233:                RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
234:      $                      CABS1( E( 1 ) )*CABS1( X( 2, J ) )
235:                DO 40 I = 2, N - 1
236:                   BI = B( I, J )
237:                   CX = E( I-1 )*X( I-1, J )
238:                   DX = D( I )*X( I, J )
239:                   EX = DCONJG( E( I ) )*X( I+1, J )
240:                   WORK( I ) = BI - CX - DX - EX
241:                   RWORK( I ) = CABS1( BI ) +
242:      $                         CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
243:      $                         CABS1( DX ) + CABS1( E( I ) )*
244:      $                         CABS1( X( I+1, J ) )
245:    40          CONTINUE
246:                BI = B( N, J )
247:                CX = E( N-1 )*X( N-1, J )
248:                DX = D( N )*X( N, J )
249:                WORK( N ) = BI - CX - DX
250:                RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
251:      $                      CABS1( X( N-1, J ) ) + CABS1( DX )
252:             END IF
253:          END IF
254: *
255: *        Compute componentwise relative backward error from formula
256: *
257: *        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
258: *
259: *        where abs(Z) is the componentwise absolute value of the matrix
260: *        or vector Z.  If the i-th component of the denominator is less
261: *        than SAFE2, then SAFE1 is added to the i-th components of the
262: *        numerator and denominator before dividing.
263: *
264:          S = ZERO
265:          DO 50 I = 1, N
266:             IF( RWORK( I ).GT.SAFE2 ) THEN
267:                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
268:             ELSE
269:                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
270:      $             ( RWORK( I )+SAFE1 ) )
271:             END IF
272:    50    CONTINUE
273:          BERR( J ) = S
274: *
275: *        Test stopping criterion. Continue iterating if
276: *           1) The residual BERR(J) is larger than machine epsilon, and
277: *           2) BERR(J) decreased by at least a factor of 2 during the
278: *              last iteration, and
279: *           3) At most ITMAX iterations tried.
280: *
281:          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
282:      $       COUNT.LE.ITMAX ) THEN
283: *
284: *           Update solution and try again.
285: *
286:             CALL ZPTTRS( UPLO, N, 1, DF, EF, WORK, N, INFO )
287:             CALL ZAXPY( N, DCMPLX( ONE ), WORK, 1, X( 1, J ), 1 )
288:             LSTRES = BERR( J )
289:             COUNT = COUNT + 1
290:             GO TO 20
291:          END IF
292: *
293: *        Bound error from formula
294: *
295: *        norm(X - XTRUE) / norm(X) .le. FERR =
296: *        norm( abs(inv(A))*
297: *           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
298: *
299: *        where
300: *          norm(Z) is the magnitude of the largest component of Z
301: *          inv(A) is the inverse of A
302: *          abs(Z) is the componentwise absolute value of the matrix or
303: *             vector Z
304: *          NZ is the maximum number of nonzeros in any row of A, plus 1
305: *          EPS is machine epsilon
306: *
307: *        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
308: *        is incremented by SAFE1 if the i-th component of
309: *        abs(A)*abs(X) + abs(B) is less than SAFE2.
310: *
311:          DO 60 I = 1, N
312:             IF( RWORK( I ).GT.SAFE2 ) THEN
313:                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
314:             ELSE
315:                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
316:      $                      SAFE1
317:             END IF
318:    60    CONTINUE
319:          IX = IDAMAX( N, RWORK, 1 )
320:          FERR( J ) = RWORK( IX )
321: *
322: *        Estimate the norm of inv(A).
323: *
324: *        Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
325: *
326: *           m(i,j) =  abs(A(i,j)), i = j,
327: *           m(i,j) = -abs(A(i,j)), i .ne. j,
328: *
329: *        and e = [ 1, 1, ..., 1 ]'.  Note M(A) = M(L)*D*M(L)'.
330: *
331: *        Solve M(L) * x = e.
332: *
333:          RWORK( 1 ) = ONE
334:          DO 70 I = 2, N
335:             RWORK( I ) = ONE + RWORK( I-1 )*ABS( EF( I-1 ) )
336:    70    CONTINUE
337: *
338: *        Solve D * M(L)' * x = b.
339: *
340:          RWORK( N ) = RWORK( N ) / DF( N )
341:          DO 80 I = N - 1, 1, -1
342:             RWORK( I ) = RWORK( I ) / DF( I ) +
343:      $                   RWORK( I+1 )*ABS( EF( I ) )
344:    80    CONTINUE
345: *
346: *        Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
347: *
348:          IX = IDAMAX( N, RWORK, 1 )
349:          FERR( J ) = FERR( J )*ABS( RWORK( IX ) )
350: *
351: *        Normalize error.
352: *
353:          LSTRES = ZERO
354:          DO 90 I = 1, N
355:             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
356:    90    CONTINUE
357:          IF( LSTRES.NE.ZERO )
358:      $      FERR( J ) = FERR( J ) / LSTRES
359: *
360:   100 CONTINUE
361: *
362:       RETURN
363: *
364: *     End of ZPTRFS
365: *
366:       END
367: