LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ctrrfs()

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.

Definition at line 180 of file ctrrfs.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  REAL BERR( * ), FERR( * ), RWORK( * )
193  COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ),
194  $ X( LDX, * )
195 * ..
196 *
197 * =====================================================================
198 *
199 * .. Parameters ..
200  REAL ZERO
201  parameter( zero = 0.0e+0 )
202  COMPLEX ONE
203  parameter( one = ( 1.0e+0, 0.0e+0 ) )
204 * ..
205 * .. Local Scalars ..
206  LOGICAL NOTRAN, NOUNIT, UPPER
207  CHARACTER TRANSN, TRANST
208  INTEGER I, J, K, KASE, NZ
209  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
210  COMPLEX ZDUM
211 * ..
212 * .. Local Arrays ..
213  INTEGER ISAVE( 3 )
214 * ..
215 * .. External Subroutines ..
216  EXTERNAL caxpy, ccopy, clacn2, ctrmv, ctrsv, xerbla
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC abs, aimag, max, real
220 * ..
221 * .. External Functions ..
222  LOGICAL LSAME
223  REAL SLAMCH
224  EXTERNAL lsame, slamch
225 * ..
226 * .. Statement Functions ..
227  REAL CABS1
228 * ..
229 * .. Statement Function definitions ..
230  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CTRRFS', -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 = slamch( 'Epsilon' )
286  safmin = slamch( '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 ccopy( n, x( 1, j ), 1, work, 1 )
298  CALL ctrmv( uplo, trans, diag, n, a, lda, work, 1 )
299  CALL caxpy( 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 CLACN2 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 clacn2( 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 ctrsv( 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 ctrsv( 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 CTRRFS
475 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine ctrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRSV
Definition: ctrsv.f:149
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:147
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:133
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: