LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 *> \date November 2011
88 *
89 *> \ingroup double_eig
90 *
91 * =====================================================================
92  SUBROUTINE dget31( RMAX, LMAX, NINFO, KNT )
93 *
94 * -- LAPACK test routine (version 3.4.0) --
95 * -- LAPACK is a software package provided by Univ. of Tennessee, --
96 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
97 * November 2011
98 *
99 * .. Scalar Arguments ..
100  INTEGER KNT, LMAX
101  DOUBLE PRECISION RMAX
102 * ..
103 * .. Array Arguments ..
104  INTEGER NINFO( 2 )
105 * ..
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110  DOUBLE PRECISION ZERO, HALF, ONE
111  parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
112  DOUBLE PRECISION TWO, THREE, FOUR
113  parameter ( two = 2.0d0, three = 3.0d0, four = 4.0d0 )
114  DOUBLE PRECISION SEVEN, TEN
115  parameter ( seven = 7.0d0, ten = 10.0d0 )
116  DOUBLE PRECISION TWNONE
117  parameter ( twnone = 21.0d0 )
118 * ..
119 * .. Local Scalars ..
120  INTEGER IA, IB, ICA, ID1, ID2, INFO, ISMIN, ITRANS,
121  $ iwi, iwr, na, nw
122  DOUBLE PRECISION BIGNUM, CA, D1, D2, DEN, EPS, RES, SCALE, SMIN,
123  $ smlnum, tmp, unfl, wi, wr, xnorm
124 * ..
125 * .. Local Arrays ..
126  LOGICAL LTRANS( 0: 1 )
127  DOUBLE PRECISION A( 2, 2 ), B( 2, 2 ), VAB( 3 ), VCA( 5 ),
128  $ vdd( 4 ), vsmin( 4 ), vwi( 4 ), vwr( 4 ),
129  $ x( 2, 2 )
130 * ..
131 * .. External Functions ..
132  DOUBLE PRECISION DLAMCH
133  EXTERNAL dlamch
134 * ..
135 * .. External Subroutines ..
136  EXTERNAL dlabad, dlaln2
137 * ..
138 * .. Intrinsic Functions ..
139  INTRINSIC abs, max, sqrt
140 * ..
141 * .. Data statements ..
142  DATA ltrans / .false., .true. /
143 * ..
144 * .. Executable Statements ..
145 *
146 * Get machine parameters
147 *
148  eps = dlamch( 'P' )
149  unfl = dlamch( 'U' )
150  smlnum = dlamch( 'S' ) / eps
151  bignum = one / smlnum
152  CALL dlabad( smlnum, bignum )
153 *
154 * Set up test case parameters
155 *
156  vsmin( 1 ) = smlnum
157  vsmin( 2 ) = eps
158  vsmin( 3 ) = one / ( ten*ten )
159  vsmin( 4 ) = one / eps
160  vab( 1 ) = sqrt( smlnum )
161  vab( 2 ) = one
162  vab( 3 ) = sqrt( bignum )
163  vwr( 1 ) = zero
164  vwr( 2 ) = half
165  vwr( 3 ) = two
166  vwr( 4 ) = one
167  vwi( 1 ) = smlnum
168  vwi( 2 ) = eps
169  vwi( 3 ) = one
170  vwi( 4 ) = two
171  vdd( 1 ) = sqrt( smlnum )
172  vdd( 2 ) = one
173  vdd( 3 ) = two
174  vdd( 4 ) = sqrt( bignum )
175  vca( 1 ) = zero
176  vca( 2 ) = sqrt( smlnum )
177  vca( 3 ) = eps
178  vca( 4 ) = half
179  vca( 5 ) = one
180 *
181  knt = 0
182  ninfo( 1 ) = 0
183  ninfo( 2 ) = 0
184  lmax = 0
185  rmax = zero
186 *
187 * Begin test loop
188 *
189  DO 190 id1 = 1, 4
190  d1 = vdd( id1 )
191  DO 180 id2 = 1, 4
192  d2 = vdd( id2 )
193  DO 170 ica = 1, 5
194  ca = vca( ica )
195  DO 160 itrans = 0, 1
196  DO 150 ismin = 1, 4
197  smin = vsmin( ismin )
198 *
199  na = 1
200  nw = 1
201  DO 30 ia = 1, 3
202  a( 1, 1 ) = vab( ia )
203  DO 20 ib = 1, 3
204  b( 1, 1 ) = vab( ib )
205  DO 10 iwr = 1, 4
206  IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
207  $ one ) THEN
208  wr = vwr( iwr )*a( 1, 1 )
209  ELSE
210  wr = vwr( iwr )
211  END IF
212  wi = zero
213  CALL dlaln2( ltrans( itrans ), na, nw,
214  $ smin, ca, a, 2, d1, d2, b, 2,
215  $ wr, wi, x, 2, scale, xnorm,
216  $ info )
217  IF( info.LT.0 )
218  $ ninfo( 1 ) = ninfo( 1 ) + 1
219  IF( info.GT.0 )
220  $ ninfo( 2 ) = ninfo( 2 ) + 1
221  res = abs( ( ca*a( 1, 1 )-wr*d1 )*
222  $ x( 1, 1 )-scale*b( 1, 1 ) )
223  IF( info.EQ.0 ) THEN
224  den = max( eps*( abs( ( ca*a( 1,
225  $ 1 )-wr*d1 )*x( 1, 1 ) ) ),
226  $ smlnum )
227  ELSE
228  den = max( smin*abs( x( 1, 1 ) ),
229  $ smlnum )
230  END IF
231  res = res / den
232  IF( abs( x( 1, 1 ) ).LT.unfl .AND.
233  $ abs( b( 1, 1 ) ).LE.smlnum*
234  $ abs( ca*a( 1, 1 )-wr*d1 ) )res = zero
235  IF( scale.GT.one )
236  $ res = res + one / eps
237  res = res + abs( xnorm-abs( x( 1, 1 ) ) )
238  $ / max( smlnum, xnorm ) / eps
239  IF( info.NE.0 .AND. info.NE.1 )
240  $ res = res + one / eps
241  knt = knt + 1
242  IF( res.GT.rmax ) THEN
243  lmax = knt
244  rmax = res
245  END IF
246  10 CONTINUE
247  20 CONTINUE
248  30 CONTINUE
249 *
250  na = 1
251  nw = 2
252  DO 70 ia = 1, 3
253  a( 1, 1 ) = vab( ia )
254  DO 60 ib = 1, 3
255  b( 1, 1 ) = vab( ib )
256  b( 1, 2 ) = -half*vab( ib )
257  DO 50 iwr = 1, 4
258  IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
259  $ one ) THEN
260  wr = vwr( iwr )*a( 1, 1 )
261  ELSE
262  wr = vwr( iwr )
263  END IF
264  DO 40 iwi = 1, 4
265  IF( d1.EQ.one .AND. d2.EQ.one .AND.
266  $ ca.EQ.one ) THEN
267  wi = vwi( iwi )*a( 1, 1 )
268  ELSE
269  wi = vwi( iwi )
270  END IF
271  CALL dlaln2( ltrans( itrans ), na, nw,
272  $ smin, ca, a, 2, d1, d2, b,
273  $ 2, wr, wi, x, 2, scale,
274  $ xnorm, info )
275  IF( info.LT.0 )
276  $ ninfo( 1 ) = ninfo( 1 ) + 1
277  IF( info.GT.0 )
278  $ ninfo( 2 ) = ninfo( 2 ) + 1
279  res = abs( ( ca*a( 1, 1 )-wr*d1 )*
280  $ x( 1, 1 )+( wi*d1 )*x( 1, 2 )-
281  $ scale*b( 1, 1 ) )
282  res = res + abs( ( -wi*d1 )*x( 1, 1 )+
283  $ ( ca*a( 1, 1 )-wr*d1 )*x( 1, 2 )-
284  $ scale*b( 1, 2 ) )
285  IF( info.EQ.0 ) THEN
286  den = max( eps*( max( abs( ca*a( 1,
287  $ 1 )-wr*d1 ), abs( d1*wi ) )*
288  $ ( abs( x( 1, 1 ) )+abs( x( 1,
289  $ 2 ) ) ) ), smlnum )
290  ELSE
291  den = max( smin*( abs( x( 1,
292  $ 1 ) )+abs( x( 1, 2 ) ) ),
293  $ smlnum )
294  END IF
295  res = res / den
296  IF( abs( x( 1, 1 ) ).LT.unfl .AND.
297  $ abs( x( 1, 2 ) ).LT.unfl .AND.
298  $ abs( b( 1, 1 ) ).LE.smlnum*
299  $ abs( ca*a( 1, 1 )-wr*d1 ) )
300  $ res = zero
301  IF( scale.GT.one )
302  $ res = res + one / eps
303  res = res + abs( xnorm-
304  $ abs( x( 1, 1 ) )-
305  $ abs( x( 1, 2 ) ) ) /
306  $ max( smlnum, xnorm ) / eps
307  IF( info.NE.0 .AND. info.NE.1 )
308  $ res = res + one / eps
309  knt = knt + 1
310  IF( res.GT.rmax ) THEN
311  lmax = knt
312  rmax = res
313  END IF
314  40 CONTINUE
315  50 CONTINUE
316  60 CONTINUE
317  70 CONTINUE
318 *
319  na = 2
320  nw = 1
321  DO 100 ia = 1, 3
322  a( 1, 1 ) = vab( ia )
323  a( 1, 2 ) = -three*vab( ia )
324  a( 2, 1 ) = -seven*vab( ia )
325  a( 2, 2 ) = twnone*vab( ia )
326  DO 90 ib = 1, 3
327  b( 1, 1 ) = vab( ib )
328  b( 2, 1 ) = -two*vab( ib )
329  DO 80 iwr = 1, 4
330  IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
331  $ one ) THEN
332  wr = vwr( iwr )*a( 1, 1 )
333  ELSE
334  wr = vwr( iwr )
335  END IF
336  wi = zero
337  CALL dlaln2( ltrans( itrans ), na, nw,
338  $ smin, ca, a, 2, d1, d2, b, 2,
339  $ wr, wi, x, 2, scale, xnorm,
340  $ info )
341  IF( info.LT.0 )
342  $ ninfo( 1 ) = ninfo( 1 ) + 1
343  IF( info.GT.0 )
344  $ ninfo( 2 ) = ninfo( 2 ) + 1
345  IF( itrans.EQ.1 ) THEN
346  tmp = a( 1, 2 )
347  a( 1, 2 ) = a( 2, 1 )
348  a( 2, 1 ) = tmp
349  END IF
350  res = abs( ( ca*a( 1, 1 )-wr*d1 )*
351  $ x( 1, 1 )+( ca*a( 1, 2 ) )*
352  $ x( 2, 1 )-scale*b( 1, 1 ) )
353  res = res + abs( ( ca*a( 2, 1 ) )*
354  $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
355  $ x( 2, 1 )-scale*b( 2, 1 ) )
356  IF( info.EQ.0 ) THEN
357  den = max( eps*( max( abs( ca*a( 1,
358  $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
359  $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
360  $ 2 )-wr*d2 ) )*max( abs( x( 1,
361  $ 1 ) ), abs( x( 2, 1 ) ) ) ),
362  $ smlnum )
363  ELSE
364  den = max( eps*( max( smin / eps,
365  $ max( abs( ca*a( 1,
366  $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
367  $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
368  $ 2 )-wr*d2 ) ) )*max( abs( x( 1,
369  $ 1 ) ), abs( x( 2, 1 ) ) ) ),
370  $ smlnum )
371  END IF
372  res = res / den
373  IF( abs( x( 1, 1 ) ).LT.unfl .AND.
374  $ abs( x( 2, 1 ) ).LT.unfl .AND.
375  $ abs( b( 1, 1 ) )+abs( b( 2, 1 ) ).LE.
376  $ smlnum*( abs( ca*a( 1,
377  $ 1 )-wr*d1 )+abs( ca*a( 1,
378  $ 2 ) )+abs( ca*a( 2,
379  $ 1 ) )+abs( ca*a( 2, 2 )-wr*d2 ) ) )
380  $ res = zero
381  IF( scale.GT.one )
382  $ res = res + one / eps
383  res = res + abs( xnorm-
384  $ max( abs( x( 1, 1 ) ), abs( x( 2,
385  $ 1 ) ) ) ) / max( smlnum, xnorm ) /
386  $ eps
387  IF( info.NE.0 .AND. info.NE.1 )
388  $ res = res + one / eps
389  knt = knt + 1
390  IF( res.GT.rmax ) THEN
391  lmax = knt
392  rmax = res
393  END IF
394  80 CONTINUE
395  90 CONTINUE
396  100 CONTINUE
397 *
398  na = 2
399  nw = 2
400  DO 140 ia = 1, 3
401  a( 1, 1 ) = vab( ia )*two
402  a( 1, 2 ) = -three*vab( ia )
403  a( 2, 1 ) = -seven*vab( ia )
404  a( 2, 2 ) = twnone*vab( ia )
405  DO 130 ib = 1, 3
406  b( 1, 1 ) = vab( ib )
407  b( 2, 1 ) = -two*vab( ib )
408  b( 1, 2 ) = four*vab( ib )
409  b( 2, 2 ) = -seven*vab( ib )
410  DO 120 iwr = 1, 4
411  IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
412  $ one ) THEN
413  wr = vwr( iwr )*a( 1, 1 )
414  ELSE
415  wr = vwr( iwr )
416  END IF
417  DO 110 iwi = 1, 4
418  IF( d1.EQ.one .AND. d2.EQ.one .AND.
419  $ ca.EQ.one ) THEN
420  wi = vwi( iwi )*a( 1, 1 )
421  ELSE
422  wi = vwi( iwi )
423  END IF
424  CALL dlaln2( ltrans( itrans ), na, nw,
425  $ smin, ca, a, 2, d1, d2, b,
426  $ 2, wr, wi, x, 2, scale,
427  $ xnorm, info )
428  IF( info.LT.0 )
429  $ ninfo( 1 ) = ninfo( 1 ) + 1
430  IF( info.GT.0 )
431  $ ninfo( 2 ) = ninfo( 2 ) + 1
432  IF( itrans.EQ.1 ) THEN
433  tmp = a( 1, 2 )
434  a( 1, 2 ) = a( 2, 1 )
435  a( 2, 1 ) = tmp
436  END IF
437  res = abs( ( ca*a( 1, 1 )-wr*d1 )*
438  $ x( 1, 1 )+( ca*a( 1, 2 ) )*
439  $ x( 2, 1 )+( wi*d1 )*x( 1, 2 )-
440  $ scale*b( 1, 1 ) )
441  res = res + abs( ( ca*a( 1,
442  $ 1 )-wr*d1 )*x( 1, 2 )+
443  $ ( ca*a( 1, 2 ) )*x( 2, 2 )-
444  $ ( wi*d1 )*x( 1, 1 )-scale*
445  $ b( 1, 2 ) )
446  res = res + abs( ( ca*a( 2, 1 ) )*
447  $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
448  $ x( 2, 1 )+( wi*d2 )*x( 2, 2 )-
449  $ scale*b( 2, 1 ) )
450  res = res + abs( ( ca*a( 2, 1 ) )*
451  $ x( 1, 2 )+( ca*a( 2, 2 )-wr*d2 )*
452  $ x( 2, 2 )-( wi*d2 )*x( 2, 1 )-
453  $ scale*b( 2, 2 ) )
454  IF( info.EQ.0 ) THEN
455  den = max( eps*( max( abs( ca*a( 1,
456  $ 1 )-wr*d1 )+abs( ca*a( 1,
457  $ 2 ) )+abs( wi*d1 ),
458  $ abs( ca*a( 2,
459  $ 1 ) )+abs( ca*a( 2,
460  $ 2 )-wr*d2 )+abs( wi*d2 ) )*
461  $ max( abs( x( 1,
462  $ 1 ) )+abs( x( 2, 1 ) ),
463  $ abs( x( 1, 2 ) )+abs( x( 2,
464  $ 2 ) ) ) ), smlnum )
465  ELSE
466  den = max( eps*( max( smin / eps,
467  $ max( abs( ca*a( 1,
468  $ 1 )-wr*d1 )+abs( ca*a( 1,
469  $ 2 ) )+abs( wi*d1 ),
470  $ abs( ca*a( 2,
471  $ 1 ) )+abs( ca*a( 2,
472  $ 2 )-wr*d2 )+abs( wi*d2 ) ) )*
473  $ max( abs( x( 1,
474  $ 1 ) )+abs( x( 2, 1 ) ),
475  $ abs( x( 1, 2 ) )+abs( x( 2,
476  $ 2 ) ) ) ), smlnum )
477  END IF
478  res = res / den
479  IF( abs( x( 1, 1 ) ).LT.unfl .AND.
480  $ abs( x( 2, 1 ) ).LT.unfl .AND.
481  $ abs( x( 1, 2 ) ).LT.unfl .AND.
482  $ abs( x( 2, 2 ) ).LT.unfl .AND.
483  $ abs( b( 1, 1 ) )+
484  $ abs( b( 2, 1 ) ).LE.smlnum*
485  $ ( abs( ca*a( 1, 1 )-wr*d1 )+
486  $ abs( ca*a( 1, 2 ) )+abs( ca*a( 2,
487  $ 1 ) )+abs( ca*a( 2,
488  $ 2 )-wr*d2 )+abs( wi*d2 )+abs( wi*
489  $ d1 ) ) )res = zero
490  IF( scale.GT.one )
491  $ res = res + one / eps
492  res = res + abs( xnorm-
493  $ max( abs( x( 1, 1 ) )+abs( x( 1,
494  $ 2 ) ), abs( x( 2,
495  $ 1 ) )+abs( x( 2, 2 ) ) ) ) /
496  $ max( smlnum, xnorm ) / eps
497  IF( info.NE.0 .AND. info.NE.1 )
498  $ res = res + one / eps
499  knt = knt + 1
500  IF( res.GT.rmax ) THEN
501  lmax = knt
502  rmax = res
503  END IF
504  110 CONTINUE
505  120 CONTINUE
506  130 CONTINUE
507  140 CONTINUE
508  150 CONTINUE
509  160 CONTINUE
510  170 CONTINUE
511  180 CONTINUE
512  190 CONTINUE
513 *
514  RETURN
515 *
516 * End of DGET31
517 *
518  END
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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:220
dget31
subroutine dget31(RMAX, LMAX, NINFO, KNT)
DGET31
Definition: dget31.f:93