LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctrrfs ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CTRRFS

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

Purpose: 
 CTRRFS 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 CTRTRS or some other
 means before entering this routine.  CTRRFS 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 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 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 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 REAL 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 REAL 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 array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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.
Date
November 2011

Definition at line 184 of file ctrrfs.f.

184 *
185 * -- LAPACK computational routine (version 3.4.0) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * November 2011
189 *
190 * .. Scalar Arguments ..
191  CHARACTER diag, trans, uplo
192  INTEGER info, lda, ldb, ldx, n, nrhs
193 * ..
194 * .. Array Arguments ..
195  REAL berr( * ), ferr( * ), rwork( * )
196  COMPLEX a( lda, * ), b( ldb, * ), work( * ),
197  $ x( ldx, * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  REAL zero
204  parameter ( zero = 0.0e+0 )
205  COMPLEX one
206  parameter ( one = ( 1.0e+0, 0.0e+0 ) )
207 * ..
208 * .. Local Scalars ..
209  LOGICAL notran, nounit, upper
210  CHARACTER transn, transt
211  INTEGER i, j, k, kase, nz
212  REAL eps, lstres, s, safe1, safe2, safmin, xk
213  COMPLEX zdum
214 * ..
215 * .. Local Arrays ..
216  INTEGER isave( 3 )
217 * ..
218 * .. External Subroutines ..
219  EXTERNAL caxpy, ccopy, clacn2, ctrmv, ctrsv, xerbla
220 * ..
221 * .. Intrinsic Functions ..
222  INTRINSIC abs, aimag, max, real
223 * ..
224 * .. External Functions ..
225  LOGICAL lsame
226  REAL slamch
227  EXTERNAL lsame, slamch
228 * ..
229 * .. Statement Functions ..
230  REAL cabs1
231 * ..
232 * .. Statement Function definitions ..
233  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
234 * ..
235 * .. Executable Statements ..
236 *
237 * Test the input parameters.
238 *
239  info = 0
240  upper = lsame( uplo, 'U' )
241  notran = lsame( trans, 'N' )
242  nounit = lsame( diag, 'N' )
243 *
244  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
245  info = -1
246  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
247  $ lsame( trans, 'C' ) ) THEN
248  info = -2
249  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
250  info = -3
251  ELSE IF( n.LT.0 ) THEN
252  info = -4
253  ELSE IF( nrhs.LT.0 ) THEN
254  info = -5
255  ELSE IF( lda.LT.max( 1, n ) ) THEN
256  info = -7
257  ELSE IF( ldb.LT.max( 1, n ) ) THEN
258  info = -9
259  ELSE IF( ldx.LT.max( 1, n ) ) THEN
260  info = -11
261  END IF
262  IF( info.NE.0 ) THEN
263  CALL xerbla( 'CTRRFS', -info )
264  RETURN
265  END IF
266 *
267 * Quick return if possible
268 *
269  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
270  DO 10 j = 1, nrhs
271  ferr( j ) = zero
272  berr( j ) = zero
273  10 CONTINUE
274  RETURN
275  END IF
276 *
277  IF( notran ) THEN
278  transn = 'N'
279  transt = 'C'
280  ELSE
281  transn = 'C'
282  transt = 'N'
283  END IF
284 *
285 * NZ = maximum number of nonzero elements in each row of A, plus 1
286 *
287  nz = n + 1
288  eps = slamch( 'Epsilon' )
289  safmin = slamch( 'Safe minimum' )
290  safe1 = nz*safmin
291  safe2 = safe1 / eps
292 *
293 * Do for each right hand side
294 *
295  DO 250 j = 1, nrhs
296 *
297 * Compute residual R = B - op(A) * X,
298 * where op(A) = A, A**T, or A**H, depending on TRANS.
299 *
300  CALL ccopy( n, x( 1, j ), 1, work, 1 )
301  CALL ctrmv( uplo, trans, diag, n, a, lda, work, 1 )
302  CALL caxpy( n, -one, b( 1, j ), 1, work, 1 )
303 *
304 * Compute componentwise relative backward error from formula
305 *
306 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
307 *
308 * where abs(Z) is the componentwise absolute value of the matrix
309 * or vector Z. If the i-th component of the denominator is less
310 * than SAFE2, then SAFE1 is added to the i-th components of the
311 * numerator and denominator before dividing.
312 *
313  DO 20 i = 1, n
314  rwork( i ) = cabs1( b( i, j ) )
315  20 CONTINUE
316 *
317  IF( notran ) THEN
318 *
319 * Compute abs(A)*abs(X) + abs(B).
320 *
321  IF( upper ) THEN
322  IF( nounit ) THEN
323  DO 40 k = 1, n
324  xk = cabs1( x( k, j ) )
325  DO 30 i = 1, k
326  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
327  30 CONTINUE
328  40 CONTINUE
329  ELSE
330  DO 60 k = 1, n
331  xk = cabs1( x( k, j ) )
332  DO 50 i = 1, k - 1
333  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
334  50 CONTINUE
335  rwork( k ) = rwork( k ) + xk
336  60 CONTINUE
337  END IF
338  ELSE
339  IF( nounit ) THEN
340  DO 80 k = 1, n
341  xk = cabs1( x( k, j ) )
342  DO 70 i = k, n
343  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
344  70 CONTINUE
345  80 CONTINUE
346  ELSE
347  DO 100 k = 1, n
348  xk = cabs1( x( k, j ) )
349  DO 90 i = k + 1, n
350  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
351  90 CONTINUE
352  rwork( k ) = rwork( k ) + xk
353  100 CONTINUE
354  END IF
355  END IF
356  ELSE
357 *
358 * Compute abs(A**H)*abs(X) + abs(B).
359 *
360  IF( upper ) THEN
361  IF( nounit ) THEN
362  DO 120 k = 1, n
363  s = zero
364  DO 110 i = 1, k
365  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
366  110 CONTINUE
367  rwork( k ) = rwork( k ) + s
368  120 CONTINUE
369  ELSE
370  DO 140 k = 1, n
371  s = cabs1( x( k, j ) )
372  DO 130 i = 1, k - 1
373  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
374  130 CONTINUE
375  rwork( k ) = rwork( k ) + s
376  140 CONTINUE
377  END IF
378  ELSE
379  IF( nounit ) THEN
380  DO 160 k = 1, n
381  s = zero
382  DO 150 i = k, n
383  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
384  150 CONTINUE
385  rwork( k ) = rwork( k ) + s
386  160 CONTINUE
387  ELSE
388  DO 180 k = 1, n
389  s = cabs1( x( k, j ) )
390  DO 170 i = k + 1, n
391  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
392  170 CONTINUE
393  rwork( k ) = rwork( k ) + s
394  180 CONTINUE
395  END IF
396  END IF
397  END IF
398  s = zero
399  DO 190 i = 1, n
400  IF( rwork( i ).GT.safe2 ) THEN
401  s = max( s, cabs1( work( i ) ) / rwork( i ) )
402  ELSE
403  s = max( s, ( cabs1( work( i ) )+safe1 ) /
404  $ ( rwork( i )+safe1 ) )
405  END IF
406  190 CONTINUE
407  berr( j ) = s
408 *
409 * Bound error from formula
410 *
411 * norm(X - XTRUE) / norm(X) .le. FERR =
412 * norm( abs(inv(op(A)))*
413 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
414 *
415 * where
416 * norm(Z) is the magnitude of the largest component of Z
417 * inv(op(A)) is the inverse of op(A)
418 * abs(Z) is the componentwise absolute value of the matrix or
419 * vector Z
420 * NZ is the maximum number of nonzeros in any row of A, plus 1
421 * EPS is machine epsilon
422 *
423 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
424 * is incremented by SAFE1 if the i-th component of
425 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
426 *
427 * Use CLACN2 to estimate the infinity-norm of the matrix
428 * inv(op(A)) * diag(W),
429 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
430 *
431  DO 200 i = 1, n
432  IF( rwork( i ).GT.safe2 ) THEN
433  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
434  ELSE
435  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
436  $ safe1
437  END IF
438  200 CONTINUE
439 *
440  kase = 0
441  210 CONTINUE
442  CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
443  IF( kase.NE.0 ) THEN
444  IF( kase.EQ.1 ) THEN
445 *
446 * Multiply by diag(W)*inv(op(A)**H).
447 *
448  CALL ctrsv( uplo, transt, diag, n, a, lda, work, 1 )
449  DO 220 i = 1, n
450  work( i ) = rwork( i )*work( i )
451  220 CONTINUE
452  ELSE
453 *
454 * Multiply by inv(op(A))*diag(W).
455 *
456  DO 230 i = 1, n
457  work( i ) = rwork( i )*work( i )
458  230 CONTINUE
459  CALL ctrsv( uplo, transn, diag, n, a, lda, work, 1 )
460  END IF
461  GO TO 210
462  END IF
463 *
464 * Normalize error.
465 *
466  lstres = zero
467  DO 240 i = 1, n
468  lstres = max( lstres, cabs1( x( i, j ) ) )
469  240 CONTINUE
470  IF( lstres.NE.zero )
471  $ ferr( j ) = ferr( j ) / lstres
472 *
473  250 CONTINUE
474 *
475  RETURN
476 *
477 * End of CTRRFS
478 *
ctrsv
subroutine ctrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRSV
Definition: ctrsv.f:151
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ctrmv
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:149
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
clacn2
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:135

Here is the call graph for this function:

Here is the caller graph for this function: