LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
170 *
171 *> \ingroup double_lin
172 *
173 * =====================================================================
174  SUBROUTINE dtpt05( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
175  $ xact, ldxact, ferr, berr, reslts )
176 *
177 * -- LAPACK test routine (version 3.4.0) --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * November 2011
181 *
182 * .. Scalar Arguments ..
183  CHARACTER diag, trans, uplo
184  INTEGER ldb, ldx, ldxact, n, nrhs
185 * ..
186 * .. Array Arguments ..
187  DOUBLE PRECISION ap( * ), b( ldb, * ), berr( * ), ferr( * ),
188  $ reslts( * ), x( ldx, * ), xact( ldxact, * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  DOUBLE PRECISION zero, one
195  parameter( zero = 0.0d+0, one = 1.0d+0 )
196 * ..
197 * .. Local Scalars ..
198  LOGICAL notran, unit, upper
199  INTEGER i, ifu, imax, j, jc, k
200  DOUBLE PRECISION axbi, diff, eps, errbnd, ovfl, tmp, unfl, xnorm
201 * ..
202 * .. External Functions ..
203  LOGICAL lsame
204  INTEGER idamax
205  DOUBLE PRECISION dlamch
206  EXTERNAL lsame, idamax, dlamch
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC abs, max, min
210 * ..
211 * .. Executable Statements ..
212 *
213 * Quick exit if N = 0 or NRHS = 0.
214 *
215  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
216  reslts( 1 ) = zero
217  reslts( 2 ) = zero
218  return
219  END IF
220 *
221  eps = dlamch( 'Epsilon' )
222  unfl = dlamch( 'Safe minimum' )
223  ovfl = one / unfl
224  upper = lsame( uplo, 'U' )
225  notran = lsame( trans, 'N' )
226  unit = lsame( diag, 'U' )
227 *
228 * Test 1: Compute the maximum of
229 * norm(X - XACT) / ( norm(X) * FERR )
230 * over all the vectors X and XACT using the infinity-norm.
231 *
232  errbnd = zero
233  DO 30 j = 1, nrhs
234  imax = idamax( n, x( 1, j ), 1 )
235  xnorm = max( abs( x( imax, j ) ), unfl )
236  diff = zero
237  DO 10 i = 1, n
238  diff = max( diff, abs( x( i, j )-xact( i, j ) ) )
239  10 continue
240 *
241  IF( xnorm.GT.one ) THEN
242  go to 20
243  ELSE IF( diff.LE.ovfl*xnorm ) THEN
244  go to 20
245  ELSE
246  errbnd = one / eps
247  go to 30
248  END IF
249 *
250  20 continue
251  IF( diff / xnorm.LE.ferr( j ) ) THEN
252  errbnd = max( errbnd, ( diff / xnorm ) / ferr( j ) )
253  ELSE
254  errbnd = one / eps
255  END IF
256  30 continue
257  reslts( 1 ) = errbnd
258 *
259 * Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
260 * (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
261 *
262  ifu = 0
263  IF( unit )
264  $ ifu = 1
265  DO 90 k = 1, nrhs
266  DO 80 i = 1, n
267  tmp = abs( b( i, k ) )
268  IF( upper ) THEN
269  jc = ( ( i-1 )*i ) / 2
270  IF( .NOT.notran ) THEN
271  DO 40 j = 1, i - ifu
272  tmp = tmp + abs( ap( jc+j ) )*abs( x( j, k ) )
273  40 continue
274  IF( unit )
275  $ tmp = tmp + abs( x( i, k ) )
276  ELSE
277  jc = jc + i
278  IF( unit ) THEN
279  tmp = tmp + abs( x( i, k ) )
280  jc = jc + i
281  END IF
282  DO 50 j = i + ifu, n
283  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
284  jc = jc + j
285  50 continue
286  END IF
287  ELSE
288  IF( notran ) THEN
289  jc = i
290  DO 60 j = 1, i - ifu
291  tmp = tmp + abs( ap( jc ) )*abs( x( j, k ) )
292  jc = jc + n - j
293  60 continue
294  IF( unit )
295  $ tmp = tmp + abs( x( i, k ) )
296  ELSE
297  jc = ( i-1 )*( n-i ) + ( i*( i+1 ) ) / 2
298  IF( unit )
299  $ tmp = tmp + abs( x( i, k ) )
300  DO 70 j = i + ifu, n
301  tmp = tmp + abs( ap( jc+j-i ) )*abs( x( j, k ) )
302  70 continue
303  END IF
304  END IF
305  IF( i.EQ.1 ) THEN
306  axbi = tmp
307  ELSE
308  axbi = min( axbi, tmp )
309  END IF
310  80 continue
311  tmp = berr( k ) / ( ( n+1 )*eps+( n+1 )*unfl /
312  $ max( axbi, ( n+1 )*unfl ) )
313  IF( k.EQ.1 ) THEN
314  reslts( 2 ) = tmp
315  ELSE
316  reslts( 2 ) = max( reslts( 2 ), tmp )
317  END IF
318  90 continue
319 *
320  return
321 *
322 * End of DTPT05
323 *
324  END