LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
sppt05.f
Go to the documentation of this file.
1 *> \brief \b SPPT05
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 SPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
12 * LDXACT, FERR, BERR, RESLTS )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER UPLO
16 * INTEGER LDB, LDX, LDXACT, N, NRHS
17 * ..
18 * .. Array Arguments ..
19 * REAL AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
20 * $ RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> SPPT05 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 matrix in packed storage format.
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] AP
70 *> \verbatim
71 *> AP is REAL array, dimension (N*(N+1)/2)
72 *> The upper or lower triangle of the symmetric matrix A, packed
73 *> columnwise in a linear array. The j-th column of A is stored
74 *> in the array AP as follows:
75 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
76 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
77 *> \endverbatim
78 *>
79 *> \param[in] B
80 *> \verbatim
81 *> B is REAL array, dimension (LDB,NRHS)
82 *> The right hand side vectors for the system of linear
83 *> equations.
84 *> \endverbatim
85 *>
86 *> \param[in] LDB
87 *> \verbatim
88 *> LDB is INTEGER
89 *> The leading dimension of the array B. LDB >= max(1,N).
90 *> \endverbatim
91 *>
92 *> \param[in] X
93 *> \verbatim
94 *> X is REAL array, dimension (LDX,NRHS)
95 *> The computed solution vectors. Each vector is stored as a
96 *> column of the matrix X.
97 *> \endverbatim
98 *>
99 *> \param[in] LDX
100 *> \verbatim
101 *> LDX is INTEGER
102 *> The leading dimension of the array X. LDX >= max(1,N).
103 *> \endverbatim
104 *>
105 *> \param[in] XACT
106 *> \verbatim
107 *> XACT is REAL array, dimension (LDX,NRHS)
108 *> The exact solution vectors. Each vector is stored as a
109 *> column of the matrix XACT.
110 *> \endverbatim
111 *>
112 *> \param[in] LDXACT
113 *> \verbatim
114 *> LDXACT is INTEGER
115 *> The leading dimension of the array XACT. LDXACT >= max(1,N).
116 *> \endverbatim
117 *>
118 *> \param[in] FERR
119 *> \verbatim
120 *> FERR is REAL array, dimension (NRHS)
121 *> The estimated forward error bounds for each solution vector
122 *> X. If XTRUE is the true solution, FERR bounds the magnitude
123 *> of the largest entry in (X - XTRUE) divided by the magnitude
124 *> of the largest entry in X.
125 *> \endverbatim
126 *>
127 *> \param[in] BERR
128 *> \verbatim
129 *> BERR is REAL array, dimension (NRHS)
130 *> The componentwise relative backward error of each solution
131 *> vector (i.e., the smallest relative change in any entry of A
132 *> or B that makes X an exact solution).
133 *> \endverbatim
134 *>
135 *> \param[out] RESLTS
136 *> \verbatim
137 *> RESLTS is REAL array, dimension (2)
138 *> The maximum over the NRHS solution vectors of the ratios:
139 *> RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
140 *> RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
141 *> \endverbatim
142 *
143 * Authors:
144 * ========
145 *
146 *> \author Univ. of Tennessee
147 *> \author Univ. of California Berkeley
148 *> \author Univ. of Colorado Denver
149 *> \author NAG Ltd.
150 *
151 *> \ingroup single_lin
152 *
153 * =====================================================================
154  SUBROUTINE sppt05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT,
155  $ LDXACT, FERR, BERR, RESLTS )
156 *
157 * -- LAPACK test routine --
158 * -- LAPACK is a software package provided by Univ. of Tennessee, --
159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160 *
161 * .. Scalar Arguments ..
162  CHARACTER UPLO
163  INTEGER LDB, LDX, LDXACT, N, NRHS
164 * ..
165 * .. Array Arguments ..
166  REAL AP( * ), B( LDB, * ), BERR( * ), FERR( * ),
167  $ reslts( * ), x( ldx, * ), xact( ldxact, * )
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..
173  REAL ZERO, ONE
174  parameter( zero = 0.0e+0, one = 1.0e+0 )
175 * ..
176 * .. Local Scalars ..
177  LOGICAL UPPER
178  INTEGER I, IMAX, J, JC, K
179  REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
180 * ..
181 * .. External Functions ..
182  LOGICAL LSAME
183  INTEGER ISAMAX
184  REAL SLAMCH
185  EXTERNAL lsame, isamax, slamch
186 * ..
187 * .. Intrinsic Functions ..
188  INTRINSIC abs, max, min
189 * ..
190 * .. Executable Statements ..
191 *
192 * Quick exit if N = 0 or NRHS = 0.
193 *
194  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
195  reslts( 1 ) = zero
196  reslts( 2 ) = zero
197  RETURN
198  END IF
199 *
200  eps = slamch( 'Epsilon' )
201  unfl = slamch( 'Safe minimum' )
202  ovfl = one / unfl
203  upper = lsame( uplo, 'U' )
204 *
205 * Test 1: Compute the maximum of
206 * norm(X - XACT) / ( norm(X) * FERR )
207 * over all the vectors X and XACT using the infinity-norm.
208 *
209  errbnd = zero
210  DO 30 j = 1, nrhs
211  imax = isamax( n, x( 1, j ), 1 )
212  xnorm = max( abs( x( imax, j ) ), unfl )
213  diff = zero
214  DO 10 i = 1, n
215  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
216  10 CONTINUE
217 *
218  IF( xnorm.GT.one ) THEN
219  GO TO 20
220  ELSE IF( diff.LE.ovfl*xnorm ) THEN
221  GO TO 20
222  ELSE
223  errbnd = one / eps
224  GO TO 30
225  END IF
226 *
227  20 CONTINUE
228  IF( diff / xnorm.LE.ferr( j ) ) THEN
229  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
230  ELSE
231  errbnd = one / eps
232  END IF
233  30 CONTINUE
234  reslts( 1 ) = errbnd
235 *
236 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
237 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
238 *
239  DO 90 k = 1, nrhs
240  DO 80 i = 1, n
241  tmp = abs( b( i, k ) )
242  IF( upper ) THEN
243  jc = ( ( i-1 )*i ) / 2
244  DO 40 j = 1, i
245  tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
246  40 CONTINUE
247  jc = jc + i
248  DO 50 j = i + 1, n
249  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
250  jc = jc + j
251  50 CONTINUE
252  ELSE
253  jc = i
254  DO 60 j = 1, i - 1
255  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
256  jc = jc + n - j
257  60 CONTINUE
258  DO 70 j = i, n
259  tmp = tmp + abs( ap( jc+j-i ) )*abs( x( j, k ) )
260  70 CONTINUE
261  END IF
262  IF( i.EQ.1 ) THEN
263  axbi = tmp
264  ELSE
265  axbi = min( axbi, tmp )
266  END IF
267  80 CONTINUE
268  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
269  $ max( axbi, ( n+1 )*unfl ) )
270  IF( k.EQ.1 ) THEN
271  reslts( 2 ) = tmp
272  ELSE
273  reslts( 2 ) = max( reslts( 2 ), tmp )
274  END IF
275  90 CONTINUE
276 *
277  RETURN
278 *
279 * End of SPPT05
280 *
281  END
subroutine sppt05(UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
SPPT05
Definition: sppt05.f:156