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