OpenRadioss 2025.1.11
OpenRadioss project
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 (3)
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 dlabad, 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 CALL dlabad( smlnum, bignum )
150*
151* Set up test case parameters
152*
153 vsmin( 1 ) = smlnum
154 vsmin( 2 ) = eps
155 vsmin( 3 ) = one / ( ten*ten )
156 vsmin( 4 ) = one / eps
157 vab( 1 ) = sqrt( smlnum )
158 vab( 2 ) = one
159 vab( 3 ) = sqrt( bignum )
160 vwr( 1 ) = zero
161 vwr( 2 ) = half
162 vwr( 3 ) = two
163 vwr( 4 ) = one
164 vwi( 1 ) = smlnum
165 vwi( 2 ) = eps
166 vwi( 3 ) = one
167 vwi( 4 ) = two
168 vdd( 1 ) = sqrt( smlnum )
169 vdd( 2 ) = one
170 vdd( 3 ) = two
171 vdd( 4 ) = sqrt( bignum )
172 vca( 1 ) = zero
173 vca( 2 ) = sqrt( smlnum )
174 vca( 3 ) = eps
175 vca( 4 ) = half
176 vca( 5 ) = one
177*
178 knt = 0
179 ninfo( 1 ) = 0
180 ninfo( 2 ) = 0
181 lmax = 0
182 rmax = zero
183*
184* Begin test loop
185*
186 DO 190 id1 = 1, 4
187 d1 = vdd( id1 )
188 DO 180 id2 = 1, 4
189 d2 = vdd( id2 )
190 DO 170 ica = 1, 5
191 ca = vca( ica )
192 DO 160 itrans = 0, 1
193 DO 150 ismin = 1, 4
194 smin = vsmin( ismin )
195*
196 na = 1
197 nw = 1
198 DO 30 ia = 1, 3
199 a( 1, 1 ) = vab( ia )
200 DO 20 ib = 1, 3
201 b( 1, 1 ) = vab( ib )
202 DO 10 iwr = 1, 4
203 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
204 $ one ) THEN
205 wr = vwr( iwr )*a( 1, 1 )
206 ELSE
207 wr = vwr( iwr )
208 END IF
209 wi = zero
210 CALL dlaln2( ltrans( itrans ), na, nw,
211 $ smin, ca, a, 2, d1, d2, b, 2,
212 $ wr, wi, x, 2, scale, xnorm,
213 $ info )
214 IF( info.LT.0 )
215 $ ninfo( 1 ) = ninfo( 1 ) + 1
216 IF( info.GT.0 )
217 $ ninfo( 2 ) = ninfo( 2 ) + 1
218 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
219 $ x( 1, 1 )-scale*b( 1, 1 ) )
220 IF( info.EQ.0 ) THEN
221 den = max( eps*( abs( ( ca*a( 1,
222 $ 1 )-wr*d1 )*x( 1, 1 ) ) ),
223 $ smlnum )
224 ELSE
225 den = max( smin*abs( x( 1, 1 ) ),
226 $ smlnum )
227 END IF
228 res = res / den
229 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
230 $ abs( b( 1, 1 ) ).LE.smlnum*
231 $ abs( ca*a( 1, 1 )-wr*d1 ) )res = zero
232 IF( scale.GT.one )
233 $ res = res + one / eps
234 res = res + abs( xnorm-abs( x( 1, 1 ) ) )
235 $ / max( smlnum, xnorm ) / eps
236 IF( info.NE.0 .AND. info.NE.1 )
237 $ res = res + one / eps
238 knt = knt + 1
239 IF( res.GT.rmax ) THEN
240 lmax = knt
241 rmax = res
242 END IF
243 10 CONTINUE
244 20 CONTINUE
245 30 CONTINUE
246*
247 na = 1
248 nw = 2
249 DO 70 ia = 1, 3
250 a( 1, 1 ) = vab( ia )
251 DO 60 ib = 1, 3
252 b( 1, 1 ) = vab( ib )
253 b( 1, 2 ) = -half*vab( ib )
254 DO 50 iwr = 1, 4
255 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
256 $ one ) THEN
257 wr = vwr( iwr )*a( 1, 1 )
258 ELSE
259 wr = vwr( iwr )
260 END IF
261 DO 40 iwi = 1, 4
262 IF( d1.EQ.one .AND. d2.EQ.one .AND.
263 $ ca.EQ.one ) THEN
264 wi = vwi( iwi )*a( 1, 1 )
265 ELSE
266 wi = vwi( iwi )
267 END IF
268 CALL dlaln2( ltrans( itrans ), na, nw,
269 $ smin, ca, a, 2, d1, d2, b,
270 $ 2, wr, wi, x, 2, scale,
271 $ xnorm, info )
272 IF( info.LT.0 )
273 $ ninfo( 1 ) = ninfo( 1 ) + 1
274 IF( info.GT.0 )
275 $ ninfo( 2 ) = ninfo( 2 ) + 1
276 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
277 $ x( 1, 1 )+( wi*d1 )*x( 1, 2 )-
278 $ scale*b( 1, 1 ) )
279 res = res + abs( ( -wi*d1 )*x( 1, 1 )+
280 $ ( ca*a( 1, 1 )-wr*d1 )*x( 1, 2 )-
281 $ scale*b( 1, 2 ) )
282 IF( info.EQ.0 ) THEN
283 den = max( eps*( max( abs( ca*a( 1,
284 $ 1 )-wr*d1 ), abs( d1*wi ) )*
285 $ ( abs( x( 1, 1 ) )+abs( x( 1,
286 $ 2 ) ) ) ), smlnum )
287 ELSE
288 den = max( smin*( abs( x( 1,
289 $ 1 ) )+abs( x( 1, 2 ) ) ),
290 $ smlnum )
291 END IF
292 res = res / den
293 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
294 $ abs( x( 1, 2 ) ).LT.unfl .AND.
295 $ abs( b( 1, 1 ) ).LE.smlnum*
296 $ abs( ca*a( 1, 1 )-wr*d1 ) )
297 $ res = zero
298 IF( scale.GT.one )
299 $ res = res + one / eps
300 res = res + abs( xnorm-
301 $ abs( x( 1, 1 ) )-
302 $ abs( x( 1, 2 ) ) ) /
303 $ max( smlnum, xnorm ) / eps
304 IF( info.NE.0 .AND. info.NE.1 )
305 $ res = res + one / eps
306 knt = knt + 1
307 IF( res.GT.rmax ) THEN
308 lmax = knt
309 rmax = res
310 END IF
311 40 CONTINUE
312 50 CONTINUE
313 60 CONTINUE
314 70 CONTINUE
315*
316 na = 2
317 nw = 1
318 DO 100 ia = 1, 3
319 a( 1, 1 ) = vab( ia )
320 a( 1, 2 ) = -three*vab( ia )
321 a( 2, 1 ) = -seven*vab( ia )
322 a( 2, 2 ) = twnone*vab( ia )
323 DO 90 ib = 1, 3
324 b( 1, 1 ) = vab( ib )
325 b( 2, 1 ) = -two*vab( ib )
326 DO 80 iwr = 1, 4
327 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
328 $ one ) THEN
329 wr = vwr( iwr )*a( 1, 1 )
330 ELSE
331 wr = vwr( iwr )
332 END IF
333 wi = zero
334 CALL dlaln2( ltrans( itrans ), na, nw,
335 $ smin, ca, a, 2, d1, d2, b, 2,
336 $ wr, wi, x, 2, scale, xnorm,
337 $ info )
338 IF( info.LT.0 )
339 $ ninfo( 1 ) = ninfo( 1 ) + 1
340 IF( info.GT.0 )
341 $ ninfo( 2 ) = ninfo( 2 ) + 1
342 IF( itrans.EQ.1 ) THEN
343 tmp = a( 1, 2 )
344 a( 1, 2 ) = a( 2, 1 )
345 a( 2, 1 ) = tmp
346 END IF
347 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
348 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
349 $ x( 2, 1 )-scale*b( 1, 1 ) )
350 res = res + abs( ( ca*a( 2, 1 ) )*
351 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
352 $ x( 2, 1 )-scale*b( 2, 1 ) )
353 IF( info.EQ.0 ) THEN
354 den = max( eps*( max( abs( ca*a( 1,
355 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
356 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
357 $ 2 )-wr*d2 ) )*max( abs( x( 1,
358 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
359 $ smlnum )
360 ELSE
361 den = max( eps*( max( smin / eps,
362 $ max( abs( ca*a( 1,
363 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
364 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
365 $ 2 )-wr*d2 ) ) )*max( abs( x( 1,
366 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
367 $ smlnum )
368 END IF
369 res = res / den
370 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
371 $ abs( x( 2, 1 ) ).LT.unfl .AND.
372 $ abs( b( 1, 1 ) )+abs( b( 2, 1 ) ).LE.
373 $ smlnum*( abs( ca*a( 1,
374 $ 1 )-wr*d1 )+abs( ca*a( 1,
375 $ 2 ) )+abs( ca*a( 2,
376 $ 1 ) )+abs( ca*a( 2, 2 )-wr*d2 ) ) )
377 $ res = zero
378 IF( scale.GT.one )
379 $ res = res + one / eps
380 res = res + abs( xnorm-
381 $ max( abs( x( 1, 1 ) ), abs( x( 2,
382 $ 1 ) ) ) ) / max( smlnum, xnorm ) /
383 $ eps
384 IF( info.NE.0 .AND. info.NE.1 )
385 $ res = res + one / eps
386 knt = knt + 1
387 IF( res.GT.rmax ) THEN
388 lmax = knt
389 rmax = res
390 END IF
391 80 CONTINUE
392 90 CONTINUE
393 100 CONTINUE
394*
395 na = 2
396 nw = 2
397 DO 140 ia = 1, 3
398 a( 1, 1 ) = vab( ia )*two
399 a( 1, 2 ) = -three*vab( ia )
400 a( 2, 1 ) = -seven*vab( ia )
401 a( 2, 2 ) = twnone*vab( ia )
402 DO 130 ib = 1, 3
403 b( 1, 1 ) = vab( ib )
404 b( 2, 1 ) = -two*vab( ib )
405 b( 1, 2 ) = four*vab( ib )
406 b( 2, 2 ) = -seven*vab( ib )
407 DO 120 iwr = 1, 4
408 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
409 $ one ) THEN
410 wr = vwr( iwr )*a( 1, 1 )
411 ELSE
412 wr = vwr( iwr )
413 END IF
414 DO 110 iwi = 1, 4
415 IF( d1.EQ.one .AND. d2.EQ.one .AND.
416 $ ca.EQ.one ) THEN
417 wi = vwi( iwi )*a( 1, 1 )
418 ELSE
419 wi = vwi( iwi )
420 END IF
421 CALL dlaln2( ltrans( itrans ), na, nw,
422 $ smin, ca, a, 2, d1, d2, b,
423 $ 2, wr, wi, x, 2, scale,
424 $ xnorm, info )
425 IF( info.LT.0 )
426 $ ninfo( 1 ) = ninfo( 1 ) + 1
427 IF( info.GT.0 )
428 $ ninfo( 2 ) = ninfo( 2 ) + 1
429 IF( itrans.EQ.1 ) THEN
430 tmp = a( 1, 2 )
431 a( 1, 2 ) = a( 2, 1 )
432 a( 2, 1 ) = tmp
433 END IF
434 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
435 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
436 $ x( 2, 1 )+( wi*d1 )*x( 1, 2 )-
437 $ scale*b( 1, 1 ) )
438 res = res + abs( ( ca*a( 1,
439 $ 1 )-wr*d1 )*x( 1, 2 )+
440 $ ( ca*a( 1, 2 ) )*x( 2, 2 )-
441 $ ( wi*d1 )*x( 1, 1 )-scale*
442 $ b( 1, 2 ) )
443 res = res + abs( ( ca*a( 2, 1 ) )*
444 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
445 $ x( 2, 1 )+( wi*d2 )*x( 2, 2 )-
446 $ scale*b( 2, 1 ) )
447 res = res + abs( ( ca*a( 2, 1 ) )*
448 $ x( 1, 2 )+( ca*a( 2, 2 )-wr*d2 )*
449 $ x( 2, 2 )-( wi*d2 )*x( 2, 1 )-
450 $ scale*b( 2, 2 ) )
451 IF( info.EQ.0 ) THEN
452 den = max( eps*( max( abs( ca*a( 1,
453 $ 1 )-wr*d1 )+abs( ca*a( 1,
454 $ 2 ) )+abs( wi*d1 ),
455 $ abs( ca*a( 2,
456 $ 1 ) )+abs( ca*a( 2,
457 $ 2 )-wr*d2 )+abs( wi*d2 ) )*
458 $ max( abs( x( 1,
459 $ 1 ) )+abs( x( 2, 1 ) ),
460 $ abs( x( 1, 2 ) )+abs( x( 2,
461 $ 2 ) ) ) ), smlnum )
462 ELSE
463 den = max( eps*( max( smin / eps,
464 $ max( abs( ca*a( 1,
465 $ 1 )-wr*d1 )+abs( ca*a( 1,
466 $ 2 ) )+abs( wi*d1 ),
467 $ abs( ca*a( 2,
468 $ 1 ) )+abs( ca*a( 2,
469 $ 2 )-wr*d2 )+abs( wi*d2 ) ) )*
470 $ max( abs( x( 1,
471 $ 1 ) )+abs( x( 2, 1 ) ),
472 $ abs( x( 1, 2 ) )+abs( x( 2,
473 $ 2 ) ) ) ), smlnum )
474 END IF
475 res = res / den
476 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
477 $ abs( x( 2, 1 ) ).LT.unfl .AND.
478 $ abs( x( 1, 2 ) ).LT.unfl .AND.
479 $ abs( x( 2, 2 ) ).LT.unfl .AND.
480 $ abs( b( 1, 1 ) )+
481 $ abs( b( 2, 1 ) ).LE.smlnum*
482 $ ( abs( ca*a( 1, 1 )-wr*d1 )+
483 $ abs( ca*a( 1, 2 ) )+abs( ca*a( 2,
484 $ 1 ) )+abs( ca*a( 2,
485 $ 2 )-wr*d2 )+abs( wi*d2 )+abs( wi*
486 $ d1 ) ) )res = zero
487 IF( scale.GT.one )
488 $ res = res + one / eps
489 res = res + abs( xnorm-
490 $ max( abs( x( 1, 1 ) )+abs( x( 1,
491 $ 2 ) ), abs( x( 2,
492 $ 1 ) )+abs( x( 2, 2 ) ) ) ) /
493 $ max( smlnum, xnorm ) / eps
494 IF( info.NE.0 .AND. info.NE.1 )
495 $ res = res + one / eps
496 knt = knt + 1
497 IF( res.GT.rmax ) THEN
498 lmax = knt
499 rmax = res
500 END IF
501 110 CONTINUE
502 120 CONTINUE
503 130 CONTINUE
504 140 CONTINUE
505 150 CONTINUE
506 160 CONTINUE
507 170 CONTINUE
508 180 CONTINUE
509 190 CONTINUE
510*
511 RETURN
512*
513* End of DGET31
514*
515 END
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
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
dget31
subroutine dget31(rmax, lmax, ninfo, knt)
DGET31
Definition dget31.f:91
max
#define max(a, b)
Definition macros.h:21