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

◆ 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)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
Here is the call graph for this function:
Here is the caller graph for this function: