LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctbrfs ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  KD,
integer  NRHS,
complex, dimension( ldab, * )  AB,
integer  LDAB,
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 
)

CTBRFS

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

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

 The solution matrix X must be computed by CTBTRS or some other
 means before entering this routine.  CTBRFS 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]KD
          KD is INTEGER
          The number of superdiagonals or subdiagonals of the
          triangular band matrix A.  KD >= 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]AB
          AB is COMPLEX array, dimension (LDAB,N)
          The upper or lower triangular band matrix A, stored in the
          first kd+1 rows of the array. The j-th column of A is stored
          in the j-th column of the array AB as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.
[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 190 of file ctbrfs.f.

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