LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
ctbt02.f
Go to the documentation of this file.
1 *> \brief \b CTBT02
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 CTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X,
12 * LDX, B, LDB, WORK, RWORK, RESID )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER DIAG, TRANS, UPLO
16 * INTEGER KD, LDAB, LDB, LDX, N, NRHS
17 * REAL RESID
18 * ..
19 * .. Array Arguments ..
20 * REAL RWORK( * )
21 * COMPLEX AB( LDAB, * ), B( LDB, * ), WORK( * ),
22 * $ X( LDX, * )
23 * ..
24 *
25 *
26 *> \par Purpose:
27 * =============
28 *>
29 *> \verbatim
30 *>
31 *> CTBT02 computes the residual for the computed solution to a
32 *> triangular system of linear equations A*x = b, A**T *x = b, or
33 *> A**H *x = b when A is a triangular band matrix. Here A**T denotes
34 *> the transpose of A, A**H denotes the conjugate transpose of A, and
35 *> x and b are N by NRHS matrices. The test ratio is the maximum over
36 *> the number of right hand sides of
37 *> norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
38 *> where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
39 *> \endverbatim
40 *
41 * Arguments:
42 * ==========
43 *
44 *> \param[in] UPLO
45 *> \verbatim
46 *> UPLO is CHARACTER*1
47 *> Specifies whether the matrix A is upper or lower triangular.
48 *> = 'U': Upper triangular
49 *> = 'L': Lower triangular
50 *> \endverbatim
51 *>
52 *> \param[in] TRANS
53 *> \verbatim
54 *> TRANS is CHARACTER*1
55 *> Specifies the operation applied to A.
56 *> = 'N': A *x = b (No transpose)
57 *> = 'T': A**T *x = b (Transpose)
58 *> = 'C': A**H *x = b (Conjugate transpose)
59 *> \endverbatim
60 *>
61 *> \param[in] DIAG
62 *> \verbatim
63 *> DIAG is CHARACTER*1
64 *> Specifies whether or not the matrix A is unit triangular.
65 *> = 'N': Non-unit triangular
66 *> = 'U': Unit triangular
67 *> \endverbatim
68 *>
69 *> \param[in] N
70 *> \verbatim
71 *> N is INTEGER
72 *> The order of the matrix A. N >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in] KD
76 *> \verbatim
77 *> KD is INTEGER
78 *> The number of superdiagonals or subdiagonals of the
79 *> triangular band matrix A. KD >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in] NRHS
83 *> \verbatim
84 *> NRHS is INTEGER
85 *> The number of right hand sides, i.e., the number of columns
86 *> of the matrices X and B. NRHS >= 0.
87 *> \endverbatim
88 *>
89 *> \param[in] AB
90 *> \verbatim
91 *> AB is COMPLEX array, dimension (LDA,N)
92 *> The upper or lower triangular band matrix A, stored in the
93 *> first kd+1 rows of the array. The j-th column of A is stored
94 *> in the j-th column of the array AB as follows:
95 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
96 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
97 *> \endverbatim
98 *>
99 *> \param[in] LDAB
100 *> \verbatim
101 *> LDAB is INTEGER
102 *> The leading dimension of the array AB. LDAB >= max(1,KD+1).
103 *> \endverbatim
104 *>
105 *> \param[in] X
106 *> \verbatim
107 *> X is COMPLEX array, dimension (LDX,NRHS)
108 *> The computed solution vectors for the system of linear
109 *> equations.
110 *> \endverbatim
111 *>
112 *> \param[in] LDX
113 *> \verbatim
114 *> LDX is INTEGER
115 *> The leading dimension of the array X. LDX >= max(1,N).
116 *> \endverbatim
117 *>
118 *> \param[in] B
119 *> \verbatim
120 *> B is COMPLEX array, dimension (LDB,NRHS)
121 *> The right hand side vectors for the system of linear
122 *> equations.
123 *> \endverbatim
124 *>
125 *> \param[in] LDB
126 *> \verbatim
127 *> LDB is INTEGER
128 *> The leading dimension of the array B. LDB >= max(1,N).
129 *> \endverbatim
130 *>
131 *> \param[out] WORK
132 *> \verbatim
133 *> WORK is COMPLEX array, dimension (N)
134 *> \endverbatim
135 *>
136 *> \param[out] RWORK
137 *> \verbatim
138 *> RWORK is REAL array, dimension (N)
139 *> \endverbatim
140 *>
141 *> \param[out] RESID
142 *> \verbatim
143 *> RESID is REAL
144 *> The maximum over the number of right hand sides of
145 *> norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
146 *> \endverbatim
147 *
148 * Authors:
149 * ========
150 *
151 *> \author Univ. of Tennessee
152 *> \author Univ. of California Berkeley
153 *> \author Univ. of Colorado Denver
154 *> \author NAG Ltd.
155 *
156 *> \date December 2016
157 *
158 *> \ingroup complex_lin
159 *
160 * =====================================================================
161  SUBROUTINE ctbt02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X,
162  $ LDX, B, LDB, WORK, RWORK, RESID )
163 *
164 * -- LAPACK test routine (version 3.7.0) --
165 * -- LAPACK is a software package provided by Univ. of Tennessee, --
166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167 * December 2016
168 *
169 * .. Scalar Arguments ..
170  CHARACTER DIAG, TRANS, UPLO
171  INTEGER KD, LDAB, LDB, LDX, N, NRHS
172  REAL RESID
173 * ..
174 * .. Array Arguments ..
175  REAL RWORK( * )
176  COMPLEX AB( ldab, * ), B( ldb, * ), WORK( * ),
177  $ x( ldx, * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  REAL ZERO, ONE
184  parameter( zero = 0.0e+0, one = 1.0e+0 )
185 * ..
186 * .. Local Scalars ..
187  INTEGER J
188  REAL ANORM, BNORM, EPS, XNORM
189 * ..
190 * .. External Functions ..
191  LOGICAL LSAME
192  REAL CLANTB, SCASUM, SLAMCH
193  EXTERNAL lsame, clantb, scasum, slamch
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL caxpy, ccopy, ctbmv
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC cmplx, max
200 * ..
201 * .. Executable Statements ..
202 *
203 * Quick exit if N = 0 or NRHS = 0
204 *
205  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
206  resid = zero
207  RETURN
208  END IF
209 *
210 * Compute the 1-norm of A or A'.
211 *
212  IF( lsame( trans, 'N' ) ) THEN
213  anorm = clantb( '1', uplo, diag, n, kd, ab, ldab, rwork )
214  ELSE
215  anorm = clantb( 'I', uplo, diag, n, kd, ab, ldab, rwork )
216  END IF
217 *
218 * Exit with RESID = 1/EPS if ANORM = 0.
219 *
220  eps = slamch( 'Epsilon' )
221  IF( anorm.LE.zero ) THEN
222  resid = one / eps
223  RETURN
224  END IF
225 *
226 * Compute the maximum over the number of right hand sides of
227 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
228 *
229  resid = zero
230  DO 10 j = 1, nrhs
231  CALL ccopy( n, x( 1, j ), 1, work, 1 )
232  CALL ctbmv( uplo, trans, diag, n, kd, ab, ldab, work, 1 )
233  CALL caxpy( n, cmplx( -one ), b( 1, j ), 1, work, 1 )
234  bnorm = scasum( n, work, 1 )
235  xnorm = scasum( n, x( 1, j ), 1 )
236  IF( xnorm.LE.zero ) THEN
237  resid = one / eps
238  ELSE
239  resid = max( resid, ( ( bnorm / anorm ) / xnorm ) / eps )
240  END IF
241  10 CONTINUE
242 *
243  RETURN
244 *
245 * End of CTBT02
246 *
247  END
subroutine ctbmv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
CTBMV
Definition: ctbmv.f:188
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine ctbt02(UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X, LDX, B, LDB, WORK, RWORK, RESID)
CTBT02
Definition: ctbt02.f:163
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:90