LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dget31.f
Go to the documentation of this file.
1*> \brief \b DGET31
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 DGET31( RMAX, LMAX, NINFO, KNT )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT, LMAX
15* DOUBLE PRECISION RMAX
16* ..
17* .. Array Arguments ..
18* INTEGER NINFO( 2 )
19* ..
20*
21*
22*> \par Purpose:
23* =============
24*>
25*> \verbatim
26*>
27*> DGET31 tests DLALN2, a routine for solving
28*>
29*> (ca A - w D)X = sB
30*>
31*> where A is an NA by NA matrix (NA=1 or 2 only), w is a real (NW=1) or
32*> complex (NW=2) constant, ca is a real constant, D is an NA by NA real
33*> diagonal matrix, and B is an NA by NW matrix (when NW=2 the second
34*> column of B contains the imaginary part of the solution). The code
35*> returns X and s, where s is a scale factor, less than or equal to 1,
36*> which is chosen to avoid overflow in X.
37*>
38*> If any singular values of ca A-w D are less than another input
39*> parameter SMIN, they are perturbed up to SMIN.
40*>
41*> The test condition is that the scaled residual
42*>
43*> norm( (ca A-w D)*X - s*B ) /
44*> ( max( ulp*norm(ca A-w D), SMIN )*norm(X) )
45*>
46*> should be on the order of 1. Here, ulp is the machine precision.
47*> Also, it is verified that SCALE is less than or equal to 1, and that
48*> XNORM = infinity-norm(X).
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[out] RMAX
55*> \verbatim
56*> RMAX is DOUBLE PRECISION
57*> Value of the largest test ratio.
58*> \endverbatim
59*>
60*> \param[out] LMAX
61*> \verbatim
62*> LMAX is INTEGER
63*> Example number where largest test ratio achieved.
64*> \endverbatim
65*>
66*> \param[out] NINFO
67*> \verbatim
68*> NINFO is INTEGER array, dimension (2)
69*> NINFO(1) = number of examples with INFO less than 0
70*> NINFO(2) = number of examples with INFO greater than 0
71*> \endverbatim
72*>
73*> \param[out] KNT
74*> \verbatim
75*> KNT is INTEGER
76*> Total number of examples tested.
77*> \endverbatim
78*
79* Authors:
80* ========
81*
82*> \author Univ. of Tennessee
83*> \author Univ. of California Berkeley
84*> \author Univ. of Colorado Denver
85*> \author NAG Ltd.
86*
87*> \ingroup double_eig
88*
89* =====================================================================
90 SUBROUTINE dget31( RMAX, LMAX, NINFO, KNT )
91*
92* -- LAPACK test routine --
93* -- LAPACK is a software package provided by Univ. of Tennessee, --
94* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95*
96* .. Scalar Arguments ..
97 INTEGER KNT, LMAX
98 DOUBLE PRECISION RMAX
99* ..
100* .. Array Arguments ..
101 INTEGER NINFO( 2 )
102* ..
103*
104* =====================================================================
105*
106* .. Parameters ..
107 DOUBLE PRECISION ZERO, HALF, ONE
108 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
109 DOUBLE PRECISION TWO, THREE, FOUR
110 parameter( two = 2.0d0, three = 3.0d0, four = 4.0d0 )
111 DOUBLE PRECISION SEVEN, TEN
112 parameter( seven = 7.0d0, ten = 10.0d0 )
113 DOUBLE PRECISION TWNONE
114 parameter( twnone = 21.0d0 )
115* ..
116* .. Local Scalars ..
117 INTEGER IA, IB, ICA, ID1, ID2, INFO, ISMIN, ITRANS,
118 $ IWI, IWR, NA, NW
119 DOUBLE PRECISION BIGNUM, CA, D1, D2, DEN, EPS, RES, SCALE, SMIN,
120 $ SMLNUM, TMP, UNFL, WI, WR, XNORM
121* ..
122* .. Local Arrays ..
123 LOGICAL LTRANS( 0: 1 )
124 DOUBLE PRECISION A( 2, 2 ), B( 2, 2 ), VAB( 3 ), VCA( 5 ),
125 $ VDD( 4 ), VSMIN( 4 ), VWI( 4 ), VWR( 4 ),
126 $ X( 2, 2 )
127* ..
128* .. External Functions ..
129 DOUBLE PRECISION DLAMCH
130 EXTERNAL dlamch
131* ..
132* .. External Subroutines ..
133 EXTERNAL dlaln2
134* ..
135* .. Intrinsic Functions ..
136 INTRINSIC abs, max, sqrt
137* ..
138* .. Data statements ..
139 DATA ltrans / .false., .true. /
140* ..
141* .. Executable Statements ..
142*
143* Get machine parameters
144*
145 eps = dlamch( 'P' )
146 unfl = dlamch( 'U' )
147 smlnum = dlamch( 'S' ) / eps
148 bignum = one / smlnum
149*
150* Set up test case parameters
151*
152 vsmin( 1 ) = smlnum
153 vsmin( 2 ) = eps
154 vsmin( 3 ) = one / ( ten*ten )
155 vsmin( 4 ) = one / eps
156 vab( 1 ) = sqrt( smlnum )
157 vab( 2 ) = one
158 vab( 3 ) = sqrt( bignum )
159 vwr( 1 ) = zero
160 vwr( 2 ) = half
161 vwr( 3 ) = two
162 vwr( 4 ) = one
163 vwi( 1 ) = smlnum
164 vwi( 2 ) = eps
165 vwi( 3 ) = one
166 vwi( 4 ) = two
167 vdd( 1 ) = sqrt( smlnum )
168 vdd( 2 ) = one
169 vdd( 3 ) = two
170 vdd( 4 ) = sqrt( bignum )
171 vca( 1 ) = zero
172 vca( 2 ) = sqrt( smlnum )
173 vca( 3 ) = eps
174 vca( 4 ) = half
175 vca( 5 ) = one
176*
177 knt = 0
178 ninfo( 1 ) = 0
179 ninfo( 2 ) = 0
180 lmax = 0
181 rmax = zero
182*
183* Begin test loop
184*
185 DO 190 id1 = 1, 4
186 d1 = vdd( id1 )
187 DO 180 id2 = 1, 4
188 d2 = vdd( id2 )
189 DO 170 ica = 1, 5
190 ca = vca( ica )
191 DO 160 itrans = 0, 1
192 DO 150 ismin = 1, 4
193 smin = vsmin( ismin )
194*
195 na = 1
196 nw = 1
197 DO 30 ia = 1, 3
198 a( 1, 1 ) = vab( ia )
199 DO 20 ib = 1, 3
200 b( 1, 1 ) = vab( ib )
201 DO 10 iwr = 1, 4
202 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
203 $ one ) THEN
204 wr = vwr( iwr )*a( 1, 1 )
205 ELSE
206 wr = vwr( iwr )
207 END IF
208 wi = zero
209 CALL dlaln2( ltrans( itrans ), na, nw,
210 $ smin, ca, a, 2, d1, d2, b, 2,
211 $ wr, wi, x, 2, scale, xnorm,
212 $ info )
213 IF( info.LT.0 )
214 $ ninfo( 1 ) = ninfo( 1 ) + 1
215 IF( info.GT.0 )
216 $ ninfo( 2 ) = ninfo( 2 ) + 1
217 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
218 $ x( 1, 1 )-scale*b( 1, 1 ) )
219 IF( info.EQ.0 ) THEN
220 den = max( eps*( abs( ( ca*a( 1,
221 $ 1 )-wr*d1 )*x( 1, 1 ) ) ),
222 $ smlnum )
223 ELSE
224 den = max( smin*abs( x( 1, 1 ) ),
225 $ smlnum )
226 END IF
227 res = res / den
228 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
229 $ abs( b( 1, 1 ) ).LE.smlnum*
230 $ abs( ca*a( 1, 1 )-wr*d1 ) )res = zero
231 IF( scale.GT.one )
232 $ res = res + one / eps
233 res = res + abs( xnorm-abs( x( 1, 1 ) ) )
234 $ / max( smlnum, xnorm ) / eps
235 IF( info.NE.0 .AND. info.NE.1 )
236 $ res = res + one / eps
237 knt = knt + 1
238 IF( res.GT.rmax ) THEN
239 lmax = knt
240 rmax = res
241 END IF
242 10 CONTINUE
243 20 CONTINUE
244 30 CONTINUE
245*
246 na = 1
247 nw = 2
248 DO 70 ia = 1, 3
249 a( 1, 1 ) = vab( ia )
250 DO 60 ib = 1, 3
251 b( 1, 1 ) = vab( ib )
252 b( 1, 2 ) = -half*vab( ib )
253 DO 50 iwr = 1, 4
254 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
255 $ one ) THEN
256 wr = vwr( iwr )*a( 1, 1 )
257 ELSE
258 wr = vwr( iwr )
259 END IF
260 DO 40 iwi = 1, 4
261 IF( d1.EQ.one .AND. d2.EQ.one .AND.
262 $ ca.EQ.one ) THEN
263 wi = vwi( iwi )*a( 1, 1 )
264 ELSE
265 wi = vwi( iwi )
266 END IF
267 CALL dlaln2( ltrans( itrans ), na, nw,
268 $ smin, ca, a, 2, d1, d2, b,
269 $ 2, wr, wi, x, 2, scale,
270 $ xnorm, info )
271 IF( info.LT.0 )
272 $ ninfo( 1 ) = ninfo( 1 ) + 1
273 IF( info.GT.0 )
274 $ ninfo( 2 ) = ninfo( 2 ) + 1
275 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
276 $ x( 1, 1 )+( wi*d1 )*x( 1, 2 )-
277 $ scale*b( 1, 1 ) )
278 res = res + abs( ( -wi*d1 )*x( 1, 1 )+
279 $ ( ca*a( 1, 1 )-wr*d1 )*x( 1, 2 )-
280 $ scale*b( 1, 2 ) )
281 IF( info.EQ.0 ) THEN
282 den = max( eps*( max( abs( ca*a( 1,
283 $ 1 )-wr*d1 ), abs( d1*wi ) )*
284 $ ( abs( x( 1, 1 ) )+abs( x( 1,
285 $ 2 ) ) ) ), smlnum )
286 ELSE
287 den = max( smin*( abs( x( 1,
288 $ 1 ) )+abs( x( 1, 2 ) ) ),
289 $ smlnum )
290 END IF
291 res = res / den
292 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
293 $ abs( x( 1, 2 ) ).LT.unfl .AND.
294 $ abs( b( 1, 1 ) ).LE.smlnum*
295 $ abs( ca*a( 1, 1 )-wr*d1 ) )
296 $ res = zero
297 IF( scale.GT.one )
298 $ res = res + one / eps
299 res = res + abs( xnorm-
300 $ abs( x( 1, 1 ) )-
301 $ abs( x( 1, 2 ) ) ) /
302 $ max( smlnum, xnorm ) / eps
303 IF( info.NE.0 .AND. info.NE.1 )
304 $ res = res + one / eps
305 knt = knt + 1
306 IF( res.GT.rmax ) THEN
307 lmax = knt
308 rmax = res
309 END IF
310 40 CONTINUE
311 50 CONTINUE
312 60 CONTINUE
313 70 CONTINUE
314*
315 na = 2
316 nw = 1
317 DO 100 ia = 1, 3
318 a( 1, 1 ) = vab( ia )
319 a( 1, 2 ) = -three*vab( ia )
320 a( 2, 1 ) = -seven*vab( ia )
321 a( 2, 2 ) = twnone*vab( ia )
322 DO 90 ib = 1, 3
323 b( 1, 1 ) = vab( ib )
324 b( 2, 1 ) = -two*vab( ib )
325 DO 80 iwr = 1, 4
326 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
327 $ one ) THEN
328 wr = vwr( iwr )*a( 1, 1 )
329 ELSE
330 wr = vwr( iwr )
331 END IF
332 wi = zero
333 CALL dlaln2( ltrans( itrans ), na, nw,
334 $ smin, ca, a, 2, d1, d2, b, 2,
335 $ wr, wi, x, 2, scale, xnorm,
336 $ info )
337 IF( info.LT.0 )
338 $ ninfo( 1 ) = ninfo( 1 ) + 1
339 IF( info.GT.0 )
340 $ ninfo( 2 ) = ninfo( 2 ) + 1
341 IF( itrans.EQ.1 ) THEN
342 tmp = a( 1, 2 )
343 a( 1, 2 ) = a( 2, 1 )
344 a( 2, 1 ) = tmp
345 END IF
346 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
347 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
348 $ x( 2, 1 )-scale*b( 1, 1 ) )
349 res = res + abs( ( ca*a( 2, 1 ) )*
350 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
351 $ x( 2, 1 )-scale*b( 2, 1 ) )
352 IF( info.EQ.0 ) THEN
353 den = max( eps*( max( abs( ca*a( 1,
354 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
355 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
356 $ 2 )-wr*d2 ) )*max( abs( x( 1,
357 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
358 $ smlnum )
359 ELSE
360 den = max( eps*( max( smin / eps,
361 $ max( abs( ca*a( 1,
362 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
363 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
364 $ 2 )-wr*d2 ) ) )*max( abs( x( 1,
365 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
366 $ smlnum )
367 END IF
368 res = res / den
369 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
370 $ abs( x( 2, 1 ) ).LT.unfl .AND.
371 $ abs( b( 1, 1 ) )+abs( b( 2, 1 ) ).LE.
372 $ smlnum*( abs( ca*a( 1,
373 $ 1 )-wr*d1 )+abs( ca*a( 1,
374 $ 2 ) )+abs( ca*a( 2,
375 $ 1 ) )+abs( ca*a( 2, 2 )-wr*d2 ) ) )
376 $ res = zero
377 IF( scale.GT.one )
378 $ res = res + one / eps
379 res = res + abs( xnorm-
380 $ max( abs( x( 1, 1 ) ), abs( x( 2,
381 $ 1 ) ) ) ) / max( smlnum, xnorm ) /
382 $ eps
383 IF( info.NE.0 .AND. info.NE.1 )
384 $ res = res + one / eps
385 knt = knt + 1
386 IF( res.GT.rmax ) THEN
387 lmax = knt
388 rmax = res
389 END IF
390 80 CONTINUE
391 90 CONTINUE
392 100 CONTINUE
393*
394 na = 2
395 nw = 2
396 DO 140 ia = 1, 3
397 a( 1, 1 ) = vab( ia )*two
398 a( 1, 2 ) = -three*vab( ia )
399 a( 2, 1 ) = -seven*vab( ia )
400 a( 2, 2 ) = twnone*vab( ia )
401 DO 130 ib = 1, 3
402 b( 1, 1 ) = vab( ib )
403 b( 2, 1 ) = -two*vab( ib )
404 b( 1, 2 ) = four*vab( ib )
405 b( 2, 2 ) = -seven*vab( ib )
406 DO 120 iwr = 1, 4
407 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
408 $ one ) THEN
409 wr = vwr( iwr )*a( 1, 1 )
410 ELSE
411 wr = vwr( iwr )
412 END IF
413 DO 110 iwi = 1, 4
414 IF( d1.EQ.one .AND. d2.EQ.one .AND.
415 $ ca.EQ.one ) THEN
416 wi = vwi( iwi )*a( 1, 1 )
417 ELSE
418 wi = vwi( iwi )
419 END IF
420 CALL dlaln2( ltrans( itrans ), na, nw,
421 $ smin, ca, a, 2, d1, d2, b,
422 $ 2, wr, wi, x, 2, scale,
423 $ xnorm, info )
424 IF( info.LT.0 )
425 $ ninfo( 1 ) = ninfo( 1 ) + 1
426 IF( info.GT.0 )
427 $ ninfo( 2 ) = ninfo( 2 ) + 1
428 IF( itrans.EQ.1 ) THEN
429 tmp = a( 1, 2 )
430 a( 1, 2 ) = a( 2, 1 )
431 a( 2, 1 ) = tmp
432 END IF
433 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
434 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
435 $ x( 2, 1 )+( wi*d1 )*x( 1, 2 )-
436 $ scale*b( 1, 1 ) )
437 res = res + abs( ( ca*a( 1,
438 $ 1 )-wr*d1 )*x( 1, 2 )+
439 $ ( ca*a( 1, 2 ) )*x( 2, 2 )-
440 $ ( wi*d1 )*x( 1, 1 )-scale*
441 $ b( 1, 2 ) )
442 res = res + abs( ( ca*a( 2, 1 ) )*
443 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
444 $ x( 2, 1 )+( wi*d2 )*x( 2, 2 )-
445 $ scale*b( 2, 1 ) )
446 res = res + abs( ( ca*a( 2, 1 ) )*
447 $ x( 1, 2 )+( ca*a( 2, 2 )-wr*d2 )*
448 $ x( 2, 2 )-( wi*d2 )*x( 2, 1 )-
449 $ scale*b( 2, 2 ) )
450 IF( info.EQ.0 ) THEN
451 den = max( eps*( max( abs( ca*a( 1,
452 $ 1 )-wr*d1 )+abs( ca*a( 1,
453 $ 2 ) )+abs( wi*d1 ),
454 $ abs( ca*a( 2,
455 $ 1 ) )+abs( ca*a( 2,
456 $ 2 )-wr*d2 )+abs( wi*d2 ) )*
457 $ max( abs( x( 1,
458 $ 1 ) )+abs( x( 2, 1 ) ),
459 $ abs( x( 1, 2 ) )+abs( x( 2,
460 $ 2 ) ) ) ), smlnum )
461 ELSE
462 den = max( eps*( max( smin / eps,
463 $ max( abs( ca*a( 1,
464 $ 1 )-wr*d1 )+abs( ca*a( 1,
465 $ 2 ) )+abs( wi*d1 ),
466 $ abs( ca*a( 2,
467 $ 1 ) )+abs( ca*a( 2,
468 $ 2 )-wr*d2 )+abs( wi*d2 ) ) )*
469 $ max( abs( x( 1,
470 $ 1 ) )+abs( x( 2, 1 ) ),
471 $ abs( x( 1, 2 ) )+abs( x( 2,
472 $ 2 ) ) ) ), smlnum )
473 END IF
474 res = res / den
475 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
476 $ abs( x( 2, 1 ) ).LT.unfl .AND.
477 $ abs( x( 1, 2 ) ).LT.unfl .AND.
478 $ abs( x( 2, 2 ) ).LT.unfl .AND.
479 $ abs( b( 1, 1 ) )+
480 $ abs( b( 2, 1 ) ).LE.smlnum*
481 $ ( abs( ca*a( 1, 1 )-wr*d1 )+
482 $ abs( ca*a( 1, 2 ) )+abs( ca*a( 2,
483 $ 1 ) )+abs( ca*a( 2,
484 $ 2 )-wr*d2 )+abs( wi*d2 )+abs( wi*
485 $ d1 ) ) )res = zero
486 IF( scale.GT.one )
487 $ res = res + one / eps
488 res = res + abs( xnorm-
489 $ max( abs( x( 1, 1 ) )+abs( x( 1,
490 $ 2 ) ), abs( x( 2,
491 $ 1 ) )+abs( x( 2, 2 ) ) ) ) /
492 $ max( smlnum, xnorm ) / eps
493 IF( info.NE.0 .AND. info.NE.1 )
494 $ res = res + one / eps
495 knt = knt + 1
496 IF( res.GT.rmax ) THEN
497 lmax = knt
498 rmax = res
499 END IF
500 110 CONTINUE
501 120 CONTINUE
502 130 CONTINUE
503 140 CONTINUE
504 150 CONTINUE
505 160 CONTINUE
506 170 CONTINUE
507 180 CONTINUE
508 190 CONTINUE
509*
510 RETURN
511*
512* End of DGET31
513*
514 END
dget31
subroutine dget31(rmax, lmax, ninfo, knt)
DGET31
Definition dget31.f:91
dlaln2
subroutine dlaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition dlaln2.f:218