LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dtrrfs()

 subroutine dtrrfs ( character uplo, character trans, character diag, integer n, integer nrhs, double precision, dimension( lda, * ) a, integer lda, double precision, dimension( ldb, * ) b, integer ldb, double precision, dimension( ldx, * ) x, integer ldx, double precision, dimension( * ) ferr, double precision, dimension( * ) berr, double precision, dimension( * ) work, integer, dimension( * ) iwork, integer info )

DTRRFS

Purpose:
``` DTRRFS provides error bounds and backward error estimates for the
solution to a system of linear equations with a triangular
coefficient matrix.

The solution matrix X must be computed by DTRTRS or some other
means before entering this routine.  DTRRFS does not do iterative
refinement because doing so cannot improve the backward error.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': A is upper triangular; = 'L': A is lower triangular.``` [in] TRANS ``` TRANS is CHARACTER*1 Specifies the form of the system of equations: = 'N': A * X = B (No transpose) = 'T': A**T * X = B (Transpose) = 'C': A**H * X = B (Conjugate transpose = Transpose)``` [in] DIAG ``` DIAG is CHARACTER*1 = 'N': A is non-unit triangular; = 'U': A is unit triangular.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] NRHS ``` NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0.``` [in] A ``` A is DOUBLE PRECISION array, dimension (LDA,N) The triangular matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [in] B ``` B is DOUBLE PRECISION array, dimension (LDB,NRHS) The right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [in] X ``` X is DOUBLE PRECISION array, dimension (LDX,NRHS) The solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] FERR ``` FERR is DOUBLE PRECISION array, dimension (NRHS) The estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error.``` [out] BERR ``` BERR is DOUBLE PRECISION array, dimension (NRHS) The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).``` [out] WORK ` WORK is DOUBLE PRECISION array, dimension (3*N)` [out] IWORK ` IWORK is INTEGER array, dimension (N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```

Definition at line 180 of file dtrrfs.f.

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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dtrmv(uplo, trans, diag, n, a, lda, x, incx)
DTRMV
Definition dtrmv.f:147
subroutine dtrsv(uplo, trans, diag, n, a, lda, x, incx)
DTRSV
Definition dtrsv.f:143
Here is the call graph for this function:
Here is the caller graph for this function: