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