1: *> \brief \b DTPRFS
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTPRFS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtprfs.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtprfs.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtprfs.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
22: * FERR, BERR, WORK, IWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER DIAG, TRANS, UPLO
26: * INTEGER INFO, LDB, LDX, N, NRHS
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IWORK( * )
30: * DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
31: * $ WORK( * ), X( LDX, * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DTPRFS provides error bounds and backward error estimates for the
41: *> solution to a system of linear equations with a triangular packed
42: *> coefficient matrix.
43: *>
44: *> The solution matrix X must be computed by DTPTRS or some other
45: *> means before entering this routine. DTPRFS does not do iterative
46: *> refinement because doing so cannot improve the backward error.
47: *> \endverbatim
48: *
49: * Arguments:
50: * ==========
51: *
52: *> \param[in] UPLO
53: *> \verbatim
54: *> UPLO is CHARACTER*1
55: *> = 'U': A is upper triangular;
56: *> = 'L': A is lower triangular.
57: *> \endverbatim
58: *>
59: *> \param[in] TRANS
60: *> \verbatim
61: *> TRANS is CHARACTER*1
62: *> Specifies the form of the system of equations:
63: *> = 'N': A * X = B (No transpose)
64: *> = 'T': A**T * X = B (Transpose)
65: *> = 'C': A**H * X = B (Conjugate transpose = Transpose)
66: *> \endverbatim
67: *>
68: *> \param[in] DIAG
69: *> \verbatim
70: *> DIAG is CHARACTER*1
71: *> = 'N': A is non-unit triangular;
72: *> = 'U': A is unit triangular.
73: *> \endverbatim
74: *>
75: *> \param[in] N
76: *> \verbatim
77: *> N is INTEGER
78: *> The order of the matrix A. N >= 0.
79: *> \endverbatim
80: *>
81: *> \param[in] NRHS
82: *> \verbatim
83: *> NRHS is INTEGER
84: *> The number of right hand sides, i.e., the number of columns
85: *> of the matrices B and X. NRHS >= 0.
86: *> \endverbatim
87: *>
88: *> \param[in] AP
89: *> \verbatim
90: *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
91: *> The upper or lower triangular matrix A, packed columnwise in
92: *> a linear array. The j-th column of A is stored in the array
93: *> AP as follows:
94: *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
95: *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
96: *> If DIAG = 'U', the diagonal elements of A are not referenced
97: *> and are assumed to be 1.
98: *> \endverbatim
99: *>
100: *> \param[in] B
101: *> \verbatim
102: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
103: *> The right hand side matrix B.
104: *> \endverbatim
105: *>
106: *> \param[in] LDB
107: *> \verbatim
108: *> LDB is INTEGER
109: *> The leading dimension of the array B. LDB >= max(1,N).
110: *> \endverbatim
111: *>
112: *> \param[in] X
113: *> \verbatim
114: *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
115: *> The solution matrix X.
116: *> \endverbatim
117: *>
118: *> \param[in] LDX
119: *> \verbatim
120: *> LDX is INTEGER
121: *> The leading dimension of the array X. LDX >= max(1,N).
122: *> \endverbatim
123: *>
124: *> \param[out] FERR
125: *> \verbatim
126: *> FERR is DOUBLE PRECISION array, dimension (NRHS)
127: *> The estimated forward error bound for each solution vector
128: *> X(j) (the j-th column of the solution matrix X).
129: *> If XTRUE is the true solution corresponding to X(j), FERR(j)
130: *> is an estimated upper bound for the magnitude of the largest
131: *> element in (X(j) - XTRUE) divided by the magnitude of the
132: *> largest element in X(j). The estimate is as reliable as
133: *> the estimate for RCOND, and is almost always a slight
134: *> overestimate of the true error.
135: *> \endverbatim
136: *>
137: *> \param[out] BERR
138: *> \verbatim
139: *> BERR is DOUBLE PRECISION array, dimension (NRHS)
140: *> The componentwise relative backward error of each solution
141: *> vector X(j) (i.e., the smallest relative change in
142: *> any element of A or B that makes X(j) an exact solution).
143: *> \endverbatim
144: *>
145: *> \param[out] WORK
146: *> \verbatim
147: *> WORK is DOUBLE PRECISION array, dimension (3*N)
148: *> \endverbatim
149: *>
150: *> \param[out] IWORK
151: *> \verbatim
152: *> IWORK is INTEGER array, dimension (N)
153: *> \endverbatim
154: *>
155: *> \param[out] INFO
156: *> \verbatim
157: *> INFO is INTEGER
158: *> = 0: successful exit
159: *> < 0: if INFO = -i, the i-th argument had an illegal value
160: *> \endverbatim
161: *
162: * Authors:
163: * ========
164: *
165: *> \author Univ. of Tennessee
166: *> \author Univ. of California Berkeley
167: *> \author Univ. of Colorado Denver
168: *> \author NAG Ltd.
169: *
170: *> \ingroup doubleOTHERcomputational
171: *
172: * =====================================================================
173: SUBROUTINE DTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
174: $ FERR, BERR, WORK, IWORK, INFO )
175: *
176: * -- LAPACK computational routine --
177: * -- LAPACK is a software package provided by Univ. of Tennessee, --
178: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179: *
180: * .. Scalar Arguments ..
181: CHARACTER DIAG, TRANS, UPLO
182: INTEGER INFO, LDB, LDX, N, NRHS
183: * ..
184: * .. Array Arguments ..
185: INTEGER IWORK( * )
186: DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
187: $ WORK( * ), X( LDX, * )
188: * ..
189: *
190: * =====================================================================
191: *
192: * .. Parameters ..
193: DOUBLE PRECISION ZERO
194: PARAMETER ( ZERO = 0.0D+0 )
195: DOUBLE PRECISION ONE
196: PARAMETER ( ONE = 1.0D+0 )
197: * ..
198: * .. Local Scalars ..
199: LOGICAL NOTRAN, NOUNIT, UPPER
200: CHARACTER TRANST
201: INTEGER I, J, K, KASE, KC, NZ
202: DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
203: * ..
204: * .. Local Arrays ..
205: INTEGER ISAVE( 3 )
206: * ..
207: * .. External Subroutines ..
208: EXTERNAL DAXPY, DCOPY, DLACN2, DTPMV, DTPSV, XERBLA
209: * ..
210: * .. Intrinsic Functions ..
211: INTRINSIC ABS, MAX
212: * ..
213: * .. External Functions ..
214: LOGICAL LSAME
215: DOUBLE PRECISION DLAMCH
216: EXTERNAL LSAME, DLAMCH
217: * ..
218: * .. Executable Statements ..
219: *
220: * Test the input parameters.
221: *
222: INFO = 0
223: UPPER = LSAME( UPLO, 'U' )
224: NOTRAN = LSAME( TRANS, 'N' )
225: NOUNIT = LSAME( DIAG, 'N' )
226: *
227: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
228: INFO = -1
229: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
230: $ LSAME( TRANS, 'C' ) ) THEN
231: INFO = -2
232: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
233: INFO = -3
234: ELSE IF( N.LT.0 ) THEN
235: INFO = -4
236: ELSE IF( NRHS.LT.0 ) THEN
237: INFO = -5
238: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
239: INFO = -8
240: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
241: INFO = -10
242: END IF
243: IF( INFO.NE.0 ) THEN
244: CALL XERBLA( 'DTPRFS', -INFO )
245: RETURN
246: END IF
247: *
248: * Quick return if possible
249: *
250: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
251: DO 10 J = 1, NRHS
252: FERR( J ) = ZERO
253: BERR( J ) = ZERO
254: 10 CONTINUE
255: RETURN
256: END IF
257: *
258: IF( NOTRAN ) THEN
259: TRANST = 'T'
260: ELSE
261: TRANST = 'N'
262: END IF
263: *
264: * NZ = maximum number of nonzero elements in each row of A, plus 1
265: *
266: NZ = N + 1
267: EPS = DLAMCH( 'Epsilon' )
268: SAFMIN = DLAMCH( 'Safe minimum' )
269: SAFE1 = NZ*SAFMIN
270: SAFE2 = SAFE1 / EPS
271: *
272: * Do for each right hand side
273: *
274: DO 250 J = 1, NRHS
275: *
276: * Compute residual R = B - op(A) * X,
277: * where op(A) = A or A**T, depending on TRANS.
278: *
279: CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 )
280: CALL DTPMV( UPLO, TRANS, DIAG, N, AP, WORK( N+1 ), 1 )
281: CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 )
282: *
283: * Compute componentwise relative backward error from formula
284: *
285: * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
286: *
287: * where abs(Z) is the componentwise absolute value of the matrix
288: * or vector Z. If the i-th component of the denominator is less
289: * than SAFE2, then SAFE1 is added to the i-th components of the
290: * numerator and denominator before dividing.
291: *
292: DO 20 I = 1, N
293: WORK( I ) = ABS( B( I, J ) )
294: 20 CONTINUE
295: *
296: IF( NOTRAN ) THEN
297: *
298: * Compute abs(A)*abs(X) + abs(B).
299: *
300: IF( UPPER ) THEN
301: KC = 1
302: IF( NOUNIT ) THEN
303: DO 40 K = 1, N
304: XK = ABS( X( K, J ) )
305: DO 30 I = 1, K
306: WORK( I ) = WORK( I ) + ABS( AP( KC+I-1 ) )*XK
307: 30 CONTINUE
308: KC = KC + K
309: 40 CONTINUE
310: ELSE
311: DO 60 K = 1, N
312: XK = ABS( X( K, J ) )
313: DO 50 I = 1, K - 1
314: WORK( I ) = WORK( I ) + ABS( AP( KC+I-1 ) )*XK
315: 50 CONTINUE
316: WORK( K ) = WORK( K ) + XK
317: KC = KC + K
318: 60 CONTINUE
319: END IF
320: ELSE
321: KC = 1
322: IF( NOUNIT ) THEN
323: DO 80 K = 1, N
324: XK = ABS( X( K, J ) )
325: DO 70 I = K, N
326: WORK( I ) = WORK( I ) + ABS( AP( KC+I-K ) )*XK
327: 70 CONTINUE
328: KC = KC + N - K + 1
329: 80 CONTINUE
330: ELSE
331: DO 100 K = 1, N
332: XK = ABS( X( K, J ) )
333: DO 90 I = K + 1, N
334: WORK( I ) = WORK( I ) + ABS( AP( KC+I-K ) )*XK
335: 90 CONTINUE
336: WORK( K ) = WORK( K ) + XK
337: KC = KC + N - K + 1
338: 100 CONTINUE
339: END IF
340: END IF
341: ELSE
342: *
343: * Compute abs(A**T)*abs(X) + abs(B).
344: *
345: IF( UPPER ) THEN
346: KC = 1
347: IF( NOUNIT ) THEN
348: DO 120 K = 1, N
349: S = ZERO
350: DO 110 I = 1, K
351: S = S + ABS( AP( KC+I-1 ) )*ABS( X( I, J ) )
352: 110 CONTINUE
353: WORK( K ) = WORK( K ) + S
354: KC = KC + K
355: 120 CONTINUE
356: ELSE
357: DO 140 K = 1, N
358: S = ABS( X( K, J ) )
359: DO 130 I = 1, K - 1
360: S = S + ABS( AP( KC+I-1 ) )*ABS( X( I, J ) )
361: 130 CONTINUE
362: WORK( K ) = WORK( K ) + S
363: KC = KC + K
364: 140 CONTINUE
365: END IF
366: ELSE
367: KC = 1
368: IF( NOUNIT ) THEN
369: DO 160 K = 1, N
370: S = ZERO
371: DO 150 I = K, N
372: S = S + ABS( AP( KC+I-K ) )*ABS( X( I, J ) )
373: 150 CONTINUE
374: WORK( K ) = WORK( K ) + S
375: KC = KC + N - K + 1
376: 160 CONTINUE
377: ELSE
378: DO 180 K = 1, N
379: S = ABS( X( K, J ) )
380: DO 170 I = K + 1, N
381: S = S + ABS( AP( KC+I-K ) )*ABS( X( I, J ) )
382: 170 CONTINUE
383: WORK( K ) = WORK( K ) + S
384: KC = KC + N - K + 1
385: 180 CONTINUE
386: END IF
387: END IF
388: END IF
389: S = ZERO
390: DO 190 I = 1, N
391: IF( WORK( I ).GT.SAFE2 ) THEN
392: S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
393: ELSE
394: S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
395: $ ( WORK( I )+SAFE1 ) )
396: END IF
397: 190 CONTINUE
398: BERR( J ) = S
399: *
400: * Bound error from formula
401: *
402: * norm(X - XTRUE) / norm(X) .le. FERR =
403: * norm( abs(inv(op(A)))*
404: * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
405: *
406: * where
407: * norm(Z) is the magnitude of the largest component of Z
408: * inv(op(A)) is the inverse of op(A)
409: * abs(Z) is the componentwise absolute value of the matrix or
410: * vector Z
411: * NZ is the maximum number of nonzeros in any row of A, plus 1
412: * EPS is machine epsilon
413: *
414: * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
415: * is incremented by SAFE1 if the i-th component of
416: * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
417: *
418: * Use DLACN2 to estimate the infinity-norm of the matrix
419: * inv(op(A)) * diag(W),
420: * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
421: *
422: DO 200 I = 1, N
423: IF( WORK( I ).GT.SAFE2 ) THEN
424: WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
425: ELSE
426: WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
427: END IF
428: 200 CONTINUE
429: *
430: KASE = 0
431: 210 CONTINUE
432: CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
433: $ KASE, ISAVE )
434: IF( KASE.NE.0 ) THEN
435: IF( KASE.EQ.1 ) THEN
436: *
437: * Multiply by diag(W)*inv(op(A)**T).
438: *
439: CALL DTPSV( UPLO, TRANST, DIAG, N, AP, WORK( N+1 ), 1 )
440: DO 220 I = 1, N
441: WORK( N+I ) = WORK( I )*WORK( N+I )
442: 220 CONTINUE
443: ELSE
444: *
445: * Multiply by inv(op(A))*diag(W).
446: *
447: DO 230 I = 1, N
448: WORK( N+I ) = WORK( I )*WORK( N+I )
449: 230 CONTINUE
450: CALL DTPSV( UPLO, TRANS, DIAG, N, AP, WORK( N+1 ), 1 )
451: END IF
452: GO TO 210
453: END IF
454: *
455: * Normalize error.
456: *
457: LSTRES = ZERO
458: DO 240 I = 1, N
459: LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
460: 240 CONTINUE
461: IF( LSTRES.NE.ZERO )
462: $ FERR( J ) = FERR( J ) / LSTRES
463: *
464: 250 CONTINUE
465: *
466: RETURN
467: *
468: * End of DTPRFS
469: *
470: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtprfs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack