LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztbrfs.f
Go to the documentation of this file.
1*> \brief \b ZTBRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZTBRFS + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztbrfs.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztbrfs.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztbrfs.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
20* LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER DIAG, TRANS, UPLO
24* INTEGER INFO, KD, LDAB, LDB, LDX, N, NRHS
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
28* COMPLEX*16 AB( LDAB, * ), B( LDB, * ), WORK( * ),
29* $ X( LDX, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZTBRFS provides error bounds and backward error estimates for the
39*> solution to a system of linear equations with a triangular band
40*> coefficient matrix.
41*>
42*> The solution matrix X must be computed by ZTBTRS or some other
43*> means before entering this routine. ZTBRFS does not do iterative
44*> refinement because doing so cannot improve the backward error.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] UPLO
51*> \verbatim
52*> UPLO is CHARACTER*1
53*> = 'U': A is upper triangular;
54*> = 'L': A is lower triangular.
55*> \endverbatim
56*>
57*> \param[in] TRANS
58*> \verbatim
59*> TRANS is CHARACTER*1
60*> Specifies the form of the system of equations:
61*> = 'N': A * X = B (No transpose)
62*> = 'T': A**T * X = B (Transpose)
63*> = 'C': A**H * X = B (Conjugate transpose)
64*> \endverbatim
65*>
66*> \param[in] DIAG
67*> \verbatim
68*> DIAG is CHARACTER*1
69*> = 'N': A is non-unit triangular;
70*> = 'U': A is unit triangular.
71*> \endverbatim
72*>
73*> \param[in] N
74*> \verbatim
75*> N is INTEGER
76*> The order of the matrix A. N >= 0.
77*> \endverbatim
78*>
79*> \param[in] KD
80*> \verbatim
81*> KD is INTEGER
82*> The number of superdiagonals or subdiagonals of the
83*> triangular band matrix A. KD >= 0.
84*> \endverbatim
85*>
86*> \param[in] NRHS
87*> \verbatim
88*> NRHS is INTEGER
89*> The number of right hand sides, i.e., the number of columns
90*> of the matrices B and X. NRHS >= 0.
91*> \endverbatim
92*>
93*> \param[in] AB
94*> \verbatim
95*> AB is COMPLEX*16 array, dimension (LDAB,N)
96*> The upper or lower triangular band matrix A, stored in the
97*> first kd+1 rows of the array. The j-th column of A is stored
98*> in the j-th column of the array AB as follows:
99*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
100*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
101*> If DIAG = 'U', the diagonal elements of A are not referenced
102*> and are assumed to be 1.
103*> \endverbatim
104*>
105*> \param[in] LDAB
106*> \verbatim
107*> LDAB is INTEGER
108*> The leading dimension of the array AB. LDAB >= KD+1.
109*> \endverbatim
110*>
111*> \param[in] B
112*> \verbatim
113*> B is COMPLEX*16 array, dimension (LDB,NRHS)
114*> The right hand side matrix B.
115*> \endverbatim
116*>
117*> \param[in] LDB
118*> \verbatim
119*> LDB is INTEGER
120*> The leading dimension of the array B. LDB >= max(1,N).
121*> \endverbatim
122*>
123*> \param[in] X
124*> \verbatim
125*> X is COMPLEX*16 array, dimension (LDX,NRHS)
126*> The solution matrix X.
127*> \endverbatim
128*>
129*> \param[in] LDX
130*> \verbatim
131*> LDX is INTEGER
132*> The leading dimension of the array X. LDX >= max(1,N).
133*> \endverbatim
134*>
135*> \param[out] FERR
136*> \verbatim
137*> FERR is DOUBLE PRECISION array, dimension (NRHS)
138*> The estimated forward error bound for each solution vector
139*> X(j) (the j-th column of the solution matrix X).
140*> If XTRUE is the true solution corresponding to X(j), FERR(j)
141*> is an estimated upper bound for the magnitude of the largest
142*> element in (X(j) - XTRUE) divided by the magnitude of the
143*> largest element in X(j). The estimate is as reliable as
144*> the estimate for RCOND, and is almost always a slight
145*> overestimate of the true error.
146*> \endverbatim
147*>
148*> \param[out] BERR
149*> \verbatim
150*> BERR is DOUBLE PRECISION array, dimension (NRHS)
151*> The componentwise relative backward error of each solution
152*> vector X(j) (i.e., the smallest relative change in
153*> any element of A or B that makes X(j) an exact solution).
154*> \endverbatim
155*>
156*> \param[out] WORK
157*> \verbatim
158*> WORK is COMPLEX*16 array, dimension (2*N)
159*> \endverbatim
160*>
161*> \param[out] RWORK
162*> \verbatim
163*> RWORK is DOUBLE PRECISION array, dimension (N)
164*> \endverbatim
165*>
166*> \param[out] INFO
167*> \verbatim
168*> INFO is INTEGER
169*> = 0: successful exit
170*> < 0: if INFO = -i, the i-th argument had an illegal value
171*> \endverbatim
172*
173* Authors:
174* ========
175*
176*> \author Univ. of Tennessee
177*> \author Univ. of California Berkeley
178*> \author Univ. of Colorado Denver
179*> \author NAG Ltd.
180*
181*> \ingroup tbrfs
182*
183* =====================================================================
184 SUBROUTINE ztbrfs( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
185 $ LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
186*
187* -- LAPACK computational routine --
188* -- LAPACK is a software package provided by Univ. of Tennessee, --
189* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190*
191* .. Scalar Arguments ..
192 CHARACTER DIAG, TRANS, UPLO
193 INTEGER INFO, KD, LDAB, LDB, LDX, N, NRHS
194* ..
195* .. Array Arguments ..
196 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
197 COMPLEX*16 AB( LDAB, * ), B( LDB, * ), WORK( * ),
198 $ x( ldx, * )
199* ..
200*
201* =====================================================================
202*
203* .. Parameters ..
204 DOUBLE PRECISION ZERO
205 parameter( zero = 0.0d+0 )
206 COMPLEX*16 ONE
207 parameter( one = ( 1.0d+0, 0.0d+0 ) )
208* ..
209* .. Local Scalars ..
210 LOGICAL NOTRAN, NOUNIT, UPPER
211 CHARACTER TRANSN, TRANST
212 INTEGER I, J, K, KASE, NZ
213 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
214 COMPLEX*16 ZDUM
215* ..
216* .. Local Arrays ..
217 INTEGER ISAVE( 3 )
218* ..
219* .. External Subroutines ..
220 EXTERNAL xerbla, zaxpy, zcopy, zlacn2, ztbmv,
221 $ ztbsv
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC abs, dble, dimag, max, min
225* ..
226* .. External Functions ..
227 LOGICAL LSAME
228 DOUBLE PRECISION DLAMCH
229 EXTERNAL lsame, dlamch
230* ..
231* .. Statement Functions ..
232 DOUBLE PRECISION CABS1
233* ..
234* .. Statement Function definitions ..
235 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
236* ..
237* .. Executable Statements ..
238*
239* Test the input parameters.
240*
241 info = 0
242 upper = lsame( uplo, 'U' )
243 notran = lsame( trans, 'N' )
244 nounit = lsame( diag, 'N' )
245*
246 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
247 info = -1
248 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
249 $ lsame( trans, 'C' ) ) THEN
250 info = -2
251 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
252 info = -3
253 ELSE IF( n.LT.0 ) THEN
254 info = -4
255 ELSE IF( kd.LT.0 ) THEN
256 info = -5
257 ELSE IF( nrhs.LT.0 ) THEN
258 info = -6
259 ELSE IF( ldab.LT.kd+1 ) THEN
260 info = -8
261 ELSE IF( ldb.LT.max( 1, n ) ) THEN
262 info = -10
263 ELSE IF( ldx.LT.max( 1, n ) ) THEN
264 info = -12
265 END IF
266 IF( info.NE.0 ) THEN
267 CALL xerbla( 'ZTBRFS', -info )
268 RETURN
269 END IF
270*
271* Quick return if possible
272*
273 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
274 DO 10 j = 1, nrhs
275 ferr( j ) = zero
276 berr( j ) = zero
277 10 CONTINUE
278 RETURN
279 END IF
280*
281 IF( notran ) THEN
282 transn = 'N'
283 transt = 'C'
284 ELSE
285 transn = 'C'
286 transt = 'N'
287 END IF
288*
289* NZ = maximum number of nonzero elements in each row of A, plus 1
290*
291 nz = kd + 2
292 eps = dlamch( 'Epsilon' )
293 safmin = dlamch( 'Safe minimum' )
294 safe1 = nz*safmin
295 safe2 = safe1 / eps
296*
297* Do for each right hand side
298*
299 DO 250 j = 1, nrhs
300*
301* Compute residual R = B - op(A) * X,
302* where op(A) = A, A**T, or A**H, depending on TRANS.
303*
304 CALL zcopy( n, x( 1, j ), 1, work, 1 )
305 CALL ztbmv( uplo, trans, diag, n, kd, ab, ldab, work, 1 )
306 CALL zaxpy( n, -one, b( 1, j ), 1, work, 1 )
307*
308* Compute componentwise relative backward error from formula
309*
310* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
311*
312* where abs(Z) is the componentwise absolute value of the matrix
313* or vector Z. If the i-th component of the denominator is less
314* than SAFE2, then SAFE1 is added to the i-th components of the
315* numerator and denominator before dividing.
316*
317 DO 20 i = 1, n
318 rwork( i ) = cabs1( b( i, j ) )
319 20 CONTINUE
320*
321 IF( notran ) THEN
322*
323* Compute abs(A)*abs(X) + abs(B).
324*
325 IF( upper ) THEN
326 IF( nounit ) THEN
327 DO 40 k = 1, n
328 xk = cabs1( x( k, j ) )
329 DO 30 i = max( 1, k-kd ), k
330 rwork( i ) = rwork( i ) +
331 $ cabs1( ab( kd+1+i-k, k ) )*xk
332 30 CONTINUE
333 40 CONTINUE
334 ELSE
335 DO 60 k = 1, n
336 xk = cabs1( x( k, j ) )
337 DO 50 i = max( 1, k-kd ), k - 1
338 rwork( i ) = rwork( i ) +
339 $ cabs1( ab( kd+1+i-k, k ) )*xk
340 50 CONTINUE
341 rwork( k ) = rwork( k ) + xk
342 60 CONTINUE
343 END IF
344 ELSE
345 IF( nounit ) THEN
346 DO 80 k = 1, n
347 xk = cabs1( x( k, j ) )
348 DO 70 i = k, min( n, k+kd )
349 rwork( i ) = rwork( i ) +
350 $ cabs1( ab( 1+i-k, k ) )*xk
351 70 CONTINUE
352 80 CONTINUE
353 ELSE
354 DO 100 k = 1, n
355 xk = cabs1( x( k, j ) )
356 DO 90 i = k + 1, min( n, k+kd )
357 rwork( i ) = rwork( i ) +
358 $ cabs1( ab( 1+i-k, k ) )*xk
359 90 CONTINUE
360 rwork( k ) = rwork( k ) + xk
361 100 CONTINUE
362 END IF
363 END IF
364 ELSE
365*
366* Compute abs(A**H)*abs(X) + abs(B).
367*
368 IF( upper ) THEN
369 IF( nounit ) THEN
370 DO 120 k = 1, n
371 s = zero
372 DO 110 i = max( 1, k-kd ), k
373 s = s + cabs1( ab( kd+1+i-k, k ) )*
374 $ cabs1( x( i, j ) )
375 110 CONTINUE
376 rwork( k ) = rwork( k ) + s
377 120 CONTINUE
378 ELSE
379 DO 140 k = 1, n
380 s = cabs1( x( k, j ) )
381 DO 130 i = max( 1, k-kd ), k - 1
382 s = s + cabs1( ab( kd+1+i-k, k ) )*
383 $ cabs1( x( i, j ) )
384 130 CONTINUE
385 rwork( k ) = rwork( k ) + s
386 140 CONTINUE
387 END IF
388 ELSE
389 IF( nounit ) THEN
390 DO 160 k = 1, n
391 s = zero
392 DO 150 i = k, min( n, k+kd )
393 s = s + cabs1( ab( 1+i-k, k ) )*
394 $ cabs1( x( i, j ) )
395 150 CONTINUE
396 rwork( k ) = rwork( k ) + s
397 160 CONTINUE
398 ELSE
399 DO 180 k = 1, n
400 s = cabs1( x( k, j ) )
401 DO 170 i = k + 1, min( n, k+kd )
402 s = s + cabs1( ab( 1+i-k, k ) )*
403 $ cabs1( x( i, j ) )
404 170 CONTINUE
405 rwork( k ) = rwork( k ) + s
406 180 CONTINUE
407 END IF
408 END IF
409 END IF
410 s = zero
411 DO 190 i = 1, n
412 IF( rwork( i ).GT.safe2 ) THEN
413 s = max( s, cabs1( work( i ) ) / rwork( i ) )
414 ELSE
415 s = max( s, ( cabs1( work( i ) )+safe1 ) /
416 $ ( rwork( i )+safe1 ) )
417 END IF
418 190 CONTINUE
419 berr( j ) = s
420*
421* Bound error from formula
422*
423* norm(X - XTRUE) / norm(X) .le. FERR =
424* norm( abs(inv(op(A)))*
425* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
426*
427* where
428* norm(Z) is the magnitude of the largest component of Z
429* inv(op(A)) is the inverse of op(A)
430* abs(Z) is the componentwise absolute value of the matrix or
431* vector Z
432* NZ is the maximum number of nonzeros in any row of A, plus 1
433* EPS is machine epsilon
434*
435* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
436* is incremented by SAFE1 if the i-th component of
437* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
438*
439* Use ZLACN2 to estimate the infinity-norm of the matrix
440* inv(op(A)) * diag(W),
441* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
442*
443 DO 200 i = 1, n
444 IF( rwork( i ).GT.safe2 ) THEN
445 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
446 ELSE
447 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
448 $ safe1
449 END IF
450 200 CONTINUE
451*
452 kase = 0
453 210 CONTINUE
454 CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
455 IF( kase.NE.0 ) THEN
456 IF( kase.EQ.1 ) THEN
457*
458* Multiply by diag(W)*inv(op(A)**H).
459*
460 CALL ztbsv( uplo, transt, diag, n, kd, ab, ldab, work,
461 $ 1 )
462 DO 220 i = 1, n
463 work( i ) = rwork( i )*work( i )
464 220 CONTINUE
465 ELSE
466*
467* Multiply by inv(op(A))*diag(W).
468*
469 DO 230 i = 1, n
470 work( i ) = rwork( i )*work( i )
471 230 CONTINUE
472 CALL ztbsv( uplo, transn, diag, n, kd, ab, ldab, work,
473 $ 1 )
474 END IF
475 GO TO 210
476 END IF
477*
478* Normalize error.
479*
480 lstres = zero
481 DO 240 i = 1, n
482 lstres = max( lstres, cabs1( x( i, j ) ) )
483 240 CONTINUE
484 IF( lstres.NE.zero )
485 $ ferr( j ) = ferr( j ) / lstres
486*
487 250 CONTINUE
488*
489 RETURN
490*
491* End of ZTBRFS
492*
493 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zlacn2
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:131
ztbmv
subroutine ztbmv(uplo, trans, diag, n, k, a, lda, x, incx)
ZTBMV
Definition ztbmv.f:186
ztbrfs
subroutine ztbrfs(uplo, trans, diag, n, kd, nrhs, ab, ldab, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZTBRFS
Definition ztbrfs.f:186
ztbsv
subroutine ztbsv(uplo, trans, diag, n, k, a, lda, x, incx)
ZTBSV
Definition ztbsv.f:189