1*> \brief \b CTBT05
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 CTBT05( 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*       REAL               BERR( * ), FERR( * ), RESLTS( * )
20*       COMPLEX            AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
21*      $                   XACT( LDXACT, * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> CTBT05 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 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 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 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 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 REAL 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 REAL 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 REAL 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*> \date November 2011
185*
186*> \ingroup complex_lin
187*
188*  =====================================================================
189      SUBROUTINE CTBT05( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
190     $                   LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS )
191*
192*  -- LAPACK test routine (version 3.4.0) --
193*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
194*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195*     November 2011
196*
197*     .. Scalar Arguments ..
198      CHARACTER          DIAG, TRANS, UPLO
199      INTEGER            KD, LDAB, LDB, LDX, LDXACT, N, NRHS
200*     ..
201*     .. Array Arguments ..
202      REAL               BERR( * ), FERR( * ), RESLTS( * )
203      COMPLEX            AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
204     $                   XACT( LDXACT, * )
205*     ..
206*
207*  =====================================================================
208*
209*     .. Parameters ..
210      REAL               ZERO, ONE
211      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
212*     ..
213*     .. Local Scalars ..
214      LOGICAL            NOTRAN, UNIT, UPPER
215      INTEGER            I, IFU, IMAX, J, K, NZ
216      REAL               AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
217      COMPLEX            ZDUM
218*     ..
219*     .. External Functions ..
220      LOGICAL            LSAME
221      INTEGER            ICAMAX
222      REAL               SLAMCH
223      EXTERNAL           LSAME, ICAMAX, SLAMCH
224*     ..
225*     .. Intrinsic Functions ..
226      INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
227*     ..
228*     .. Statement Functions ..
229      REAL               CABS1
230*     ..
231*     .. Statement Function definitions ..
232      CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
233*     ..
234*     .. Executable Statements ..
235*
236*     Quick exit if N = 0 or NRHS = 0.
237*
238      IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
239         RESLTS( 1 ) = ZERO
240         RESLTS( 2 ) = ZERO
241         RETURN
242      END IF
243*
244      EPS = SLAMCH( 'Epsilon' )
245      UNFL = SLAMCH( 'Safe minimum' )
246      OVFL = ONE / UNFL
247      UPPER = LSAME( UPLO, 'U' )
248      NOTRAN = LSAME( TRANS, 'N' )
249      UNIT = LSAME( DIAG, 'U' )
250      NZ = MIN( KD, N-1 ) + 1
251*
252*     Test 1:  Compute the maximum of
253*        norm(X - XACT) / ( norm(X) * FERR )
254*     over all the vectors X and XACT using the infinity-norm.
255*
256      ERRBND = ZERO
257      DO 30 J = 1, NRHS
258         IMAX = ICAMAX( N, X( 1, J ), 1 )
259         XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
260         DIFF = ZERO
261         DO 10 I = 1, N
262            DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
263   10    CONTINUE
264*
265         IF( XNORM.GT.ONE ) THEN
266            GO TO 20
267         ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
268            GO TO 20
269         ELSE
270            ERRBND = ONE / EPS
271            GO TO 30
272         END IF
273*
274   20    CONTINUE
275         IF( DIFF / XNORM.LE.FERR( J ) ) THEN
276            ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
277         ELSE
278            ERRBND = ONE / EPS
279         END IF
280   30 CONTINUE
281      RESLTS( 1 ) = ERRBND
282*
283*     Test 2:  Compute the maximum of BERR / ( NZ*EPS + (*) ), where
284*     (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
285*
286      IFU = 0
287      IF( UNIT )
288     $   IFU = 1
289      DO 90 K = 1, NRHS
290         DO 80 I = 1, N
291            TMP = CABS1( B( I, K ) )
292            IF( UPPER ) THEN
293               IF( .NOT.NOTRAN ) THEN
294                  DO 40 J = MAX( I-KD, 1 ), I - IFU
295                     TMP = TMP + CABS1( AB( KD+1-I+J, I ) )*
296     $                     CABS1( X( J, K ) )
297   40             CONTINUE
298                  IF( UNIT )
299     $               TMP = TMP + CABS1( X( I, K ) )
300               ELSE
301                  IF( UNIT )
302     $               TMP = TMP + CABS1( X( I, K ) )
303                  DO 50 J = I + IFU, MIN( I+KD, N )
304                     TMP = TMP + CABS1( AB( KD+1+I-J, J ) )*
305     $                     CABS1( X( J, K ) )
306   50             CONTINUE
307               END IF
308            ELSE
309               IF( NOTRAN ) THEN
310                  DO 60 J = MAX( I-KD, 1 ), I - IFU
311                     TMP = TMP + CABS1( AB( 1+I-J, J ) )*
312     $                     CABS1( X( J, K ) )
313   60             CONTINUE
314                  IF( UNIT )
315     $               TMP = TMP + CABS1( X( I, K ) )
316               ELSE
317                  IF( UNIT )
318     $               TMP = TMP + CABS1( X( I, K ) )
319                  DO 70 J = I + IFU, MIN( I+KD, N )
320                     TMP = TMP + CABS1( AB( 1+J-I, I ) )*
321     $                     CABS1( X( J, K ) )
322   70             CONTINUE
323               END IF
324            END IF
325            IF( I.EQ.1 ) THEN
326               AXBI = TMP
327            ELSE
328               AXBI = MIN( AXBI, TMP )
329            END IF
330   80    CONTINUE
331         TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
332         IF( K.EQ.1 ) THEN
333            RESLTS( 2 ) = TMP
334         ELSE
335            RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
336         END IF
337   90 CONTINUE
338*
339      RETURN
340*
341*     End of CTBT05
342*
343      END
344