LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
dpot05.f
Go to the documentation of this file.
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
subroutine dpot05(UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
DPOT05
Definition: dpot05.f:166