LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtrrfs.f
Go to the documentation of this file.
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRSV
Definition: dtrsv.f:143
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
Definition: dtrmv.f:147
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:136
subroutine dtrrfs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
DTRRFS
Definition: dtrrfs.f:182