LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
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

Download DTRRFS + dependencies [TGZ] [ZIP] [TXT]

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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
Here is the call graph for this function:
Here is the caller graph for this function: