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