LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
ctbrfs.f
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)
XERBLA
Definition: xerbla.f:60
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
ctbsv
subroutine ctbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
CTBSV
Definition: ctbsv.f:189
ctbmv
subroutine ctbmv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
CTBMV
Definition: ctbmv.f:186
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
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