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