1*> \brief \b DPOT05
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 DPOT05( 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*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
20*      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> DPOT05 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*> \date November 2011
160*
161*> \ingroup double_lin
162*
163*  =====================================================================
164      SUBROUTINE DPOT05( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, XACT,
165     $                   LDXACT, FERR, BERR, RESLTS )
166*
167*  -- LAPACK test routine (version 3.4.0) --
168*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
169*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170*     November 2011
171*
172*     .. Scalar Arguments ..
173      CHARACTER          UPLO
174      INTEGER            LDA, LDB, LDX, LDXACT, N, NRHS
175*     ..
176*     .. Array Arguments ..
177      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
178     $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
179*     ..
180*
181*  =====================================================================
182*
183*     .. Parameters ..
184      DOUBLE PRECISION   ZERO, ONE
185      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
186*     ..
187*     .. Local Scalars ..
188      LOGICAL            UPPER
189      INTEGER            I, IMAX, J, K
190      DOUBLE PRECISION   AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
191*     ..
192*     .. External Functions ..
193      LOGICAL            LSAME
194      INTEGER            IDAMAX
195      DOUBLE PRECISION   DLAMCH
196      EXTERNAL           LSAME, IDAMAX, DLAMCH
197*     ..
198*     .. Intrinsic Functions ..
199      INTRINSIC          ABS, MAX, MIN
200*     ..
201*     .. Executable Statements ..
202*
203*     Quick exit if N = 0 or NRHS = 0.
204*
205      IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
206         RESLTS( 1 ) = ZERO
207         RESLTS( 2 ) = ZERO
208         RETURN
209      END IF
210*
211      EPS = DLAMCH( 'Epsilon' )
212      UNFL = DLAMCH( 'Safe minimum' )
213      OVFL = ONE / UNFL
214      UPPER = LSAME( UPLO, 'U' )
215*
216*     Test 1:  Compute the maximum of
217*        norm(X - XACT) / ( norm(X) * FERR )
218*     over all the vectors X and XACT using the infinity-norm.
219*
220      ERRBND = ZERO
221      DO 30 J = 1, NRHS
222         IMAX = IDAMAX( N, X( 1, J ), 1 )
223         XNORM = MAX( ABS( X( IMAX, J ) ), UNFL )
224         DIFF = ZERO
225         DO 10 I = 1, N
226            DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
227   10    CONTINUE
228*
229         IF( XNORM.GT.ONE ) THEN
230            GO TO 20
231         ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
232            GO TO 20
233         ELSE
234            ERRBND = ONE / EPS
235            GO TO 30
236         END IF
237*
238   20    CONTINUE
239         IF( DIFF / XNORM.LE.FERR( J ) ) THEN
240            ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
241         ELSE
242            ERRBND = ONE / EPS
243         END IF
244   30 CONTINUE
245      RESLTS( 1 ) = ERRBND
246*
247*     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
248*     (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
249*
250      DO 90 K = 1, NRHS
251         DO 80 I = 1, N
252            TMP = ABS( B( I, K ) )
253            IF( UPPER ) THEN
254               DO 40 J = 1, I
255                  TMP = TMP + ABS( A( J, I ) )*ABS( X( J, K ) )
256   40          CONTINUE
257               DO 50 J = I + 1, N
258                  TMP = TMP + ABS( A( I, J ) )*ABS( X( J, K ) )
259   50          CONTINUE
260            ELSE
261               DO 60 J = 1, I - 1
262                  TMP = TMP + ABS( A( I, J ) )*ABS( X( J, K ) )
263   60          CONTINUE
264               DO 70 J = I, N
265                  TMP = TMP + ABS( A( J, I ) )*ABS( X( J, K ) )
266   70          CONTINUE
267            END IF
268            IF( I.EQ.1 ) THEN
269               AXBI = TMP
270            ELSE
271               AXBI = MIN( AXBI, TMP )
272            END IF
273   80    CONTINUE
274         TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL /
275     $         MAX( AXBI, ( N+1 )*UNFL ) )
276         IF( K.EQ.1 ) THEN
277            RESLTS( 2 ) = TMP
278         ELSE
279            RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
280         END IF
281   90 CONTINUE
282*
283      RETURN
284*
285*     End of DPOT05
286*
287      END
288