LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
strt05.f
Go to the documentation of this file.
1 *> \brief \b STRT05
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 STRT05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
12 * LDX, XACT, LDXACT, FERR, BERR, RESLTS )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER DIAG, TRANS, 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 *> STRT05 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 *> triangular 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 matrix A is upper or lower triangular.
50 *> = 'U': Upper triangular
51 *> = 'L': Lower triangular
52 *> \endverbatim
53 *>
54 *> \param[in] TRANS
55 *> \verbatim
56 *> TRANS is CHARACTER*1
57 *> Specifies the form of the system of equations.
58 *> = 'N': A * X = B (No transpose)
59 *> = 'T': A'* X = B (Transpose)
60 *> = 'C': A'* X = B (Conjugate transpose = Transpose)
61 *> \endverbatim
62 *>
63 *> \param[in] DIAG
64 *> \verbatim
65 *> DIAG is CHARACTER*1
66 *> Specifies whether or not the matrix A is unit triangular.
67 *> = 'N': Non-unit triangular
68 *> = 'U': Unit triangular
69 *> \endverbatim
70 *>
71 *> \param[in] N
72 *> \verbatim
73 *> N is INTEGER
74 *> The number of rows of the matrices X, B, and XACT, and the
75 *> order of the matrix A. N >= 0.
76 *> \endverbatim
77 *>
78 *> \param[in] NRHS
79 *> \verbatim
80 *> NRHS is INTEGER
81 *> The number of columns of the matrices X, B, and XACT.
82 *> NRHS >= 0.
83 *> \endverbatim
84 *>
85 *> \param[in] A
86 *> \verbatim
87 *> A is REAL array, dimension (LDA,N)
88 *> The triangular matrix A. If UPLO = 'U', the leading n by n
89 *> upper triangular part of the array A contains the upper
90 *> triangular matrix, and the strictly lower triangular part of
91 *> A is not referenced. If UPLO = 'L', the leading n by n lower
92 *> triangular part of the array A contains the lower triangular
93 *> matrix, and the strictly upper triangular part of A is not
94 *> referenced. If DIAG = 'U', the diagonal elements of A are
95 *> also not referenced and are assumed to be 1.
96 *> \endverbatim
97 *>
98 *> \param[in] LDA
99 *> \verbatim
100 *> LDA is INTEGER
101 *> The leading dimension of the array A. LDA >= max(1,N).
102 *> \endverbatim
103 *>
104 *> \param[in] B
105 *> \verbatim
106 *> B is REAL array, dimension (LDB,NRHS)
107 *> The right hand side vectors for the system of linear
108 *> equations.
109 *> \endverbatim
110 *>
111 *> \param[in] LDB
112 *> \verbatim
113 *> LDB is INTEGER
114 *> The leading dimension of the array B. LDB >= max(1,N).
115 *> \endverbatim
116 *>
117 *> \param[in] X
118 *> \verbatim
119 *> X is REAL array, dimension (LDX,NRHS)
120 *> The computed solution vectors. Each vector is stored as a
121 *> column of the matrix X.
122 *> \endverbatim
123 *>
124 *> \param[in] LDX
125 *> \verbatim
126 *> LDX is INTEGER
127 *> The leading dimension of the array X. LDX >= max(1,N).
128 *> \endverbatim
129 *>
130 *> \param[in] XACT
131 *> \verbatim
132 *> XACT is REAL array, dimension (LDX,NRHS)
133 *> The exact solution vectors. Each vector is stored as a
134 *> column of the matrix XACT.
135 *> \endverbatim
136 *>
137 *> \param[in] LDXACT
138 *> \verbatim
139 *> LDXACT is INTEGER
140 *> The leading dimension of the array XACT. LDXACT >= max(1,N).
141 *> \endverbatim
142 *>
143 *> \param[in] FERR
144 *> \verbatim
145 *> FERR is REAL array, dimension (NRHS)
146 *> The estimated forward error bounds for each solution vector
147 *> X. If XTRUE is the true solution, FERR bounds the magnitude
148 *> of the largest entry in (X - XTRUE) divided by the magnitude
149 *> of the largest entry in X.
150 *> \endverbatim
151 *>
152 *> \param[in] BERR
153 *> \verbatim
154 *> BERR is REAL array, dimension (NRHS)
155 *> The componentwise relative backward error of each solution
156 *> vector (i.e., the smallest relative change in any entry of A
157 *> or B that makes X an exact solution).
158 *> \endverbatim
159 *>
160 *> \param[out] RESLTS
161 *> \verbatim
162 *> RESLTS is REAL array, dimension (2)
163 *> The maximum over the NRHS solution vectors of the ratios:
164 *> RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
165 *> RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
166 *> \endverbatim
167 *
168 * Authors:
169 * ========
170 *
171 *> \author Univ. of Tennessee
172 *> \author Univ. of California Berkeley
173 *> \author Univ. of Colorado Denver
174 *> \author NAG Ltd.
175 *
176 *> \date December 2016
177 *
178 *> \ingroup single_lin
179 *
180 * =====================================================================
181  SUBROUTINE strt05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
182  $ LDX, XACT, LDXACT, FERR, BERR, RESLTS )
183 *
184 * -- LAPACK test routine (version 3.7.0) --
185 * -- LAPACK is a software package provided by Univ. of Tennessee, --
186 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187 * December 2016
188 *
189 * .. Scalar Arguments ..
190  CHARACTER DIAG, TRANS, UPLO
191  INTEGER LDA, LDB, LDX, LDXACT, N, NRHS
192 * ..
193 * .. Array Arguments ..
194  REAL A( lda, * ), B( ldb, * ), BERR( * ), FERR( * ),
195  $ reslts( * ), x( ldx, * ), xact( ldxact, * )
196 * ..
197 *
198 * =====================================================================
199 *
200 * .. Parameters ..
201  REAL ZERO, ONE
202  parameter( zero = 0.0e+0, one = 1.0e+0 )
203 * ..
204 * .. Local Scalars ..
205  LOGICAL NOTRAN, UNIT, UPPER
206  INTEGER I, IFU, IMAX, J, K
207  REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
208 * ..
209 * .. External Functions ..
210  LOGICAL LSAME
211  INTEGER ISAMAX
212  REAL SLAMCH
213  EXTERNAL lsame, isamax, slamch
214 * ..
215 * .. Intrinsic Functions ..
216  INTRINSIC abs, max, min
217 * ..
218 * .. Executable Statements ..
219 *
220 * Quick exit if N = 0 or NRHS = 0.
221 *
222  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
223  reslts( 1 ) = zero
224  reslts( 2 ) = zero
225  RETURN
226  END IF
227 *
228  eps = slamch( 'Epsilon' )
229  unfl = slamch( 'Safe minimum' )
230  ovfl = one / unfl
231  upper = lsame( uplo, 'U' )
232  notran = lsame( trans, 'N' )
233  unit = lsame( diag, 'U' )
234 *
235 * Test 1: Compute the maximum of
236 * norm(X - XACT) / ( norm(X) * FERR )
237 * over all the vectors X and XACT using the infinity-norm.
238 *
239  errbnd = zero
240  DO 30 j = 1, nrhs
241  imax = isamax( n, x( 1, j ), 1 )
242  xnorm = max( abs( x( imax, j ) ), unfl )
243  diff = zero
244  DO 10 i = 1, n
245  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
246  10 CONTINUE
247 *
248  IF( xnorm.GT.one ) THEN
249  GO TO 20
250  ELSE IF( diff.LE.ovfl*xnorm ) THEN
251  GO TO 20
252  ELSE
253  errbnd = one / eps
254  GO TO 30
255  END IF
256 *
257  20 CONTINUE
258  IF( diff / xnorm.LE.ferr( j ) ) THEN
259  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
260  ELSE
261  errbnd = one / eps
262  END IF
263  30 CONTINUE
264  reslts( 1 ) = errbnd
265 *
266 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
267 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
268 *
269  ifu = 0
270  IF( unit )
271  $ ifu = 1
272  DO 90 k = 1, nrhs
273  DO 80 i = 1, n
274  tmp = abs( b( i, k ) )
275  IF( upper ) THEN
276  IF( .NOT.notran ) THEN
277  DO 40 j = 1, i - ifu
278  tmp = tmp + abs( a( j, i ) )*abs( x( j, k ) )
279  40 CONTINUE
280  IF( unit )
281  $ tmp = tmp + abs( x( i, k ) )
282  ELSE
283  IF( unit )
284  $ tmp = tmp + abs( x( i, k ) )
285  DO 50 j = i + ifu, n
286  tmp = tmp + abs( a( i, j ) )*abs( x( j, k ) )
287  50 CONTINUE
288  END IF
289  ELSE
290  IF( notran ) THEN
291  DO 60 j = 1, i - ifu
292  tmp = tmp + abs( a( i, j ) )*abs( x( j, k ) )
293  60 CONTINUE
294  IF( unit )
295  $ tmp = tmp + abs( x( i, k ) )
296  ELSE
297  IF( unit )
298  $ tmp = tmp + abs( x( i, k ) )
299  DO 70 j = i + ifu, n
300  tmp = tmp + abs( a( j, i ) )*abs( x( j, k ) )
301  70 CONTINUE
302  END IF
303  END IF
304  IF( i.EQ.1 ) THEN
305  axbi = tmp
306  ELSE
307  axbi = min( axbi, tmp )
308  END IF
309  80 CONTINUE
310  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
311  $ max( axbi, ( n+1 )*unfl ) )
312  IF( k.EQ.1 ) THEN
313  reslts( 2 ) = tmp
314  ELSE
315  reslts( 2 ) = max( reslts( 2 ), tmp )
316  END IF
317  90 CONTINUE
318 *
319  RETURN
320 *
321 * End of STRT05
322 *
323  END
subroutine strt05(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
STRT05
Definition: strt05.f:183