LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dget37.f
Go to the documentation of this file.
1 *> \brief \b DGET37
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 DGET37( RMAX, LMAX, NINFO, KNT, NIN )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER KNT, NIN
15 * ..
16 * .. Array Arguments ..
17 * INTEGER LMAX( 3 ), NINFO( 3 )
18 * DOUBLE PRECISION RMAX( 3 )
19 * ..
20 *
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> DGET37 tests DTRSNA, a routine for estimating condition numbers of
28 *> eigenvalues and/or right eigenvectors of a matrix.
29 *>
30 *> The test matrices are read from a file with logical unit number NIN.
31 *> \endverbatim
32 *
33 * Arguments:
34 * ==========
35 *
36 *> \param[out] RMAX
37 *> \verbatim
38 *> RMAX is DOUBLE PRECISION array, dimension (3)
39 *> Value of the largest test ratio.
40 *> RMAX(1) = largest ratio comparing different calls to DTRSNA
41 *> RMAX(2) = largest error in reciprocal condition
42 *> numbers taking their conditioning into account
43 *> RMAX(3) = largest error in reciprocal condition
44 *> numbers not taking their conditioning into
45 *> account (may be larger than RMAX(2))
46 *> \endverbatim
47 *>
48 *> \param[out] LMAX
49 *> \verbatim
50 *> LMAX is INTEGER array, dimension (3)
51 *> LMAX(i) is example number where largest test ratio
52 *> RMAX(i) is achieved. Also:
53 *> If DGEHRD returns INFO nonzero on example i, LMAX(1)=i
54 *> If DHSEQR returns INFO nonzero on example i, LMAX(2)=i
55 *> If DTRSNA returns INFO nonzero on example i, LMAX(3)=i
56 *> \endverbatim
57 *>
58 *> \param[out] NINFO
59 *> \verbatim
60 *> NINFO is INTEGER array, dimension (3)
61 *> NINFO(1) = No. of times DGEHRD returned INFO nonzero
62 *> NINFO(2) = No. of times DHSEQR returned INFO nonzero
63 *> NINFO(3) = No. of times DTRSNA returned INFO nonzero
64 *> \endverbatim
65 *>
66 *> \param[out] KNT
67 *> \verbatim
68 *> KNT is INTEGER
69 *> Total number of examples tested.
70 *> \endverbatim
71 *>
72 *> \param[in] NIN
73 *> \verbatim
74 *> NIN is INTEGER
75 *> Input logical unit number
76 *> \endverbatim
77 *
78 * Authors:
79 * ========
80 *
81 *> \author Univ. of Tennessee
82 *> \author Univ. of California Berkeley
83 *> \author Univ. of Colorado Denver
84 *> \author NAG Ltd.
85 *
86 *> \date November 2011
87 *
88 *> \ingroup double_eig
89 *
90 * =====================================================================
91  SUBROUTINE dget37( RMAX, LMAX, NINFO, KNT, NIN )
92 *
93 * -- LAPACK test routine (version 3.4.0) --
94 * -- LAPACK is a software package provided by Univ. of Tennessee, --
95 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96 * November 2011
97 *
98 * .. Scalar Arguments ..
99  INTEGER knt, nin
100 * ..
101 * .. Array Arguments ..
102  INTEGER lmax( 3 ), ninfo( 3 )
103  DOUBLE PRECISION rmax( 3 )
104 * ..
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109  DOUBLE PRECISION zero, one, two
110  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
111  DOUBLE PRECISION epsin
112  parameter( epsin = 5.9605d-8 )
113  INTEGER ldt, lwork
114  parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
115 * ..
116 * .. Local Scalars ..
117  INTEGER i, icmp, ifnd, info, iscl, j, kmin, m, n
118  DOUBLE PRECISION bignum, eps, smlnum, tnrm, tol, tolin, v,
119  $ vimin, vmax, vmul, vrmin
120 * ..
121 * .. Local Arrays ..
122  LOGICAL select( ldt )
123  INTEGER iwork( 2*ldt ), lcmp( 3 )
124  DOUBLE PRECISION dum( 1 ), le( ldt, ldt ), re( ldt, ldt ),
125  $ s( ldt ), sep( ldt ), sepin( ldt ),
126  $ septmp( ldt ), sin( ldt ), stmp( ldt ),
127  $ t( ldt, ldt ), tmp( ldt, ldt ), val( 3 ),
128  $ wi( ldt ), wiin( ldt ), witmp( ldt ),
129  $ work( lwork ), wr( ldt ), wrin( ldt ),
130  $ wrtmp( ldt )
131 * ..
132 * .. External Functions ..
133  DOUBLE PRECISION dlamch, dlange
134  EXTERNAL dlamch, dlange
135 * ..
136 * .. External Subroutines ..
137  EXTERNAL dcopy, dgehrd, dhseqr, dlabad, dlacpy, dscal,
138  $ dtrevc, dtrsna
139 * ..
140 * .. Intrinsic Functions ..
141  INTRINSIC dble, max, sqrt
142 * ..
143 * .. Executable Statements ..
144 *
145  eps = dlamch( 'P' )
146  smlnum = dlamch( 'S' ) / eps
147  bignum = one / smlnum
148  CALL dlabad( smlnum, bignum )
149 *
150 * EPSIN = 2**(-24) = precision to which input data computed
151 *
152  eps = max( eps, epsin )
153  rmax( 1 ) = zero
154  rmax( 2 ) = zero
155  rmax( 3 ) = zero
156  lmax( 1 ) = 0
157  lmax( 2 ) = 0
158  lmax( 3 ) = 0
159  knt = 0
160  ninfo( 1 ) = 0
161  ninfo( 2 ) = 0
162  ninfo( 3 ) = 0
163 *
164  val( 1 ) = sqrt( smlnum )
165  val( 2 ) = one
166  val( 3 ) = sqrt( bignum )
167 *
168 * Read input data until N=0. Assume input eigenvalues are sorted
169 * lexicographically (increasing by real part, then decreasing by
170 * imaginary part)
171 *
172  10 continue
173  READ( nin, fmt = * )n
174  IF( n.EQ.0 )
175  $ return
176  DO 20 i = 1, n
177  READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
178  20 continue
179  DO 30 i = 1, n
180  READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
181  30 continue
182  tnrm = dlange( 'M', n, n, tmp, ldt, work )
183 *
184 * Begin test
185 *
186  DO 240 iscl = 1, 3
187 *
188 * Scale input matrix
189 *
190  knt = knt + 1
191  CALL dlacpy( 'F', n, n, tmp, ldt, t, ldt )
192  vmul = val( iscl )
193  DO 40 i = 1, n
194  CALL dscal( n, vmul, t( 1, i ), 1 )
195  40 continue
196  IF( tnrm.EQ.zero )
197  $ vmul = one
198 *
199 * Compute eigenvalues and eigenvectors
200 *
201  CALL dgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
202  $ info )
203  IF( info.NE.0 ) THEN
204  lmax( 1 ) = knt
205  ninfo( 1 ) = ninfo( 1 ) + 1
206  go to 240
207  END IF
208  DO 60 j = 1, n - 2
209  DO 50 i = j + 2, n
210  t( i, j ) = zero
211  50 continue
212  60 continue
213 *
214 * Compute Schur form
215 *
216  CALL dhseqr( 'S', 'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
217  $ lwork, info )
218  IF( info.NE.0 ) THEN
219  lmax( 2 ) = knt
220  ninfo( 2 ) = ninfo( 2 ) + 1
221  go to 240
222  END IF
223 *
224 * Compute eigenvectors
225 *
226  CALL dtrevc( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
227  $ ldt, n, m, work, info )
228 *
229 * Compute condition numbers
230 *
231  CALL dtrsna( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
232  $ ldt, s, sep, n, m, work, n, iwork, info )
233  IF( info.NE.0 ) THEN
234  lmax( 3 ) = knt
235  ninfo( 3 ) = ninfo( 3 ) + 1
236  go to 240
237  END IF
238 *
239 * Sort eigenvalues and condition numbers lexicographically
240 * to compare with inputs
241 *
242  CALL dcopy( n, wr, 1, wrtmp, 1 )
243  CALL dcopy( n, wi, 1, witmp, 1 )
244  CALL dcopy( n, s, 1, stmp, 1 )
245  CALL dcopy( n, sep, 1, septmp, 1 )
246  CALL dscal( n, one / vmul, septmp, 1 )
247  DO 80 i = 1, n - 1
248  kmin = i
249  vrmin = wrtmp( i )
250  vimin = witmp( i )
251  DO 70 j = i + 1, n
252  IF( wrtmp( j ).LT.vrmin ) THEN
253  kmin = j
254  vrmin = wrtmp( j )
255  vimin = witmp( j )
256  END IF
257  70 continue
258  wrtmp( kmin ) = wrtmp( i )
259  witmp( kmin ) = witmp( i )
260  wrtmp( i ) = vrmin
261  witmp( i ) = vimin
262  vrmin = stmp( kmin )
263  stmp( kmin ) = stmp( i )
264  stmp( i ) = vrmin
265  vrmin = septmp( kmin )
266  septmp( kmin ) = septmp( i )
267  septmp( i ) = vrmin
268  80 continue
269 *
270 * Compare condition numbers for eigenvalues
271 * taking their condition numbers into account
272 *
273  v = max( two*dble( n )*eps*tnrm, smlnum )
274  IF( tnrm.EQ.zero )
275  $ v = one
276  DO 90 i = 1, n
277  IF( v.GT.septmp( i ) ) THEN
278  tol = one
279  ELSE
280  tol = v / septmp( i )
281  END IF
282  IF( v.GT.sepin( i ) ) THEN
283  tolin = one
284  ELSE
285  tolin = v / sepin( i )
286  END IF
287  tol = max( tol, smlnum / eps )
288  tolin = max( tolin, smlnum / eps )
289  IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
290  vmax = one / eps
291  ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
292  vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
293  ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
294  vmax = one / eps
295  ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
296  vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
297  ELSE
298  vmax = one
299  END IF
300  IF( vmax.GT.rmax( 2 ) ) THEN
301  rmax( 2 ) = vmax
302  IF( ninfo( 2 ).EQ.0 )
303  $ lmax( 2 ) = knt
304  END IF
305  90 continue
306 *
307 * Compare condition numbers for eigenvectors
308 * taking their condition numbers into account
309 *
310  DO 100 i = 1, n
311  IF( v.GT.septmp( i )*stmp( i ) ) THEN
312  tol = septmp( i )
313  ELSE
314  tol = v / stmp( i )
315  END IF
316  IF( v.GT.sepin( i )*sin( i ) ) THEN
317  tolin = sepin( i )
318  ELSE
319  tolin = v / sin( i )
320  END IF
321  tol = max( tol, smlnum / eps )
322  tolin = max( tolin, smlnum / eps )
323  IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
324  vmax = one / eps
325  ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
326  vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
327  ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
328  vmax = one / eps
329  ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
330  vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
331  ELSE
332  vmax = one
333  END IF
334  IF( vmax.GT.rmax( 2 ) ) THEN
335  rmax( 2 ) = vmax
336  IF( ninfo( 2 ).EQ.0 )
337  $ lmax( 2 ) = knt
338  END IF
339  100 continue
340 *
341 * Compare condition numbers for eigenvalues
342 * without taking their condition numbers into account
343 *
344  DO 110 i = 1, n
345  IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
346  $ dble( 2*n )*eps ) THEN
347  vmax = one
348  ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
349  vmax = one / eps
350  ELSE IF( sin( i ).GT.stmp( i ) ) THEN
351  vmax = sin( i ) / stmp( i )
352  ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
353  vmax = one / eps
354  ELSE IF( sin( i ).LT.stmp( i ) ) THEN
355  vmax = stmp( i ) / sin( i )
356  ELSE
357  vmax = one
358  END IF
359  IF( vmax.GT.rmax( 3 ) ) THEN
360  rmax( 3 ) = vmax
361  IF( ninfo( 3 ).EQ.0 )
362  $ lmax( 3 ) = knt
363  END IF
364  110 continue
365 *
366 * Compare condition numbers for eigenvectors
367 * without taking their condition numbers into account
368 *
369  DO 120 i = 1, n
370  IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
371  vmax = one
372  ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
373  vmax = one / eps
374  ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
375  vmax = sepin( i ) / septmp( i )
376  ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
377  vmax = one / eps
378  ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
379  vmax = septmp( i ) / sepin( i )
380  ELSE
381  vmax = one
382  END IF
383  IF( vmax.GT.rmax( 3 ) ) THEN
384  rmax( 3 ) = vmax
385  IF( ninfo( 3 ).EQ.0 )
386  $ lmax( 3 ) = knt
387  END IF
388  120 continue
389 *
390 * Compute eigenvalue condition numbers only and compare
391 *
392  vmax = zero
393  dum( 1 ) = -one
394  CALL dcopy( n, dum, 0, stmp, 1 )
395  CALL dcopy( n, dum, 0, septmp, 1 )
396  CALL dtrsna( 'Eigcond', 'All', SELECT, n, t, ldt, le, ldt, re,
397  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
398  IF( info.NE.0 ) THEN
399  lmax( 3 ) = knt
400  ninfo( 3 ) = ninfo( 3 ) + 1
401  go to 240
402  END IF
403  DO 130 i = 1, n
404  IF( stmp( i ).NE.s( i ) )
405  $ vmax = one / eps
406  IF( septmp( i ).NE.dum( 1 ) )
407  $ vmax = one / eps
408  130 continue
409 *
410 * Compute eigenvector condition numbers only and compare
411 *
412  CALL dcopy( n, dum, 0, stmp, 1 )
413  CALL dcopy( n, dum, 0, septmp, 1 )
414  CALL dtrsna( 'Veccond', 'All', SELECT, n, t, ldt, le, ldt, re,
415  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
416  IF( info.NE.0 ) THEN
417  lmax( 3 ) = knt
418  ninfo( 3 ) = ninfo( 3 ) + 1
419  go to 240
420  END IF
421  DO 140 i = 1, n
422  IF( stmp( i ).NE.dum( 1 ) )
423  $ vmax = one / eps
424  IF( septmp( i ).NE.sep( i ) )
425  $ vmax = one / eps
426  140 continue
427 *
428 * Compute all condition numbers using SELECT and compare
429 *
430  DO 150 i = 1, n
431  SELECT( i ) = .true.
432  150 continue
433  CALL dcopy( n, dum, 0, stmp, 1 )
434  CALL dcopy( n, dum, 0, septmp, 1 )
435  CALL dtrsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
436  $ re, ldt, stmp, septmp, n, m, work, n, iwork,
437  $ info )
438  IF( info.NE.0 ) THEN
439  lmax( 3 ) = knt
440  ninfo( 3 ) = ninfo( 3 ) + 1
441  go to 240
442  END IF
443  DO 160 i = 1, n
444  IF( septmp( i ).NE.sep( i ) )
445  $ vmax = one / eps
446  IF( stmp( i ).NE.s( i ) )
447  $ vmax = one / eps
448  160 continue
449 *
450 * Compute eigenvalue condition numbers using SELECT and compare
451 *
452  CALL dcopy( n, dum, 0, stmp, 1 )
453  CALL dcopy( n, dum, 0, septmp, 1 )
454  CALL dtrsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
455  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
456  IF( info.NE.0 ) THEN
457  lmax( 3 ) = knt
458  ninfo( 3 ) = ninfo( 3 ) + 1
459  go to 240
460  END IF
461  DO 170 i = 1, n
462  IF( stmp( i ).NE.s( i ) )
463  $ vmax = one / eps
464  IF( septmp( i ).NE.dum( 1 ) )
465  $ vmax = one / eps
466  170 continue
467 *
468 * Compute eigenvector condition numbers using SELECT and compare
469 *
470  CALL dcopy( n, dum, 0, stmp, 1 )
471  CALL dcopy( n, dum, 0, septmp, 1 )
472  CALL dtrsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
473  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
474  IF( info.NE.0 ) THEN
475  lmax( 3 ) = knt
476  ninfo( 3 ) = ninfo( 3 ) + 1
477  go to 240
478  END IF
479  DO 180 i = 1, n
480  IF( stmp( i ).NE.dum( 1 ) )
481  $ vmax = one / eps
482  IF( septmp( i ).NE.sep( i ) )
483  $ vmax = one / eps
484  180 continue
485  IF( vmax.GT.rmax( 1 ) ) THEN
486  rmax( 1 ) = vmax
487  IF( ninfo( 1 ).EQ.0 )
488  $ lmax( 1 ) = knt
489  END IF
490 *
491 * Select first real and first complex eigenvalue
492 *
493  IF( wi( 1 ).EQ.zero ) THEN
494  lcmp( 1 ) = 1
495  ifnd = 0
496  DO 190 i = 2, n
497  IF( ifnd.EQ.1 .OR. wi( i ).EQ.zero ) THEN
498  SELECT( i ) = .false.
499  ELSE
500  ifnd = 1
501  lcmp( 2 ) = i
502  lcmp( 3 ) = i + 1
503  CALL dcopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
504  CALL dcopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
505  CALL dcopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
506  CALL dcopy( n, le( 1, i+1 ), 1, le( 1, 3 ), 1 )
507  END IF
508  190 continue
509  IF( ifnd.EQ.0 ) THEN
510  icmp = 1
511  ELSE
512  icmp = 3
513  END IF
514  ELSE
515  lcmp( 1 ) = 1
516  lcmp( 2 ) = 2
517  ifnd = 0
518  DO 200 i = 3, n
519  IF( ifnd.EQ.1 .OR. wi( i ).NE.zero ) THEN
520  SELECT( i ) = .false.
521  ELSE
522  lcmp( 3 ) = i
523  ifnd = 1
524  CALL dcopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
525  CALL dcopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
526  END IF
527  200 continue
528  IF( ifnd.EQ.0 ) THEN
529  icmp = 2
530  ELSE
531  icmp = 3
532  END IF
533  END IF
534 *
535 * Compute all selected condition numbers
536 *
537  CALL dcopy( icmp, dum, 0, stmp, 1 )
538  CALL dcopy( icmp, dum, 0, septmp, 1 )
539  CALL dtrsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
540  $ re, ldt, stmp, septmp, n, m, work, n, iwork,
541  $ info )
542  IF( info.NE.0 ) THEN
543  lmax( 3 ) = knt
544  ninfo( 3 ) = ninfo( 3 ) + 1
545  go to 240
546  END IF
547  DO 210 i = 1, icmp
548  j = lcmp( i )
549  IF( septmp( i ).NE.sep( j ) )
550  $ vmax = one / eps
551  IF( stmp( i ).NE.s( j ) )
552  $ vmax = one / eps
553  210 continue
554 *
555 * Compute selected eigenvalue condition numbers
556 *
557  CALL dcopy( icmp, dum, 0, stmp, 1 )
558  CALL dcopy( icmp, dum, 0, septmp, 1 )
559  CALL dtrsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
560  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
561  IF( info.NE.0 ) THEN
562  lmax( 3 ) = knt
563  ninfo( 3 ) = ninfo( 3 ) + 1
564  go to 240
565  END IF
566  DO 220 i = 1, icmp
567  j = lcmp( i )
568  IF( stmp( i ).NE.s( j ) )
569  $ vmax = one / eps
570  IF( septmp( i ).NE.dum( 1 ) )
571  $ vmax = one / eps
572  220 continue
573 *
574 * Compute selected eigenvector condition numbers
575 *
576  CALL dcopy( icmp, dum, 0, stmp, 1 )
577  CALL dcopy( icmp, dum, 0, septmp, 1 )
578  CALL dtrsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
579  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
580  IF( info.NE.0 ) THEN
581  lmax( 3 ) = knt
582  ninfo( 3 ) = ninfo( 3 ) + 1
583  go to 240
584  END IF
585  DO 230 i = 1, icmp
586  j = lcmp( i )
587  IF( stmp( i ).NE.dum( 1 ) )
588  $ vmax = one / eps
589  IF( septmp( i ).NE.sep( j ) )
590  $ vmax = one / eps
591  230 continue
592  IF( vmax.GT.rmax( 1 ) ) THEN
593  rmax( 1 ) = vmax
594  IF( ninfo( 1 ).EQ.0 )
595  $ lmax( 1 ) = knt
596  END IF
597  240 continue
598  go to 10
599 *
600 * End of DGET37
601 *
602  END