LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \htmlonly
9 *> Download ZTBRFS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztbrfs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztbrfs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztbrfs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZTBRFS( 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 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
30 * COMPLEX*16 AB( LDAB, * ), B( LDB, * ), WORK( * ),
31 * $ X( LDX, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZTBRFS 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 ZTBTRS or some other
45 *> means before entering this routine. ZTBRFS 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (2*N)
161 *> \endverbatim
162 *>
163 *> \param[out] RWORK
164 *> \verbatim
165 *> RWORK is DOUBLE PRECISION 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 *> \date November 2011
184 *
185 *> \ingroup complex16OTHERcomputational
186 *
187 * =====================================================================
188  SUBROUTINE ztbrfs( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
189  $ ldb, x, ldx, ferr, berr, work, rwork, info )
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  DOUBLE PRECISION berr( * ), ferr( * ), rwork( * )
202  COMPLEX*16 ab( ldab, * ), b( ldb, * ), work( * ),
203  $ x( ldx, * )
204 * ..
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209  DOUBLE PRECISION zero
210  parameter( zero = 0.0d+0 )
211  COMPLEX*16 one
212  parameter( one = ( 1.0d+0, 0.0d+0 ) )
213 * ..
214 * .. Local Scalars ..
215  LOGICAL notran, nounit, upper
216  CHARACTER transn, transt
217  INTEGER i, j, k, kase, nz
218  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
219  COMPLEX*16 zdum
220 * ..
221 * .. Local Arrays ..
222  INTEGER isave( 3 )
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL xerbla, zaxpy, zcopy, zlacn2, ztbmv, ztbsv
226 * ..
227 * .. Intrinsic Functions ..
228  INTRINSIC abs, dble, dimag, max, min
229 * ..
230 * .. External Functions ..
231  LOGICAL lsame
232  DOUBLE PRECISION dlamch
233  EXTERNAL lsame, dlamch
234 * ..
235 * .. Statement Functions ..
236  DOUBLE PRECISION cabs1
237 * ..
238 * .. Statement Function definitions ..
239  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZTBRFS', -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 = dlamch( 'Epsilon' )
297  safmin = dlamch( '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 zcopy( n, x( 1, j ), 1, work, 1 )
309  CALL ztbmv( uplo, trans, diag, n, kd, ab, ldab, work, 1 )
310  CALL zaxpy( 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 ZLACN2 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 zlacn2( 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 ztbsv( 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 ztbsv( 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 ZTBRFS
496 *
497  END