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