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