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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, LDB, LDX, N, NRHS
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   BERR( * ), D( * ), DF( * ), FERR( * ),
015:      $                   RWORK( * )
016:       COMPLEX*16         B( LDB, * ), E( * ), EF( * ), WORK( * ),
017:      $                   X( LDX, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  ZPTRFS improves the computed solution to a system of linear
024: *  equations when the coefficient matrix is Hermitian positive definite
025: *  and tridiagonal, and provides error bounds and backward error
026: *  estimates for the solution.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  UPLO    (input) CHARACTER*1
032: *          Specifies whether the superdiagonal or the subdiagonal of the
033: *          tridiagonal matrix A is stored and the form of the
034: *          factorization:
035: *          = 'U':  E is the superdiagonal of A, and A = U**H*D*U;
036: *          = 'L':  E is the subdiagonal of A, and A = L*D*L**H.
037: *          (The two forms are equivalent if A is real.)
038: *
039: *  N       (input) INTEGER
040: *          The order of the matrix A.  N >= 0.
041: *
042: *  NRHS    (input) INTEGER
043: *          The number of right hand sides, i.e., the number of columns
044: *          of the matrix B.  NRHS >= 0.
045: *
046: *  D       (input) DOUBLE PRECISION array, dimension (N)
047: *          The n real diagonal elements of the tridiagonal matrix A.
048: *
049: *  E       (input) COMPLEX*16 array, dimension (N-1)
050: *          The (n-1) off-diagonal elements of the tridiagonal matrix A
051: *          (see UPLO).
052: *
053: *  DF      (input) DOUBLE PRECISION array, dimension (N)
054: *          The n diagonal elements of the diagonal matrix D from
055: *          the factorization computed by ZPTTRF.
056: *
057: *  EF      (input) COMPLEX*16 array, dimension (N-1)
058: *          The (n-1) off-diagonal elements of the unit bidiagonal
059: *          factor U or L from the factorization computed by ZPTTRF
060: *          (see UPLO).
061: *
062: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
063: *          The right hand side matrix B.
064: *
065: *  LDB     (input) INTEGER
066: *          The leading dimension of the array B.  LDB >= max(1,N).
067: *
068: *  X       (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
069: *          On entry, the solution matrix X, as computed by ZPTTRS.
070: *          On exit, the improved solution matrix X.
071: *
072: *  LDX     (input) INTEGER
073: *          The leading dimension of the array X.  LDX >= max(1,N).
074: *
075: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
076: *          The forward error bound for each solution vector
077: *          X(j) (the j-th column of the solution matrix X).
078: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
079: *          is an estimated upper bound for the magnitude of the largest
080: *          element in (X(j) - XTRUE) divided by the magnitude of the
081: *          largest element in X(j).
082: *
083: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
084: *          The componentwise relative backward error of each solution
085: *          vector X(j) (i.e., the smallest relative change in
086: *          any element of A or B that makes X(j) an exact solution).
087: *
088: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
089: *
090: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
091: *
092: *  INFO    (output) INTEGER
093: *          = 0:  successful exit
094: *          < 0:  if INFO = -i, the i-th argument had an illegal value
095: *
096: *  Internal Parameters
097: *  ===================
098: *
099: *  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
368: