LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctbrfs.f
Go to the documentation of this file.
1*> \brief \b CTBRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTBRFS + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctbrfs.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctbrfs.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctbrfs.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
22* LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER DIAG, TRANS, UPLO
26* INTEGER INFO, KD, LDAB, LDB, LDX, N, NRHS
27* ..
28* .. Array Arguments ..
29* REAL BERR( * ), FERR( * ), RWORK( * )
30* COMPLEX AB( LDAB, * ), B( LDB, * ), WORK( * ),
31* $ X( LDX, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> CTBRFS provides error bounds and backward error estimates for the
41*> solution to a system of linear equations with a triangular band
42*> coefficient matrix.
43*>
44*> The solution matrix X must be computed by CTBTRS or some other
45*> means before entering this routine. CTBRFS does not do iterative
46*> refinement because doing so cannot improve the backward error.
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] UPLO
53*> \verbatim
54*> UPLO is CHARACTER*1
55*> = 'U': A is upper triangular;
56*> = 'L': A is lower triangular.
57*> \endverbatim
58*>
59*> \param[in] TRANS
60*> \verbatim
61*> TRANS is CHARACTER*1
62*> Specifies the form of the system of equations:
63*> = 'N': A * X = B (No transpose)
64*> = 'T': A**T * X = B (Transpose)
65*> = 'C': A**H * X = B (Conjugate transpose)
66*> \endverbatim
67*>
68*> \param[in] DIAG
69*> \verbatim
70*> DIAG is CHARACTER*1
71*> = 'N': A is non-unit triangular;
72*> = 'U': A is unit triangular.
73*> \endverbatim
74*>
75*> \param[in] N
76*> \verbatim
77*> N is INTEGER
78*> The order of the matrix A. N >= 0.
79*> \endverbatim
80*>
81*> \param[in] KD
82*> \verbatim
83*> KD is INTEGER
84*> The number of superdiagonals or subdiagonals of the
85*> triangular band matrix A. KD >= 0.
86*> \endverbatim
87*>
88*> \param[in] NRHS
89*> \verbatim
90*> NRHS is INTEGER
91*> The number of right hand sides, i.e., the number of columns
92*> of the matrices B and X. NRHS >= 0.
93*> \endverbatim
94*>
95*> \param[in] AB
96*> \verbatim
97*> AB is COMPLEX array, dimension (LDAB,N)
98*> The upper or lower triangular band matrix A, stored in the
99*> first kd+1 rows of the array. The j-th column of A is stored
100*> in the j-th column of the array AB as follows:
101*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
102*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
103*> If DIAG = 'U', the diagonal elements of A are not referenced
104*> and are assumed to be 1.
105*> \endverbatim
106*>
107*> \param[in] LDAB
108*> \verbatim
109*> LDAB is INTEGER
110*> The leading dimension of the array AB. LDAB >= KD+1.
111*> \endverbatim
112*>
113*> \param[in] B
114*> \verbatim
115*> B is COMPLEX array, dimension (LDB,NRHS)
116*> The right hand side matrix B.
117*> \endverbatim
118*>
119*> \param[in] LDB
120*> \verbatim
121*> LDB is INTEGER
122*> The leading dimension of the array B. LDB >= max(1,N).
123*> \endverbatim
124*>
125*> \param[in] X
126*> \verbatim
127*> X is COMPLEX array, dimension (LDX,NRHS)
128*> The solution matrix X.
129*> \endverbatim
130*>
131*> \param[in] LDX
132*> \verbatim
133*> LDX is INTEGER
134*> The leading dimension of the array X. LDX >= max(1,N).
135*> \endverbatim
136*>
137*> \param[out] FERR
138*> \verbatim
139*> FERR is REAL array, dimension (NRHS)
140*> The estimated forward error bound for each solution vector
141*> X(j) (the j-th column of the solution matrix X).
142*> If XTRUE is the true solution corresponding to X(j), FERR(j)
143*> is an estimated upper bound for the magnitude of the largest
144*> element in (X(j) - XTRUE) divided by the magnitude of the
145*> largest element in X(j). The estimate is as reliable as
146*> the estimate for RCOND, and is almost always a slight
147*> overestimate of the true error.
148*> \endverbatim
149*>
150*> \param[out] BERR
151*> \verbatim
152*> BERR is REAL array, dimension (NRHS)
153*> The componentwise relative backward error of each solution
154*> vector X(j) (i.e., the smallest relative change in
155*> any element of A or B that makes X(j) an exact solution).
156*> \endverbatim
157*>
158*> \param[out] WORK
159*> \verbatim
160*> WORK is COMPLEX array, dimension (2*N)
161*> \endverbatim
162*>
163*> \param[out] RWORK
164*> \verbatim
165*> RWORK is REAL array, dimension (N)
166*> \endverbatim
167*>
168*> \param[out] INFO
169*> \verbatim
170*> INFO is INTEGER
171*> = 0: successful exit
172*> < 0: if INFO = -i, the i-th argument had an illegal value
173*> \endverbatim
174*
175* Authors:
176* ========
177*
178*> \author Univ. of Tennessee
179*> \author Univ. of California Berkeley
180*> \author Univ. of Colorado Denver
181*> \author NAG Ltd.
182*
183*> \ingroup tbrfs
184*
185* =====================================================================
186 SUBROUTINE ctbrfs( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
187 $ LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
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*
494 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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:133
ctbmv
subroutine ctbmv(uplo, trans, diag, n, k, a, lda, x, incx)
CTBMV
Definition ctbmv.f:186
ctbrfs
subroutine ctbrfs(uplo, trans, diag, n, kd, nrhs, ab, ldab, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CTBRFS
Definition ctbrfs.f:188
ctbsv
subroutine ctbsv(uplo, trans, diag, n, k, a, lda, x, incx)
CTBSV
Definition ctbsv.f:189