LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cget37.f
Go to the documentation of this file.
1 *> \brief \b CGET37
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 CGET37( RMAX, LMAX, NINFO, KNT, NIN )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER KNT, NIN
15 * ..
16 * .. Array Arguments ..
17 * INTEGER LMAX( 3 ), NINFO( 3 )
18 * REAL RMAX( 3 )
19 * ..
20 *
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> CGET37 tests CTRSNA, 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 REAL array, dimension (3)
39 *> Value of the largest test ratio.
40 *> RMAX(1) = largest ratio comparing different calls to CTRSNA
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 CGEHRD returns INFO nonzero on example i, LMAX(1)=i
54 *> If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
55 *> If CTRSNA 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 CGEHRD returned INFO nonzero
62 *> NINFO(2) = No. of times CHSEQR returned INFO nonzero
63 *> NINFO(3) = No. of times CTRSNA 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 complex_eig
89 *
90 * =====================================================================
91  SUBROUTINE cget37( 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  REAL rmax( 3 )
104 * ..
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109  REAL zero, one, two
110  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
111  REAL epsin
112  parameter( epsin = 5.9605e-8 )
113  INTEGER ldt, lwork
114  parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
115 * ..
116 * .. Local Scalars ..
117  INTEGER i, icmp, info, iscl, isrt, j, kmin, m, n
118  REAL bignum, eps, smlnum, tnrm, tol, tolin, v,
119  $ vcmin, vmax, vmin, vmul
120 * ..
121 * .. Local Arrays ..
122  LOGICAL select( ldt )
123  INTEGER lcmp( 3 )
124  REAL dum( 1 ), rwork( 2*ldt ), s( ldt ), sep( ldt ),
125  $ sepin( ldt ), septmp( ldt ), sin( ldt ),
126  $ stmp( ldt ), val( 3 ), wiin( ldt ),
127  $ wrin( ldt ), wsrt( ldt )
128  COMPLEX cdum( 1 ), le( ldt, ldt ), re( ldt, ldt ),
129  $ t( ldt, ldt ), tmp( ldt, ldt ), w( ldt ),
130  $ work( lwork ), wtmp( ldt )
131 * ..
132 * .. External Functions ..
133  REAL clange, slamch
134  EXTERNAL clange, slamch
135 * ..
136 * .. External Subroutines ..
137  EXTERNAL ccopy, cgehrd, chseqr, clacpy, csscal, ctrevc,
138  $ ctrsna, scopy, slabad, sscal
139 * ..
140 * .. Intrinsic Functions ..
141  INTRINSIC aimag, max, REAL, sqrt
142 * ..
143 * .. Executable Statements ..
144 *
145  eps = slamch( 'P' )
146  smlnum = slamch( 'S' ) / eps
147  bignum = one / smlnum
148  CALL slabad( 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  val( 1 ) = sqrt( smlnum )
164  val( 2 ) = one
165  val( 3 ) = sqrt( bignum )
166 *
167 * Read input data until N=0. Assume input eigenvalues are sorted
168 * lexicographically (increasing by real part if ISRT = 0,
169 * increasing by imaginary part if ISRT = 1)
170 *
171  10 continue
172  READ( nin, fmt = * )n, isrt
173  IF( n.EQ.0 )
174  $ return
175  DO 20 i = 1, n
176  READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
177  20 continue
178  DO 30 i = 1, n
179  READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
180  30 continue
181  tnrm = clange( 'M', n, n, tmp, ldt, rwork )
182  DO 260 iscl = 1, 3
183 *
184 * Scale input matrix
185 *
186  knt = knt + 1
187  CALL clacpy( 'F', n, n, tmp, ldt, t, ldt )
188  vmul = val( iscl )
189  DO 40 i = 1, n
190  CALL csscal( n, vmul, t( 1, i ), 1 )
191  40 continue
192  IF( tnrm.EQ.zero )
193  $ vmul = one
194 *
195 * Compute eigenvalues and eigenvectors
196 *
197  CALL cgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
198  $ info )
199  IF( info.NE.0 ) THEN
200  lmax( 1 ) = knt
201  ninfo( 1 ) = ninfo( 1 ) + 1
202  go to 260
203  END IF
204  DO 60 j = 1, n - 2
205  DO 50 i = j + 2, n
206  t( i, j ) = zero
207  50 continue
208  60 continue
209 *
210 * Compute Schur form
211 *
212  CALL chseqr( 'S', 'N', n, 1, n, t, ldt, w, cdum, 1, work,
213  $ lwork, info )
214  IF( info.NE.0 ) THEN
215  lmax( 2 ) = knt
216  ninfo( 2 ) = ninfo( 2 ) + 1
217  go to 260
218  END IF
219 *
220 * Compute eigenvectors
221 *
222  DO 70 i = 1, n
223  SELECT( i ) = .true.
224  70 continue
225  CALL ctrevc( 'B', 'A', SELECT, n, t, ldt, le, ldt, re, ldt, n,
226  $ m, work, rwork, info )
227 *
228 * Compute condition numbers
229 *
230  CALL ctrsna( 'B', 'A', SELECT, n, t, ldt, le, ldt, re, ldt, s,
231  $ sep, n, m, work, n, rwork, info )
232  IF( info.NE.0 ) THEN
233  lmax( 3 ) = knt
234  ninfo( 3 ) = ninfo( 3 ) + 1
235  go to 260
236  END IF
237 *
238 * Sort eigenvalues and condition numbers lexicographically
239 * to compare with inputs
240 *
241  CALL ccopy( n, w, 1, wtmp, 1 )
242  IF( isrt.EQ.0 ) THEN
243 *
244 * Sort by increasing real part
245 *
246  DO 80 i = 1, n
247  wsrt( i ) = REAL( W( I ) )
248  80 continue
249  ELSE
250 *
251 * Sort by increasing imaginary part
252 *
253  DO 90 i = 1, n
254  wsrt( i ) = aimag( w( i ) )
255  90 continue
256  END IF
257  CALL scopy( n, s, 1, stmp, 1 )
258  CALL scopy( n, sep, 1, septmp, 1 )
259  CALL sscal( n, one / vmul, septmp, 1 )
260  DO 110 i = 1, n - 1
261  kmin = i
262  vmin = wsrt( i )
263  DO 100 j = i + 1, n
264  IF( wsrt( j ).LT.vmin ) THEN
265  kmin = j
266  vmin = wsrt( j )
267  END IF
268  100 continue
269  wsrt( kmin ) = wsrt( i )
270  wsrt( i ) = vmin
271  vcmin = wtmp( i )
272  wtmp( i ) = w( kmin )
273  wtmp( kmin ) = vcmin
274  vmin = stmp( kmin )
275  stmp( kmin ) = stmp( i )
276  stmp( i ) = vmin
277  vmin = septmp( kmin )
278  septmp( kmin ) = septmp( i )
279  septmp( i ) = vmin
280  110 continue
281 *
282 * Compare condition numbers for eigenvalues
283 * taking their condition numbers into account
284 *
285  v = max( two*REAL( n )*eps*tnrm, smlnum )
286  IF( tnrm.EQ.zero )
287  $ v = one
288  DO 120 i = 1, n
289  IF( v.GT.septmp( i ) ) THEN
290  tol = one
291  ELSE
292  tol = v / septmp( i )
293  END IF
294  IF( v.GT.sepin( i ) ) THEN
295  tolin = one
296  ELSE
297  tolin = v / sepin( i )
298  END IF
299  tol = max( tol, smlnum / eps )
300  tolin = max( tolin, smlnum / eps )
301  IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
302  vmax = one / eps
303  ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
304  vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
305  ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
306  vmax = one / eps
307  ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
308  vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
309  ELSE
310  vmax = one
311  END IF
312  IF( vmax.GT.rmax( 2 ) ) THEN
313  rmax( 2 ) = vmax
314  IF( ninfo( 2 ).EQ.0 )
315  $ lmax( 2 ) = knt
316  END IF
317  120 continue
318 *
319 * Compare condition numbers for eigenvectors
320 * taking their condition numbers into account
321 *
322  DO 130 i = 1, n
323  IF( v.GT.septmp( i )*stmp( i ) ) THEN
324  tol = septmp( i )
325  ELSE
326  tol = v / stmp( i )
327  END IF
328  IF( v.GT.sepin( i )*sin( i ) ) THEN
329  tolin = sepin( i )
330  ELSE
331  tolin = v / sin( i )
332  END IF
333  tol = max( tol, smlnum / eps )
334  tolin = max( tolin, smlnum / eps )
335  IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
336  vmax = one / eps
337  ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
338  vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
339  ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
340  vmax = one / eps
341  ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
342  vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
343  ELSE
344  vmax = one
345  END IF
346  IF( vmax.GT.rmax( 2 ) ) THEN
347  rmax( 2 ) = vmax
348  IF( ninfo( 2 ).EQ.0 )
349  $ lmax( 2 ) = knt
350  END IF
351  130 continue
352 *
353 * Compare condition numbers for eigenvalues
354 * without taking their condition numbers into account
355 *
356  DO 140 i = 1, n
357  IF( sin( i ).LE.REAL( 2*n )*eps .AND. stmp( i ).LE.
358  $ REAL( 2*n )*eps ) then
359  vmax = one
360  ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
361  vmax = one / eps
362  ELSE IF( sin( i ).GT.stmp( i ) ) THEN
363  vmax = sin( i ) / stmp( i )
364  ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
365  vmax = one / eps
366  ELSE IF( sin( i ).LT.stmp( i ) ) THEN
367  vmax = stmp( i ) / sin( i )
368  ELSE
369  vmax = one
370  END IF
371  IF( vmax.GT.rmax( 3 ) ) THEN
372  rmax( 3 ) = vmax
373  IF( ninfo( 3 ).EQ.0 )
374  $ lmax( 3 ) = knt
375  END IF
376  140 continue
377 *
378 * Compare condition numbers for eigenvectors
379 * without taking their condition numbers into account
380 *
381  DO 150 i = 1, n
382  IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
383  vmax = one
384  ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
385  vmax = one / eps
386  ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
387  vmax = sepin( i ) / septmp( i )
388  ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
389  vmax = one / eps
390  ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
391  vmax = septmp( i ) / sepin( i )
392  ELSE
393  vmax = one
394  END IF
395  IF( vmax.GT.rmax( 3 ) ) THEN
396  rmax( 3 ) = vmax
397  IF( ninfo( 3 ).EQ.0 )
398  $ lmax( 3 ) = knt
399  END IF
400  150 continue
401 *
402 * Compute eigenvalue condition numbers only and compare
403 *
404  vmax = zero
405  dum( 1 ) = -one
406  CALL scopy( n, dum, 0, stmp, 1 )
407  CALL scopy( n, dum, 0, septmp, 1 )
408  CALL ctrsna( 'E', 'A', SELECT, n, t, ldt, le, ldt, re, ldt,
409  $ stmp, septmp, n, m, work, n, rwork, info )
410  IF( info.NE.0 ) THEN
411  lmax( 3 ) = knt
412  ninfo( 3 ) = ninfo( 3 ) + 1
413  go to 260
414  END IF
415  DO 160 i = 1, n
416  IF( stmp( i ).NE.s( i ) )
417  $ vmax = one / eps
418  IF( septmp( i ).NE.dum( 1 ) )
419  $ vmax = one / eps
420  160 continue
421 *
422 * Compute eigenvector condition numbers only and compare
423 *
424  CALL scopy( n, dum, 0, stmp, 1 )
425  CALL scopy( n, dum, 0, septmp, 1 )
426  CALL ctrsna( 'V', 'A', SELECT, n, t, ldt, le, ldt, re, ldt,
427  $ stmp, septmp, n, m, work, n, rwork, info )
428  IF( info.NE.0 ) THEN
429  lmax( 3 ) = knt
430  ninfo( 3 ) = ninfo( 3 ) + 1
431  go to 260
432  END IF
433  DO 170 i = 1, n
434  IF( stmp( i ).NE.dum( 1 ) )
435  $ vmax = one / eps
436  IF( septmp( i ).NE.sep( i ) )
437  $ vmax = one / eps
438  170 continue
439 *
440 * Compute all condition numbers using SELECT and compare
441 *
442  DO 180 i = 1, n
443  SELECT( i ) = .true.
444  180 continue
445  CALL scopy( n, dum, 0, stmp, 1 )
446  CALL scopy( n, dum, 0, septmp, 1 )
447  CALL ctrsna( 'B', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
448  $ stmp, septmp, n, m, work, n, rwork, info )
449  IF( info.NE.0 ) THEN
450  lmax( 3 ) = knt
451  ninfo( 3 ) = ninfo( 3 ) + 1
452  go to 260
453  END IF
454  DO 190 i = 1, n
455  IF( septmp( i ).NE.sep( i ) )
456  $ vmax = one / eps
457  IF( stmp( i ).NE.s( i ) )
458  $ vmax = one / eps
459  190 continue
460 *
461 * Compute eigenvalue condition numbers using SELECT and compare
462 *
463  CALL scopy( n, dum, 0, stmp, 1 )
464  CALL scopy( n, dum, 0, septmp, 1 )
465  CALL ctrsna( 'E', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
466  $ stmp, septmp, n, m, work, n, rwork, info )
467  IF( info.NE.0 ) THEN
468  lmax( 3 ) = knt
469  ninfo( 3 ) = ninfo( 3 ) + 1
470  go to 260
471  END IF
472  DO 200 i = 1, n
473  IF( stmp( i ).NE.s( i ) )
474  $ vmax = one / eps
475  IF( septmp( i ).NE.dum( 1 ) )
476  $ vmax = one / eps
477  200 continue
478 *
479 * Compute eigenvector condition numbers using SELECT and compare
480 *
481  CALL scopy( n, dum, 0, stmp, 1 )
482  CALL scopy( n, dum, 0, septmp, 1 )
483  CALL ctrsna( 'V', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
484  $ stmp, septmp, n, m, work, n, rwork, info )
485  IF( info.NE.0 ) THEN
486  lmax( 3 ) = knt
487  ninfo( 3 ) = ninfo( 3 ) + 1
488  go to 260
489  END IF
490  DO 210 i = 1, n
491  IF( stmp( i ).NE.dum( 1 ) )
492  $ vmax = one / eps
493  IF( septmp( i ).NE.sep( i ) )
494  $ vmax = one / eps
495  210 continue
496  IF( vmax.GT.rmax( 1 ) ) THEN
497  rmax( 1 ) = vmax
498  IF( ninfo( 1 ).EQ.0 )
499  $ lmax( 1 ) = knt
500  END IF
501 *
502 * Select second and next to last eigenvalues
503 *
504  DO 220 i = 1, n
505  SELECT( i ) = .false.
506  220 continue
507  icmp = 0
508  IF( n.GT.1 ) THEN
509  icmp = 1
510  lcmp( 1 ) = 2
511  SELECT( 2 ) = .true.
512  CALL ccopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
513  CALL ccopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
514  END IF
515  IF( n.GT.3 ) THEN
516  icmp = 2
517  lcmp( 2 ) = n - 1
518  SELECT( n-1 ) = .true.
519  CALL ccopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
520  CALL ccopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
521  END IF
522 *
523 * Compute all selected condition numbers
524 *
525  CALL scopy( icmp, dum, 0, stmp, 1 )
526  CALL scopy( icmp, dum, 0, septmp, 1 )
527  CALL ctrsna( 'B', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
528  $ stmp, septmp, n, m, work, n, rwork, info )
529  IF( info.NE.0 ) THEN
530  lmax( 3 ) = knt
531  ninfo( 3 ) = ninfo( 3 ) + 1
532  go to 260
533  END IF
534  DO 230 i = 1, icmp
535  j = lcmp( i )
536  IF( septmp( i ).NE.sep( j ) )
537  $ vmax = one / eps
538  IF( stmp( i ).NE.s( j ) )
539  $ vmax = one / eps
540  230 continue
541 *
542 * Compute selected eigenvalue condition numbers
543 *
544  CALL scopy( icmp, dum, 0, stmp, 1 )
545  CALL scopy( icmp, dum, 0, septmp, 1 )
546  CALL ctrsna( 'E', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
547  $ stmp, septmp, n, m, work, n, rwork, info )
548  IF( info.NE.0 ) THEN
549  lmax( 3 ) = knt
550  ninfo( 3 ) = ninfo( 3 ) + 1
551  go to 260
552  END IF
553  DO 240 i = 1, icmp
554  j = lcmp( i )
555  IF( stmp( i ).NE.s( j ) )
556  $ vmax = one / eps
557  IF( septmp( i ).NE.dum( 1 ) )
558  $ vmax = one / eps
559  240 continue
560 *
561 * Compute selected eigenvector condition numbers
562 *
563  CALL scopy( icmp, dum, 0, stmp, 1 )
564  CALL scopy( icmp, dum, 0, septmp, 1 )
565  CALL ctrsna( 'V', 'S', SELECT, n, t, ldt, le, ldt, re, ldt,
566  $ stmp, septmp, n, m, work, n, rwork, info )
567  IF( info.NE.0 ) THEN
568  lmax( 3 ) = knt
569  ninfo( 3 ) = ninfo( 3 ) + 1
570  go to 260
571  END IF
572  DO 250 i = 1, icmp
573  j = lcmp( i )
574  IF( stmp( i ).NE.dum( 1 ) )
575  $ vmax = one / eps
576  IF( septmp( i ).NE.sep( j ) )
577  $ vmax = one / eps
578  250 continue
579  IF( vmax.GT.rmax( 1 ) ) THEN
580  rmax( 1 ) = vmax
581  IF( ninfo( 1 ).EQ.0 )
582  $ lmax( 1 ) = knt
583  END IF
584  260 continue
585  go to 10
586 *
587 * End of CGET37
588 *
589  END