LAPACK 3.11.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, dlabad, 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 CALL dlabad( smlnum, bignum )
163*
164* Set up test case parameters
165*
166 vm1( 1 ) = one
167 vm1( 2 ) = sqrt( smlnum )
168 vm1( 3 ) = sqrt( vm1( 2 ) )
169 vm1( 4 ) = sqrt( bignum )
170 vm1( 5 ) = sqrt( vm1( 4 ) )
171*
172 vm2( 1 ) = one
173 vm2( 2 ) = sqrt( smlnum )
174 vm2( 3 ) = sqrt( vm2( 2 ) )
175 vm2( 4 ) = sqrt( bignum )
176 vm2( 5 ) = sqrt( vm2( 4 ) )
177*
178 vm3( 1 ) = one
179 vm3( 2 ) = sqrt( smlnum )
180 vm3( 3 ) = sqrt( vm3( 2 ) )
181 vm3( 4 ) = sqrt( bignum )
182 vm3( 5 ) = sqrt( vm3( 4 ) )
183*
184 vm4( 1 ) = one
185 vm4( 2 ) = sqrt( smlnum )
186 vm4( 3 ) = sqrt( vm4( 2 ) )
187 vm4( 4 ) = sqrt( bignum )
188 vm4( 5 ) = sqrt( vm4( 4 ) )
189*
190 vm5( 1 ) = one
191 vm5( 2 ) = eps
192 vm5( 3 ) = sqrt( smlnum )
193*
194* Initialization
195*
196 knt = 0
197 rmax = zero
198 ninfo = 0
199 smlnum = smlnum / eps
200*
201* Begin test loop
202*
203 DO 140 ivm5 = 1, 3
204 DO 130 ivm4 = 1, 5
205 DO 120 ivm3 = 1, 5
206 DO 110 ivm2 = 1, 5
207 DO 100 ivm1 = 1, 5
208 DO 90 ndim = 1, 6
209*
210 n = idim( ndim )
211 DO 20 i = 1, n
212 DO 10 j = 1, n
213 t( i, j ) = dble( ival( i, j, ndim ) )*
214 $ vm1( ivm1 )
215 IF( i.GE.j )
216 $ t( i, j ) = t( i, j )*vm5( ivm5 )
217 10 CONTINUE
218 20 CONTINUE
219*
220 w = one*vm2( ivm2 )
221*
222 DO 30 i = 1, n
223 b( i ) = cos( dble( i ) )*vm3( ivm3 )
224 30 CONTINUE
225*
226 DO 40 i = 1, 2*n
227 d( i ) = sin( dble( i ) )*vm4( ivm4 )
228 40 CONTINUE
229*
230 norm = dlange( '1', n, n, t, ldt, work )
231 k = idamax( n, b, 1 )
232 normtb = norm + abs( b( k ) ) + abs( w )
233*
234 CALL dcopy( n, d, 1, x, 1 )
235 knt = knt + 1
236 CALL dlaqtr( .false., .true., n, t, ldt, dum,
237 $ dumm, scale, x, work, info )
238 IF( info.NE.0 )
239 $ ninfo = ninfo + 1
240*
241* || T*x - scale*d || /
242* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
243*
244 CALL dcopy( n, d, 1, y, 1 )
245 CALL dgemv( 'No transpose', n, n, one, t, ldt,
246 $ x, 1, -scale, y, 1 )
247 xnorm = dasum( n, x, 1 )
248 resid = dasum( n, y, 1 )
249 domin = max( smlnum, ( smlnum / eps )*norm,
250 $ ( norm*eps )*xnorm )
251 resid = resid / domin
252 IF( resid.GT.rmax ) THEN
253 rmax = resid
254 lmax = knt
255 END IF
256*
257 CALL dcopy( n, d, 1, x, 1 )
258 knt = knt + 1
259 CALL dlaqtr( .true., .true., n, t, ldt, dum,
260 $ dumm, scale, x, work, info )
261 IF( info.NE.0 )
262 $ ninfo = ninfo + 1
263*
264* || T*x - scale*d || /
265* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
266*
267 CALL dcopy( n, d, 1, y, 1 )
268 CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
269 $ 1, -scale, y, 1 )
270 xnorm = dasum( n, x, 1 )
271 resid = dasum( n, y, 1 )
272 domin = max( smlnum, ( smlnum / eps )*norm,
273 $ ( norm*eps )*xnorm )
274 resid = resid / domin
275 IF( resid.GT.rmax ) THEN
276 rmax = resid
277 lmax = knt
278 END IF
279*
280 CALL dcopy( 2*n, d, 1, x, 1 )
281 knt = knt + 1
282 CALL dlaqtr( .false., .false., n, t, ldt, b, w,
283 $ scale, x, work, info )
284 IF( info.NE.0 )
285 $ ninfo = ninfo + 1
286*
287* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
288* max(ulp*(||T||+||B||)*(||x1||+||x2||),
289* smlnum/ulp * (||T||+||B||), smlnum )
290*
291*
292 CALL dcopy( 2*n, d, 1, y, 1 )
293 y( 1 ) = ddot( n, b, 1, x( 1+n ), 1 ) +
294 $ scale*y( 1 )
295 DO 50 i = 2, n
296 y( i ) = w*x( i+n ) + scale*y( i )
297 50 CONTINUE
298 CALL dgemv( 'No transpose', n, n, one, t, ldt,
299 $ x, 1, -one, y, 1 )
300*
301 y( 1+n ) = ddot( n, b, 1, x, 1 ) -
302 $ scale*y( 1+n )
303 DO 60 i = 2, n
304 y( i+n ) = w*x( i ) - scale*y( i+n )
305 60 CONTINUE
306 CALL dgemv( 'No transpose', n, n, one, t, ldt,
307 $ x( 1+n ), 1, one, y( 1+n ), 1 )
308*
309 resid = dasum( 2*n, y, 1 )
310 domin = max( smlnum, ( smlnum / eps )*normtb,
311 $ eps*( normtb*dasum( 2*n, x, 1 ) ) )
312 resid = resid / domin
313 IF( resid.GT.rmax ) THEN
314 rmax = resid
315 lmax = knt
316 END IF
317*
318 CALL dcopy( 2*n, d, 1, x, 1 )
319 knt = knt + 1
320 CALL dlaqtr( .true., .false., n, t, ldt, b, w,
321 $ scale, x, work, info )
322 IF( info.NE.0 )
323 $ ninfo = ninfo + 1
324*
325* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
326* max(ulp*(||T||+||B||)*(||x1||+||x2||),
327* smlnum/ulp * (||T||+||B||), smlnum )
328*
329 CALL dcopy( 2*n, d, 1, y, 1 )
330 y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
331 DO 70 i = 2, n
332 y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
333 $ scale*y( i )
334 70 CONTINUE
335 CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
336 $ 1, one, y, 1 )
337*
338 y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
339 DO 80 i = 2, n
340 y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
341 $ scale*y( i+n )
342 80 CONTINUE
343 CALL dgemv( 'Transpose', n, n, one, t, ldt,
344 $ x( 1+n ), 1, -one, y( 1+n ), 1 )
345*
346 resid = dasum( 2*n, y, 1 )
347 domin = max( smlnum, ( smlnum / eps )*normtb,
348 $ eps*( normtb*dasum( 2*n, x, 1 ) ) )
349 resid = resid / domin
350 IF( resid.GT.rmax ) THEN
351 rmax = resid
352 lmax = knt
353 END IF
354*
355 90 CONTINUE
356 100 CONTINUE
357 110 CONTINUE
358 120 CONTINUE
359 130 CONTINUE
360 140 CONTINUE
361*
362 RETURN
363*
364* End of DGET39
365*
366 END
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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:156
dget39
subroutine dget39(RMAX, LMAX, NINFO, KNT)
DGET39
Definition: dget39.f:103
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