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

◆ ztrrfs()

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

ZTRRFS

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

Purpose:
 ZTRRFS 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 ZTRTRS or some other
 means before entering this routine.  ZTRRFS 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)
[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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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 ztrrfs.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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
193 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ),
194 $ X( LDX, * )
195* ..
196*
197* =====================================================================
198*
199* .. Parameters ..
200 DOUBLE PRECISION ZERO
201 parameter( zero = 0.0d+0 )
202 COMPLEX*16 ONE
203 parameter( one = ( 1.0d+0, 0.0d+0 ) )
204* ..
205* .. Local Scalars ..
206 LOGICAL NOTRAN, NOUNIT, UPPER
207 CHARACTER TRANSN, TRANST
208 INTEGER I, J, K, KASE, NZ
209 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
210 COMPLEX*16 ZDUM
211* ..
212* .. Local Arrays ..
213 INTEGER ISAVE( 3 )
214* ..
215* .. External Subroutines ..
216 EXTERNAL xerbla, zaxpy, zcopy, zlacn2, ztrmv, ztrsv
217* ..
218* .. Intrinsic Functions ..
219 INTRINSIC abs, dble, dimag, max
220* ..
221* .. External Functions ..
222 LOGICAL LSAME
223 DOUBLE PRECISION DLAMCH
224 EXTERNAL lsame, dlamch
225* ..
226* .. Statement Functions ..
227 DOUBLE PRECISION CABS1
228* ..
229* .. Statement Function definitions ..
230 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
231* ..
232* .. Executable Statements ..
233*
234* Test the input parameters.
235*
236 info = 0
237 upper = lsame( uplo, 'U' )
238 notran = lsame( trans, 'N' )
239 nounit = lsame( diag, 'N' )
240*
241 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
242 info = -1
243 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
244 $ lsame( trans, 'C' ) ) THEN
245 info = -2
246 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
247 info = -3
248 ELSE IF( n.LT.0 ) THEN
249 info = -4
250 ELSE IF( nrhs.LT.0 ) THEN
251 info = -5
252 ELSE IF( lda.LT.max( 1, n ) ) THEN
253 info = -7
254 ELSE IF( ldb.LT.max( 1, n ) ) THEN
255 info = -9
256 ELSE IF( ldx.LT.max( 1, n ) ) THEN
257 info = -11
258 END IF
259 IF( info.NE.0 ) THEN
260 CALL xerbla( 'ZTRRFS', -info )
261 RETURN
262 END IF
263*
264* Quick return if possible
265*
266 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
267 DO 10 j = 1, nrhs
268 ferr( j ) = zero
269 berr( j ) = zero
270 10 CONTINUE
271 RETURN
272 END IF
273*
274 IF( notran ) THEN
275 transn = 'N'
276 transt = 'C'
277 ELSE
278 transn = 'C'
279 transt = 'N'
280 END IF
281*
282* NZ = maximum number of nonzero elements in each row of A, plus 1
283*
284 nz = n + 1
285 eps = dlamch( 'Epsilon' )
286 safmin = dlamch( 'Safe minimum' )
287 safe1 = nz*safmin
288 safe2 = safe1 / eps
289*
290* Do for each right hand side
291*
292 DO 250 j = 1, nrhs
293*
294* Compute residual R = B - op(A) * X,
295* where op(A) = A, A**T, or A**H, depending on TRANS.
296*
297 CALL zcopy( n, x( 1, j ), 1, work, 1 )
298 CALL ztrmv( uplo, trans, diag, n, a, lda, work, 1 )
299 CALL zaxpy( n, -one, b( 1, j ), 1, work, 1 )
300*
301* Compute componentwise relative backward error from formula
302*
303* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
304*
305* where abs(Z) is the componentwise absolute value of the matrix
306* or vector Z. If the i-th component of the denominator is less
307* than SAFE2, then SAFE1 is added to the i-th components of the
308* numerator and denominator before dividing.
309*
310 DO 20 i = 1, n
311 rwork( i ) = cabs1( b( i, j ) )
312 20 CONTINUE
313*
314 IF( notran ) THEN
315*
316* Compute abs(A)*abs(X) + abs(B).
317*
318 IF( upper ) THEN
319 IF( nounit ) THEN
320 DO 40 k = 1, n
321 xk = cabs1( x( k, j ) )
322 DO 30 i = 1, k
323 rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
324 30 CONTINUE
325 40 CONTINUE
326 ELSE
327 DO 60 k = 1, n
328 xk = cabs1( x( k, j ) )
329 DO 50 i = 1, k - 1
330 rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
331 50 CONTINUE
332 rwork( k ) = rwork( k ) + xk
333 60 CONTINUE
334 END IF
335 ELSE
336 IF( nounit ) THEN
337 DO 80 k = 1, n
338 xk = cabs1( x( k, j ) )
339 DO 70 i = k, n
340 rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
341 70 CONTINUE
342 80 CONTINUE
343 ELSE
344 DO 100 k = 1, n
345 xk = cabs1( x( k, j ) )
346 DO 90 i = k + 1, n
347 rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
348 90 CONTINUE
349 rwork( k ) = rwork( k ) + xk
350 100 CONTINUE
351 END IF
352 END IF
353 ELSE
354*
355* Compute abs(A**H)*abs(X) + abs(B).
356*
357 IF( upper ) THEN
358 IF( nounit ) THEN
359 DO 120 k = 1, n
360 s = zero
361 DO 110 i = 1, k
362 s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
363 110 CONTINUE
364 rwork( k ) = rwork( k ) + s
365 120 CONTINUE
366 ELSE
367 DO 140 k = 1, n
368 s = cabs1( x( k, j ) )
369 DO 130 i = 1, k - 1
370 s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
371 130 CONTINUE
372 rwork( k ) = rwork( k ) + s
373 140 CONTINUE
374 END IF
375 ELSE
376 IF( nounit ) THEN
377 DO 160 k = 1, n
378 s = zero
379 DO 150 i = k, n
380 s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
381 150 CONTINUE
382 rwork( k ) = rwork( k ) + s
383 160 CONTINUE
384 ELSE
385 DO 180 k = 1, n
386 s = cabs1( x( k, j ) )
387 DO 170 i = k + 1, n
388 s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
389 170 CONTINUE
390 rwork( k ) = rwork( k ) + s
391 180 CONTINUE
392 END IF
393 END IF
394 END IF
395 s = zero
396 DO 190 i = 1, n
397 IF( rwork( i ).GT.safe2 ) THEN
398 s = max( s, cabs1( work( i ) ) / rwork( i ) )
399 ELSE
400 s = max( s, ( cabs1( work( i ) )+safe1 ) /
401 $ ( rwork( i )+safe1 ) )
402 END IF
403 190 CONTINUE
404 berr( j ) = s
405*
406* Bound error from formula
407*
408* norm(X - XTRUE) / norm(X) .le. FERR =
409* norm( abs(inv(op(A)))*
410* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
411*
412* where
413* norm(Z) is the magnitude of the largest component of Z
414* inv(op(A)) is the inverse of op(A)
415* abs(Z) is the componentwise absolute value of the matrix or
416* vector Z
417* NZ is the maximum number of nonzeros in any row of A, plus 1
418* EPS is machine epsilon
419*
420* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
421* is incremented by SAFE1 if the i-th component of
422* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
423*
424* Use ZLACN2 to estimate the infinity-norm of the matrix
425* inv(op(A)) * diag(W),
426* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
427*
428 DO 200 i = 1, n
429 IF( rwork( i ).GT.safe2 ) THEN
430 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
431 ELSE
432 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
433 $ safe1
434 END IF
435 200 CONTINUE
436*
437 kase = 0
438 210 CONTINUE
439 CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
440 IF( kase.NE.0 ) THEN
441 IF( kase.EQ.1 ) THEN
442*
443* Multiply by diag(W)*inv(op(A)**H).
444*
445 CALL ztrsv( uplo, transt, diag, n, a, lda, work, 1 )
446 DO 220 i = 1, n
447 work( i ) = rwork( i )*work( i )
448 220 CONTINUE
449 ELSE
450*
451* Multiply by inv(op(A))*diag(W).
452*
453 DO 230 i = 1, n
454 work( i ) = rwork( i )*work( i )
455 230 CONTINUE
456 CALL ztrsv( uplo, transn, diag, n, a, lda, work, 1 )
457 END IF
458 GO TO 210
459 END IF
460*
461* Normalize error.
462*
463 lstres = zero
464 DO 240 i = 1, n
465 lstres = max( lstres, cabs1( x( i, j ) ) )
466 240 CONTINUE
467 IF( lstres.NE.zero )
468 $ ferr( j ) = ferr( j ) / lstres
469*
470 250 CONTINUE
471*
472 RETURN
473*
474* End of ZTRRFS
475*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zlacn2
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:133
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztrmv
subroutine ztrmv(uplo, trans, diag, n, a, lda, x, incx)
ZTRMV
Definition ztrmv.f:147
ztrsv
subroutine ztrsv(uplo, trans, diag, n, a, lda, x, incx)
ZTRSV
Definition ztrsv.f:149
Here is the call graph for this function:
Here is the caller graph for this function: