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