1*> \brief \b ZTBT05
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 ZTBT05( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
12*                          LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          DIAG, TRANS, UPLO
16*       INTEGER            KD, LDAB, LDB, LDX, LDXACT, N, NRHS
17*       ..
18*       .. Array Arguments ..
19*       DOUBLE PRECISION   BERR( * ), FERR( * ), RESLTS( * )
20*       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
21*      $                   XACT( LDXACT, * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> ZTBT05 tests the error bounds from iterative refinement for the
31*> computed solution to a system of equations A*X = B, where A is a
32*> triangular band matrix.
33*>
34*> RESLTS(1) = test of the error bound
35*>           = norm(X - XACT) / ( norm(X) * FERR )
36*>
37*> A large value is returned if this ratio is not less than one.
38*>
39*> RESLTS(2) = residual from the iterative refinement routine
40*>           = the maximum of BERR / ( NZ*EPS + (*) ), where
41*>             (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
42*>             and NZ = max. number of nonzeros in any row of A, plus 1
43*> \endverbatim
44*
45*  Arguments:
46*  ==========
47*
48*> \param[in] UPLO
49*> \verbatim
50*>          UPLO is CHARACTER*1
51*>          Specifies whether the matrix A is upper or lower triangular.
52*>          = 'U':  Upper triangular
53*>          = 'L':  Lower triangular
54*> \endverbatim
55*>
56*> \param[in] TRANS
57*> \verbatim
58*>          TRANS is CHARACTER*1
59*>          Specifies the form of the system of equations.
60*>          = 'N':  A * X = B  (No transpose)
61*>          = 'T':  A'* X = B  (Transpose)
62*>          = 'C':  A'* X = B  (Conjugate transpose = Transpose)
63*> \endverbatim
64*>
65*> \param[in] DIAG
66*> \verbatim
67*>          DIAG is CHARACTER*1
68*>          Specifies whether or not the matrix A is unit triangular.
69*>          = 'N':  Non-unit triangular
70*>          = 'U':  Unit triangular
71*> \endverbatim
72*>
73*> \param[in] N
74*> \verbatim
75*>          N is INTEGER
76*>          The number of rows of the matrices X, B, and XACT, and the
77*>          order of the matrix A.  N >= 0.
78*> \endverbatim
79*>
80*> \param[in] KD
81*> \verbatim
82*>          KD is INTEGER
83*>          The number of super-diagonals of the matrix A if UPLO = 'U',
84*>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
85*> \endverbatim
86*>
87*> \param[in] NRHS
88*> \verbatim
89*>          NRHS is INTEGER
90*>          The number of columns of the matrices X, B, and XACT.
91*>          NRHS >= 0.
92*> \endverbatim
93*>
94*> \param[in] AB
95*> \verbatim
96*>          AB is COMPLEX*16 array, dimension (LDAB,N)
97*>          The upper or lower triangular band matrix A, stored in the
98*>          first kd+1 rows of the array. The j-th column of A is stored
99*>          in the j-th column of the array AB as follows:
100*>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
101*>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
102*>          If DIAG = 'U', the diagonal elements of A are not referenced
103*>          and are assumed to be 1.
104*> \endverbatim
105*>
106*> \param[in] LDAB
107*> \verbatim
108*>          LDAB is INTEGER
109*>          The leading dimension of the array AB.  LDAB >= KD+1.
110*> \endverbatim
111*>
112*> \param[in] B
113*> \verbatim
114*>          B is COMPLEX*16 array, dimension (LDB,NRHS)
115*>          The right hand side vectors for the system of linear
116*>          equations.
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 computed solution vectors.  Each vector is stored as a
129*>          column of the matrix X.
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] XACT
139*> \verbatim
140*>          XACT is COMPLEX*16 array, dimension (LDX,NRHS)
141*>          The exact solution vectors.  Each vector is stored as a
142*>          column of the matrix XACT.
143*> \endverbatim
144*>
145*> \param[in] LDXACT
146*> \verbatim
147*>          LDXACT is INTEGER
148*>          The leading dimension of the array XACT.  LDXACT >= max(1,N).
149*> \endverbatim
150*>
151*> \param[in] FERR
152*> \verbatim
153*>          FERR is DOUBLE PRECISION array, dimension (NRHS)
154*>          The estimated forward error bounds for each solution vector
155*>          X.  If XTRUE is the true solution, FERR bounds the magnitude
156*>          of the largest entry in (X - XTRUE) divided by the magnitude
157*>          of the largest entry in X.
158*> \endverbatim
159*>
160*> \param[in] BERR
161*> \verbatim
162*>          BERR is DOUBLE PRECISION array, dimension (NRHS)
163*>          The componentwise relative backward error of each solution
164*>          vector (i.e., the smallest relative change in any entry of A
165*>          or B that makes X an exact solution).
166*> \endverbatim
167*>
168*> \param[out] RESLTS
169*> \verbatim
170*>          RESLTS is DOUBLE PRECISION array, dimension (2)
171*>          The maximum over the NRHS solution vectors of the ratios:
172*>          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
173*>          RESLTS(2) = BERR / ( NZ*EPS + (*) )
174*> \endverbatim
175*
176*  Authors:
177*  ========
178*
179*> \author Univ. of Tennessee
180*> \author Univ. of California Berkeley
181*> \author Univ. of Colorado Denver
182*> \author NAG Ltd.
183*
184*> \ingroup complex16_lin
185*
186*  =====================================================================
187      SUBROUTINE ZTBT05( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
188     $                   LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS )
189*
190*  -- LAPACK test routine --
191*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
192*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
193*
194*     .. Scalar Arguments ..
195      CHARACTER          DIAG, TRANS, UPLO
196      INTEGER            KD, LDAB, LDB, LDX, LDXACT, N, NRHS
197*     ..
198*     .. Array Arguments ..
199      DOUBLE PRECISION   BERR( * ), FERR( * ), RESLTS( * )
200      COMPLEX*16         AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
201     $                   XACT( LDXACT, * )
202*     ..
203*
204*  =====================================================================
205*
206*     .. Parameters ..
207      DOUBLE PRECISION   ZERO, ONE
208      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
209*     ..
210*     .. Local Scalars ..
211      LOGICAL            NOTRAN, UNIT, UPPER
212      INTEGER            I, IFU, IMAX, J, K, NZ
213      DOUBLE PRECISION   AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
214      COMPLEX*16         ZDUM
215*     ..
216*     .. External Functions ..
217      LOGICAL            LSAME
218      INTEGER            IZAMAX
219      DOUBLE PRECISION   DLAMCH
220      EXTERNAL           LSAME, IZAMAX, DLAMCH
221*     ..
222*     .. Intrinsic Functions ..
223      INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
224*     ..
225*     .. Statement Functions ..
226      DOUBLE PRECISION   CABS1
227*     ..
228*     .. Statement Function definitions ..
229      CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
230*     ..
231*     .. Executable Statements ..
232*
233*     Quick exit if N = 0 or NRHS = 0.
234*
235      IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
236         RESLTS( 1 ) = ZERO
237         RESLTS( 2 ) = ZERO
238         RETURN
239      END IF
240*
241      EPS = DLAMCH( 'Epsilon' )
242      UNFL = DLAMCH( 'Safe minimum' )
243      OVFL = ONE / UNFL
244      UPPER = LSAME( UPLO, 'U' )
245      NOTRAN = LSAME( TRANS, 'N' )
246      UNIT = LSAME( DIAG, 'U' )
247      NZ = MIN( KD, N-1 ) + 1
248*
249*     Test 1:  Compute the maximum of
250*        norm(X - XACT) / ( norm(X) * FERR )
251*     over all the vectors X and XACT using the infinity-norm.
252*
253      ERRBND = ZERO
254      DO 30 J = 1, NRHS
255         IMAX = IZAMAX( N, X( 1, J ), 1 )
256         XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
257         DIFF = ZERO
258         DO 10 I = 1, N
259            DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
260   10    CONTINUE
261*
262         IF( XNORM.GT.ONE ) THEN
263            GO TO 20
264         ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
265            GO TO 20
266         ELSE
267            ERRBND = ONE / EPS
268            GO TO 30
269         END IF
270*
271   20    CONTINUE
272         IF( DIFF / XNORM.LE.FERR( J ) ) THEN
273            ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
274         ELSE
275            ERRBND = ONE / EPS
276         END IF
277   30 CONTINUE
278      RESLTS( 1 ) = ERRBND
279*
280*     Test 2:  Compute the maximum of BERR / ( NZ*EPS + (*) ), where
281*     (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
282*
283      IFU = 0
284      IF( UNIT )
285     $   IFU = 1
286      DO 90 K = 1, NRHS
287         DO 80 I = 1, N
288            TMP = CABS1( B( I, K ) )
289            IF( UPPER ) THEN
290               IF( .NOT.NOTRAN ) THEN
291                  DO 40 J = MAX( I-KD, 1 ), I - IFU
292                     TMP = TMP + CABS1( AB( KD+1-I+J, I ) )*
293     $                     CABS1( X( J, K ) )
294   40             CONTINUE
295                  IF( UNIT )
296     $               TMP = TMP + CABS1( X( I, K ) )
297               ELSE
298                  IF( UNIT )
299     $               TMP = TMP + CABS1( X( I, K ) )
300                  DO 50 J = I + IFU, MIN( I+KD, N )
301                     TMP = TMP + CABS1( AB( KD+1+I-J, J ) )*
302     $                     CABS1( X( J, K ) )
303   50             CONTINUE
304               END IF
305            ELSE
306               IF( NOTRAN ) THEN
307                  DO 60 J = MAX( I-KD, 1 ), I - IFU
308                     TMP = TMP + CABS1( AB( 1+I-J, J ) )*
309     $                     CABS1( X( J, K ) )
310   60             CONTINUE
311                  IF( UNIT )
312     $               TMP = TMP + CABS1( X( I, K ) )
313               ELSE
314                  IF( UNIT )
315     $               TMP = TMP + CABS1( X( I, K ) )
316                  DO 70 J = I + IFU, MIN( I+KD, N )
317                     TMP = TMP + CABS1( AB( 1+J-I, I ) )*
318     $                     CABS1( X( J, K ) )
319   70             CONTINUE
320               END IF
321            END IF
322            IF( I.EQ.1 ) THEN
323               AXBI = TMP
324            ELSE
325               AXBI = MIN( AXBI, TMP )
326            END IF
327   80    CONTINUE
328         TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
329         IF( K.EQ.1 ) THEN
330            RESLTS( 2 ) = TMP
331         ELSE
332            RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
333         END IF
334   90 CONTINUE
335*
336      RETURN
337*
338*     End of ZTBT05
339*
340      END
341