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