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 December 2016
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.7.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*     December 2016
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
289