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