LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dget39.f
Go to the documentation of this file.
1*> \brief \b DGET39
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 DGET39( RMAX, LMAX, NINFO, KNT )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT, LMAX, NINFO
15* DOUBLE PRECISION RMAX
16* ..
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> DGET39 tests DLAQTR, a routine for solving the real or
25*> special complex quasi upper triangular system
26*>
27*> op(T)*p = scale*c,
28*> or
29*> op(T + iB)*(p+iq) = scale*(c+id),
30*>
31*> in real arithmetic. T is upper quasi-triangular.
32*> If it is complex, then the first diagonal block of T must be
33*> 1 by 1, B has the special structure
34*>
35*> B = [ b(1) b(2) ... b(n) ]
36*> [ w ]
37*> [ w ]
38*> [ . ]
39*> [ w ]
40*>
41*> op(A) = A or A', where A' denotes the conjugate transpose of
42*> the matrix A.
43*>
44*> On input, X = [ c ]. On output, X = [ p ].
45*> [ d ] [ q ]
46*>
47*> Scale is an output less than or equal to 1, chosen to avoid
48*> overflow in X.
49*> This subroutine is specially designed for the condition number
50*> estimation in the eigenproblem routine DTRSNA.
51*>
52*> The test code verifies that the following residual is order 1:
53*>
54*> ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
55*> -----------------------------------------
56*> max(ulp*(||T||+||B||)*(||x1||+||x2||),
57*> (||T||+||B||)*smlnum/ulp,
58*> smlnum)
59*>
60*> (The (||T||+||B||)*smlnum/ulp term accounts for possible
61*> (gradual or nongradual) underflow in x1 and x2.)
62*> \endverbatim
63*
64* Arguments:
65* ==========
66*
67*> \param[out] RMAX
68*> \verbatim
69*> RMAX is DOUBLE PRECISION
70*> Value of the largest test ratio.
71*> \endverbatim
72*>
73*> \param[out] LMAX
74*> \verbatim
75*> LMAX is INTEGER
76*> Example number where largest test ratio achieved.
77*> \endverbatim
78*>
79*> \param[out] NINFO
80*> \verbatim
81*> NINFO is INTEGER
82*> Number of examples where INFO is nonzero.
83*> \endverbatim
84*>
85*> \param[out] KNT
86*> \verbatim
87*> KNT is INTEGER
88*> Total number of examples tested.
89*> \endverbatim
90*
91* Authors:
92* ========
93*
94*> \author Univ. of Tennessee
95*> \author Univ. of California Berkeley
96*> \author Univ. of Colorado Denver
97*> \author NAG Ltd.
98*
99*> \ingroup double_eig
100*
101* =====================================================================
102 SUBROUTINE dget39( RMAX, LMAX, NINFO, KNT )
103*
104* -- LAPACK test routine --
105* -- LAPACK is a software package provided by Univ. of Tennessee, --
106* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
107*
108* .. Scalar Arguments ..
109 INTEGER KNT, LMAX, NINFO
110 DOUBLE PRECISION RMAX
111* ..
112*
113* =====================================================================
114*
115* .. Parameters ..
116 INTEGER LDT, LDT2
117 parameter( ldt = 10, ldt2 = 2*ldt )
118 DOUBLE PRECISION ZERO, ONE
119 parameter( zero = 0.0d0, one = 1.0d0 )
120* ..
121* .. Local Scalars ..
122 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
123 $ NDIM
124 DOUBLE PRECISION BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
125 $ SCALE, SMLNUM, W, XNORM
126* ..
127* .. External Functions ..
128 INTEGER IDAMAX
129 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
130 EXTERNAL idamax, dasum, ddot, dlamch, dlange
131* ..
132* .. External Subroutines ..
133 EXTERNAL dcopy, dgemv, dlaqtr
134* ..
135* .. Intrinsic Functions ..
136 INTRINSIC abs, cos, dble, max, sin, sqrt
137* ..
138* .. Local Arrays ..
139 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
140 DOUBLE PRECISION B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
141 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
142 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
143* ..
144* .. Data statements ..
145 DATA idim / 4, 5*5 /
146 DATA ival / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
147 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
148 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
149 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
150 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
151 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
152 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
153 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
154* ..
155* .. Executable Statements ..
156*
157* Get machine parameters
158*
159 eps = dlamch( 'P' )
160 smlnum = dlamch( 'S' )
161 bignum = one / smlnum
162*
163* Set up test case parameters
164*
165 vm1( 1 ) = one
166 vm1( 2 ) = sqrt( smlnum )
167 vm1( 3 ) = sqrt( vm1( 2 ) )
168 vm1( 4 ) = sqrt( bignum )
169 vm1( 5 ) = sqrt( vm1( 4 ) )
170*
171 vm2( 1 ) = one
172 vm2( 2 ) = sqrt( smlnum )
173 vm2( 3 ) = sqrt( vm2( 2 ) )
174 vm2( 4 ) = sqrt( bignum )
175 vm2( 5 ) = sqrt( vm2( 4 ) )
176*
177 vm3( 1 ) = one
178 vm3( 2 ) = sqrt( smlnum )
179 vm3( 3 ) = sqrt( vm3( 2 ) )
180 vm3( 4 ) = sqrt( bignum )
181 vm3( 5 ) = sqrt( vm3( 4 ) )
182*
183 vm4( 1 ) = one
184 vm4( 2 ) = sqrt( smlnum )
185 vm4( 3 ) = sqrt( vm4( 2 ) )
186 vm4( 4 ) = sqrt( bignum )
187 vm4( 5 ) = sqrt( vm4( 4 ) )
188*
189 vm5( 1 ) = one
190 vm5( 2 ) = eps
191 vm5( 3 ) = sqrt( smlnum )
192*
193* Initialization
194*
195 knt = 0
196 rmax = zero
197 ninfo = 0
198 smlnum = smlnum / eps
199*
200* Begin test loop
201*
202 DO 140 ivm5 = 1, 3
203 DO 130 ivm4 = 1, 5
204 DO 120 ivm3 = 1, 5
205 DO 110 ivm2 = 1, 5
206 DO 100 ivm1 = 1, 5
207 DO 90 ndim = 1, 6
208*
209 n = idim( ndim )
210 DO 20 i = 1, n
211 DO 10 j = 1, n
212 t( i, j ) = dble( ival( i, j, ndim ) )*
213 $ vm1( ivm1 )
214 IF( i.GE.j )
215 $ t( i, j ) = t( i, j )*vm5( ivm5 )
216 10 CONTINUE
217 20 CONTINUE
218*
219 w = one*vm2( ivm2 )
220*
221 DO 30 i = 1, n
222 b( i ) = cos( dble( i ) )*vm3( ivm3 )
223 30 CONTINUE
224*
225 DO 40 i = 1, 2*n
226 d( i ) = sin( dble( i ) )*vm4( ivm4 )
227 40 CONTINUE
228*
229 norm = dlange( '1', n, n, t, ldt, work )
230 k = idamax( n, b, 1 )
231 normtb = norm + abs( b( k ) ) + abs( w )
232*
233 CALL dcopy( n, d, 1, x, 1 )
234 knt = knt + 1
235 CALL dlaqtr( .false., .true., n, t, ldt, dum,
236 $ dumm, scale, x, work, info )
237 IF( info.NE.0 )
238 $ ninfo = ninfo + 1
239*
240* || T*x - scale*d || /
241* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
242*
243 CALL dcopy( n, d, 1, y, 1 )
244 CALL dgemv( 'No transpose', n, n, one, t, ldt,
245 $ x, 1, -scale, y, 1 )
246 xnorm = dasum( n, x, 1 )
247 resid = dasum( n, y, 1 )
248 domin = max( smlnum, ( smlnum / eps )*norm,
249 $ ( norm*eps )*xnorm )
250 resid = resid / domin
251 IF( resid.GT.rmax ) THEN
252 rmax = resid
253 lmax = knt
254 END IF
255*
256 CALL dcopy( n, d, 1, x, 1 )
257 knt = knt + 1
258 CALL dlaqtr( .true., .true., n, t, ldt, dum,
259 $ dumm, scale, x, work, info )
260 IF( info.NE.0 )
261 $ ninfo = ninfo + 1
262*
263* || T*x - scale*d || /
264* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
265*
266 CALL dcopy( n, d, 1, y, 1 )
267 CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
268 $ 1, -scale, y, 1 )
269 xnorm = dasum( n, x, 1 )
270 resid = dasum( n, y, 1 )
271 domin = max( smlnum, ( smlnum / eps )*norm,
272 $ ( norm*eps )*xnorm )
273 resid = resid / domin
274 IF( resid.GT.rmax ) THEN
275 rmax = resid
276 lmax = knt
277 END IF
278*
279 CALL dcopy( 2*n, d, 1, x, 1 )
280 knt = knt + 1
281 CALL dlaqtr( .false., .false., n, t, ldt, b, w,
282 $ scale, x, work, info )
283 IF( info.NE.0 )
284 $ ninfo = ninfo + 1
285*
286* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
287* max(ulp*(||T||+||B||)*(||x1||+||x2||),
288* smlnum/ulp * (||T||+||B||), smlnum )
289*
290*
291 CALL dcopy( 2*n, d, 1, y, 1 )
292 y( 1 ) = ddot( n, b, 1, x( 1+n ), 1 ) +
293 $ scale*y( 1 )
294 DO 50 i = 2, n
295 y( i ) = w*x( i+n ) + scale*y( i )
296 50 CONTINUE
297 CALL dgemv( 'No transpose', n, n, one, t, ldt,
298 $ x, 1, -one, y, 1 )
299*
300 y( 1+n ) = ddot( n, b, 1, x, 1 ) -
301 $ scale*y( 1+n )
302 DO 60 i = 2, n
303 y( i+n ) = w*x( i ) - scale*y( i+n )
304 60 CONTINUE
305 CALL dgemv( 'No transpose', n, n, one, t, ldt,
306 $ x( 1+n ), 1, one, y( 1+n ), 1 )
307*
308 resid = dasum( 2*n, y, 1 )
309 domin = max( smlnum, ( smlnum / eps )*normtb,
310 $ eps*( normtb*dasum( 2*n, x, 1 ) ) )
311 resid = resid / domin
312 IF( resid.GT.rmax ) THEN
313 rmax = resid
314 lmax = knt
315 END IF
316*
317 CALL dcopy( 2*n, d, 1, x, 1 )
318 knt = knt + 1
319 CALL dlaqtr( .true., .false., n, t, ldt, b, w,
320 $ scale, x, work, info )
321 IF( info.NE.0 )
322 $ ninfo = ninfo + 1
323*
324* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
325* max(ulp*(||T||+||B||)*(||x1||+||x2||),
326* smlnum/ulp * (||T||+||B||), smlnum )
327*
328 CALL dcopy( 2*n, d, 1, y, 1 )
329 y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
330 DO 70 i = 2, n
331 y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
332 $ scale*y( i )
333 70 CONTINUE
334 CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
335 $ 1, one, y, 1 )
336*
337 y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
338 DO 80 i = 2, n
339 y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
340 $ scale*y( i+n )
341 80 CONTINUE
342 CALL dgemv( 'Transpose', n, n, one, t, ldt,
343 $ x( 1+n ), 1, -one, y( 1+n ), 1 )
344*
345 resid = dasum( 2*n, y, 1 )
346 domin = max( smlnum, ( smlnum / eps )*normtb,
347 $ eps*( normtb*dasum( 2*n, x, 1 ) ) )
348 resid = resid / domin
349 IF( resid.GT.rmax ) THEN
350 rmax = resid
351 lmax = knt
352 END IF
353*
354 90 CONTINUE
355 100 CONTINUE
356 110 CONTINUE
357 120 CONTINUE
358 130 CONTINUE
359 140 CONTINUE
360*
361 RETURN
362*
363* End of DGET39
364*
365 END
dget39
subroutine dget39(rmax, lmax, ninfo, knt)
DGET39
Definition dget39.f:103
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dlaqtr
subroutine dlaqtr(ltran, lreal, n, t, ldt, b, w, scale, x, work, info)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
Definition dlaqtr.f:165