LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ctbrfs()

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.

Definition at line 186 of file ctbrfs.f.

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