LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ztbt03.f
Go to the documentation of this file.
1 *> \brief \b ZTBT03
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZTBT03( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB,
12 * SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK,
13 * RESID )
14 *
15 * .. Scalar Arguments ..
16 * CHARACTER DIAG, TRANS, UPLO
17 * INTEGER KD, LDAB, LDB, LDX, N, NRHS
18 * DOUBLE PRECISION RESID, SCALE, TSCAL
19 * ..
20 * .. Array Arguments ..
21 * DOUBLE PRECISION CNORM( * )
22 * COMPLEX*16 AB( LDAB, * ), B( LDB, * ), WORK( * ),
23 * $ X( LDX, * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> ZTBT03 computes the residual for the solution to a scaled triangular
33 *> system of equations A*x = s*b, A**T *x = s*b, or A**H *x = s*b
34 *> when A is a triangular band matrix. Here A**T denotes the transpose
35 *> of A, A**H denotes the conjugate transpose of A, s is a scalar, and
36 *> x and b are N by NRHS matrices. The test ratio is the maximum over
37 *> the number of right hand sides of
38 *> norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
39 *> where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
40 *> \endverbatim
41 *
42 * Arguments:
43 * ==========
44 *
45 *> \param[in] UPLO
46 *> \verbatim
47 *> UPLO is CHARACTER*1
48 *> Specifies whether the matrix A is upper or lower triangular.
49 *> = 'U': Upper triangular
50 *> = 'L': Lower triangular
51 *> \endverbatim
52 *>
53 *> \param[in] TRANS
54 *> \verbatim
55 *> TRANS is CHARACTER*1
56 *> Specifies the operation applied to A.
57 *> = 'N': A *x = s*b (No transpose)
58 *> = 'T': A**T *x = s*b (Transpose)
59 *> = 'C': A**H *x = s*b (Conjugate transpose)
60 *> \endverbatim
61 *>
62 *> \param[in] DIAG
63 *> \verbatim
64 *> DIAG is CHARACTER*1
65 *> Specifies whether or not the matrix A is unit triangular.
66 *> = 'N': Non-unit triangular
67 *> = 'U': Unit triangular
68 *> \endverbatim
69 *>
70 *> \param[in] N
71 *> \verbatim
72 *> N is INTEGER
73 *> The order of the matrix A. N >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in] KD
77 *> \verbatim
78 *> KD is INTEGER
79 *> The number of superdiagonals or subdiagonals of the
80 *> triangular band matrix A. KD >= 0.
81 *> \endverbatim
82 *>
83 *> \param[in] NRHS
84 *> \verbatim
85 *> NRHS is INTEGER
86 *> The number of right hand sides, i.e., the number of columns
87 *> of the matrices X and B. NRHS >= 0.
88 *> \endverbatim
89 *>
90 *> \param[in] AB
91 *> \verbatim
92 *> AB is COMPLEX*16 array, dimension (LDAB,N)
93 *> The upper or lower triangular band matrix A, stored in the
94 *> first kd+1 rows of the array. The j-th column of A is stored
95 *> in the j-th column of the array AB as follows:
96 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
97 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
98 *> \endverbatim
99 *>
100 *> \param[in] LDAB
101 *> \verbatim
102 *> LDAB is INTEGER
103 *> The leading dimension of the array AB. LDAB >= KD+1.
104 *> \endverbatim
105 *>
106 *> \param[in] SCALE
107 *> \verbatim
108 *> SCALE is DOUBLE PRECISION
109 *> The scaling factor s used in solving the triangular system.
110 *> \endverbatim
111 *>
112 *> \param[in] CNORM
113 *> \verbatim
114 *> CNORM is DOUBLE PRECISION array, dimension (N)
115 *> The 1-norms of the columns of A, not counting the diagonal.
116 *> \endverbatim
117 *>
118 *> \param[in] TSCAL
119 *> \verbatim
120 *> TSCAL is DOUBLE PRECISION
121 *> The scaling factor used in computing the 1-norms in CNORM.
122 *> CNORM actually contains the column norms of TSCAL*A.
123 *> \endverbatim
124 *>
125 *> \param[in] X
126 *> \verbatim
127 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
128 *> The computed solution vectors for the system of linear
129 *> equations.
130 *> \endverbatim
131 *>
132 *> \param[in] LDX
133 *> \verbatim
134 *> LDX is INTEGER
135 *> The leading dimension of the array X. LDX >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[in] B
139 *> \verbatim
140 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
141 *> The right hand side vectors for the system of linear
142 *> equations.
143 *> \endverbatim
144 *>
145 *> \param[in] LDB
146 *> \verbatim
147 *> LDB is INTEGER
148 *> The leading dimension of the array B. LDB >= max(1,N).
149 *> \endverbatim
150 *>
151 *> \param[out] WORK
152 *> \verbatim
153 *> WORK is COMPLEX*16 array, dimension (N)
154 *> \endverbatim
155 *>
156 *> \param[out] RESID
157 *> \verbatim
158 *> RESID is DOUBLE PRECISION
159 *> The maximum over the number of right hand sides of
160 *> norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
161 *> \endverbatim
162 *
163 * Authors:
164 * ========
165 *
166 *> \author Univ. of Tennessee
167 *> \author Univ. of California Berkeley
168 *> \author Univ. of Colorado Denver
169 *> \author NAG Ltd.
170 *
171 *> \date November 2011
172 *
173 *> \ingroup complex16_lin
174 *
175 * =====================================================================
176  SUBROUTINE ztbt03( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB,
177  $ scale, cnorm, tscal, x, ldx, b, ldb, work,
178  $ resid )
179 *
180 * -- LAPACK test routine (version 3.4.0) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * November 2011
184 *
185 * .. Scalar Arguments ..
186  CHARACTER diag, trans, uplo
187  INTEGER kd, ldab, ldb, ldx, n, nrhs
188  DOUBLE PRECISION resid, scale, tscal
189 * ..
190 * .. Array Arguments ..
191  DOUBLE PRECISION cnorm( * )
192  COMPLEX*16 ab( ldab, * ), b( ldb, * ), work( * ),
193  $ x( ldx, * )
194 * ..
195 *
196 * =====================================================================
197 *
198 *
199 * .. Parameters ..
200  DOUBLE PRECISION one, zero
201  parameter( one = 1.0d+0, zero = 0.0d+0 )
202 * ..
203 * .. Local Scalars ..
204  INTEGER ix, j
205  DOUBLE PRECISION eps, err, smlnum, tnorm, xnorm, xscal
206 * ..
207 * .. External Functions ..
208  LOGICAL lsame
209  INTEGER izamax
210  DOUBLE PRECISION dlamch
211  EXTERNAL lsame, izamax, dlamch
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL zaxpy, zcopy, zdscal, ztbmv
215 * ..
216 * .. Intrinsic Functions ..
217  INTRINSIC abs, dble, dcmplx, max
218 * ..
219 * .. Executable Statements ..
220 *
221 * Quick exit if N = 0
222 *
223  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
224  resid = zero
225  return
226  END IF
227  eps = dlamch( 'Epsilon' )
228  smlnum = dlamch( 'Safe minimum' )
229 *
230 * Compute the norm of the triangular matrix A using the column
231 * norms already computed by ZLATBS.
232 *
233  tnorm = zero
234  IF( lsame( diag, 'N' ) ) THEN
235  IF( lsame( uplo, 'U' ) ) THEN
236  DO 10 j = 1, n
237  tnorm = max( tnorm, tscal*abs( ab( kd+1, j ) )+
238  $ cnorm( j ) )
239  10 continue
240  ELSE
241  DO 20 j = 1, n
242  tnorm = max( tnorm, tscal*abs( ab( 1, j ) )+cnorm( j ) )
243  20 continue
244  END IF
245  ELSE
246  DO 30 j = 1, n
247  tnorm = max( tnorm, tscal+cnorm( j ) )
248  30 continue
249  END IF
250 *
251 * Compute the maximum over the number of right hand sides of
252 * norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
253 *
254  resid = zero
255  DO 40 j = 1, nrhs
256  CALL zcopy( n, x( 1, j ), 1, work, 1 )
257  ix = izamax( n, work, 1 )
258  xnorm = max( one, abs( x( ix, j ) ) )
259  xscal = ( one / xnorm ) / dble( kd+1 )
260  CALL zdscal( n, xscal, work, 1 )
261  CALL ztbmv( uplo, trans, diag, n, kd, ab, ldab, work, 1 )
262  CALL zaxpy( n, dcmplx( -scale*xscal ), b( 1, j ), 1, work, 1 )
263  ix = izamax( n, work, 1 )
264  err = tscal*abs( work( ix ) )
265  ix = izamax( n, x( 1, j ), 1 )
266  xnorm = abs( x( ix, j ) )
267  IF( err*smlnum.LE.xnorm ) THEN
268  IF( xnorm.GT.zero )
269  $ err = err / xnorm
270  ELSE
271  IF( err.GT.zero )
272  $ err = one / eps
273  END IF
274  IF( err*smlnum.LE.tnorm ) THEN
275  IF( tnorm.GT.zero )
276  $ err = err / tnorm
277  ELSE
278  IF( err.GT.zero )
279  $ err = one / eps
280  END IF
281  resid = max( resid, err )
282  40 continue
283 *
284  return
285 *
286 * End of ZTBT03
287 *
288  END