LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cget38.f
Go to the documentation of this file.
1 *> \brief \b CGET38
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 CGET38( 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 *> CGET38 tests CTRSEN, a routine for estimating condition numbers of a
28 *> cluster of eigenvalues and/or its associated right invariant subspace
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 *> Values of the largest test ratios.
40 *> RMAX(1) = largest residuals from CHST01 or comparing
41 *> different calls to CTRSEN
42 *> RMAX(2) = largest error in reciprocal condition
43 *> numbers taking their conditioning into account
44 *> RMAX(3) = largest error in reciprocal condition
45 *> numbers not taking their conditioning into
46 *> account (may be larger than RMAX(2))
47 *> \endverbatim
48 *>
49 *> \param[out] LMAX
50 *> \verbatim
51 *> LMAX is INTEGER array, dimension (3)
52 *> LMAX(i) is example number where largest test ratio
53 *> RMAX(i) is achieved. Also:
54 *> If CGEHRD returns INFO nonzero on example i, LMAX(1)=i
55 *> If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
56 *> If CTRSEN returns INFO nonzero on example i, LMAX(3)=i
57 *> \endverbatim
58 *>
59 *> \param[out] NINFO
60 *> \verbatim
61 *> NINFO is INTEGER array, dimension (3)
62 *> NINFO(1) = No. of times CGEHRD returned INFO nonzero
63 *> NINFO(2) = No. of times CHSEQR returned INFO nonzero
64 *> NINFO(3) = No. of times CTRSEN returned INFO nonzero
65 *> \endverbatim
66 *>
67 *> \param[out] KNT
68 *> \verbatim
69 *> KNT is INTEGER
70 *> Total number of examples tested.
71 *> \endverbatim
72 *>
73 *> \param[in] NIN
74 *> \verbatim
75 *> NIN is INTEGER
76 *> Input logical unit number.
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 complex_eig
90 *
91 * =====================================================================
92  SUBROUTINE cget38( RMAX, LMAX, NINFO, KNT, NIN )
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, nin
101 * ..
102 * .. Array Arguments ..
103  INTEGER lmax( 3 ), ninfo( 3 )
104  REAL rmax( 3 )
105 * ..
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110  INTEGER ldt, lwork
111  parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
112  REAL zero, one, two
113  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
114  REAL epsin
115  parameter( epsin = 5.9605e-8 )
116  COMPLEX czero
117  parameter( czero = ( 0.0e+0, 0.0e+0 ) )
118 * ..
119 * .. Local Scalars ..
120  INTEGER i, info, iscl, isrt, itmp, j, kmin, m, n, ndim
121  REAL bignum, eps, s, sep, sepin, septmp, sin,
122  $ smlnum, stmp, tnrm, tol, tolin, v, vmax, vmin,
123  $ vmul
124 * ..
125 * .. Local Arrays ..
126  LOGICAL select( ldt )
127  INTEGER ipnt( ldt ), iselec( ldt )
128  REAL result( 2 ), rwork( ldt ), val( 3 ),
129  $ wsrt( ldt )
130  COMPLEX q( ldt, ldt ), qsav( ldt, ldt ),
131  $ qtmp( ldt, ldt ), t( ldt, ldt ),
132  $ tmp( ldt, ldt ), tsav( ldt, ldt ),
133  $ tsav1( ldt, ldt ), ttmp( ldt, ldt ), w( ldt ),
134  $ work( lwork ), wtmp( ldt )
135 * ..
136 * .. External Functions ..
137  REAL clange, slamch
138  EXTERNAL clange, slamch
139 * ..
140 * .. External Subroutines ..
141  EXTERNAL cgehrd, chseqr, chst01, clacpy, csscal, ctrsen,
142  $ cunghr, slabad
143 * ..
144 * .. Intrinsic Functions ..
145  INTRINSIC aimag, max, REAL, sqrt
146 * ..
147 * .. Executable Statements ..
148 *
149  eps = slamch( 'P' )
150  smlnum = slamch( 'S' ) / eps
151  bignum = one / smlnum
152  CALL slabad( smlnum, bignum )
153 *
154 * EPSIN = 2**(-24) = precision to which input data computed
155 *
156  eps = max( eps, epsin )
157  rmax( 1 ) = zero
158  rmax( 2 ) = zero
159  rmax( 3 ) = zero
160  lmax( 1 ) = 0
161  lmax( 2 ) = 0
162  lmax( 3 ) = 0
163  knt = 0
164  ninfo( 1 ) = 0
165  ninfo( 2 ) = 0
166  ninfo( 3 ) = 0
167  val( 1 ) = sqrt( smlnum )
168  val( 2 ) = one
169  val( 3 ) = sqrt( sqrt( bignum ) )
170 *
171 * Read input data until N=0. Assume input eigenvalues are sorted
172 * lexicographically (increasing by real part, then decreasing by
173 * imaginary part)
174 *
175  10 continue
176  READ( nin, fmt = * )n, ndim, isrt
177  IF( n.EQ.0 )
178  $ return
179  READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
180  DO 20 i = 1, n
181  READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
182  20 continue
183  READ( nin, fmt = * )sin, sepin
184 *
185  tnrm = clange( 'M', n, n, tmp, ldt, rwork )
186  DO 200 iscl = 1, 3
187 *
188 * Scale input matrix
189 *
190  knt = knt + 1
191  CALL clacpy( 'F', n, n, tmp, ldt, t, ldt )
192  vmul = val( iscl )
193  DO 30 i = 1, n
194  CALL csscal( n, vmul, t( 1, i ), 1 )
195  30 continue
196  IF( tnrm.EQ.zero )
197  $ vmul = one
198  CALL clacpy( 'F', n, n, t, ldt, tsav, ldt )
199 *
200 * Compute Schur form
201 *
202  CALL cgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
203  $ info )
204  IF( info.NE.0 ) THEN
205  lmax( 1 ) = knt
206  ninfo( 1 ) = ninfo( 1 ) + 1
207  go to 200
208  END IF
209 *
210 * Generate unitary matrix
211 *
212  CALL clacpy( 'L', n, n, t, ldt, q, ldt )
213  CALL cunghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
214  $ info )
215 *
216 * Compute Schur form
217 *
218  DO 50 j = 1, n - 2
219  DO 40 i = j + 2, n
220  t( i, j ) = czero
221  40 continue
222  50 continue
223  CALL chseqr( 'S', 'V', n, 1, n, t, ldt, w, q, ldt, work, lwork,
224  $ info )
225  IF( info.NE.0 ) THEN
226  lmax( 2 ) = knt
227  ninfo( 2 ) = ninfo( 2 ) + 1
228  go to 200
229  END IF
230 *
231 * Sort, select eigenvalues
232 *
233  DO 60 i = 1, n
234  ipnt( i ) = i
235  SELECT( i ) = .false.
236  60 continue
237  IF( isrt.EQ.0 ) THEN
238  DO 70 i = 1, n
239  wsrt( i ) = REAL( W( I ) )
240  70 continue
241  ELSE
242  DO 80 i = 1, n
243  wsrt( i ) = aimag( w( i ) )
244  80 continue
245  END IF
246  DO 100 i = 1, n - 1
247  kmin = i
248  vmin = wsrt( i )
249  DO 90 j = i + 1, n
250  IF( wsrt( j ).LT.vmin ) THEN
251  kmin = j
252  vmin = wsrt( j )
253  END IF
254  90 continue
255  wsrt( kmin ) = wsrt( i )
256  wsrt( i ) = vmin
257  itmp = ipnt( i )
258  ipnt( i ) = ipnt( kmin )
259  ipnt( kmin ) = itmp
260  100 continue
261  DO 110 i = 1, ndim
262  SELECT( ipnt( iselec( i ) ) ) = .true.
263  110 continue
264 *
265 * Compute condition numbers
266 *
267  CALL clacpy( 'F', n, n, q, ldt, qsav, ldt )
268  CALL clacpy( 'F', n, n, t, ldt, tsav1, ldt )
269  CALL ctrsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wtmp, m, s,
270  $ sep, work, lwork, info )
271  IF( info.NE.0 ) THEN
272  lmax( 3 ) = knt
273  ninfo( 3 ) = ninfo( 3 ) + 1
274  go to 200
275  END IF
276  septmp = sep / vmul
277  stmp = s
278 *
279 * Compute residuals
280 *
281  CALL chst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
282  $ rwork, result )
283  vmax = max( result( 1 ), result( 2 ) )
284  IF( vmax.GT.rmax( 1 ) ) THEN
285  rmax( 1 ) = vmax
286  IF( ninfo( 1 ).EQ.0 )
287  $ lmax( 1 ) = knt
288  END IF
289 *
290 * Compare condition number for eigenvalue cluster
291 * taking its condition number into account
292 *
293  v = max( two*REAL( n )*eps*tnrm, smlnum )
294  IF( tnrm.EQ.zero )
295  $ v = one
296  IF( v.GT.septmp ) THEN
297  tol = one
298  ELSE
299  tol = v / septmp
300  END IF
301  IF( v.GT.sepin ) THEN
302  tolin = one
303  ELSE
304  tolin = v / sepin
305  END IF
306  tol = max( tol, smlnum / eps )
307  tolin = max( tolin, smlnum / eps )
308  IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
309  vmax = one / eps
310  ELSE IF( sin-tolin.GT.stmp+tol ) THEN
311  vmax = ( sin-tolin ) / ( stmp+tol )
312  ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
313  vmax = one / eps
314  ELSE IF( sin+tolin.LT.stmp-tol ) THEN
315  vmax = ( stmp-tol ) / ( sin+tolin )
316  ELSE
317  vmax = one
318  END IF
319  IF( vmax.GT.rmax( 2 ) ) THEN
320  rmax( 2 ) = vmax
321  IF( ninfo( 2 ).EQ.0 )
322  $ lmax( 2 ) = knt
323  END IF
324 *
325 * Compare condition numbers for invariant subspace
326 * taking its condition number into account
327 *
328  IF( v.GT.septmp*stmp ) THEN
329  tol = septmp
330  ELSE
331  tol = v / stmp
332  END IF
333  IF( v.GT.sepin*sin ) THEN
334  tolin = sepin
335  ELSE
336  tolin = v / sin
337  END IF
338  tol = max( tol, smlnum / eps )
339  tolin = max( tolin, smlnum / eps )
340  IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
341  vmax = one / eps
342  ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
343  vmax = ( sepin-tolin ) / ( septmp+tol )
344  ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
345  vmax = one / eps
346  ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
347  vmax = ( septmp-tol ) / ( sepin+tolin )
348  ELSE
349  vmax = one
350  END IF
351  IF( vmax.GT.rmax( 2 ) ) THEN
352  rmax( 2 ) = vmax
353  IF( ninfo( 2 ).EQ.0 )
354  $ lmax( 2 ) = knt
355  END IF
356 *
357 * Compare condition number for eigenvalue cluster
358 * without taking its condition number into account
359 *
360  IF( sin.LE.REAL( 2*n )*eps .AND. stmp.LE.REAL( 2*n )*eps ) then
361  vmax = one
362  ELSE IF( eps*sin.GT.stmp ) THEN
363  vmax = one / eps
364  ELSE IF( sin.GT.stmp ) THEN
365  vmax = sin / stmp
366  ELSE IF( sin.LT.eps*stmp ) THEN
367  vmax = one / eps
368  ELSE IF( sin.LT.stmp ) THEN
369  vmax = stmp / sin
370  ELSE
371  vmax = one
372  END IF
373  IF( vmax.GT.rmax( 3 ) ) THEN
374  rmax( 3 ) = vmax
375  IF( ninfo( 3 ).EQ.0 )
376  $ lmax( 3 ) = knt
377  END IF
378 *
379 * Compare condition numbers for invariant subspace
380 * without taking its condition number into account
381 *
382  IF( sepin.LE.v .AND. septmp.LE.v ) THEN
383  vmax = one
384  ELSE IF( eps*sepin.GT.septmp ) THEN
385  vmax = one / eps
386  ELSE IF( sepin.GT.septmp ) THEN
387  vmax = sepin / septmp
388  ELSE IF( sepin.LT.eps*septmp ) THEN
389  vmax = one / eps
390  ELSE IF( sepin.LT.septmp ) THEN
391  vmax = septmp / sepin
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 *
401 * Compute eigenvalue condition number only and compare
402 * Update Q
403 *
404  vmax = zero
405  CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
406  CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
407  septmp = -one
408  stmp = -one
409  CALL ctrsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
410  $ m, stmp, septmp, work, lwork, info )
411  IF( info.NE.0 ) THEN
412  lmax( 3 ) = knt
413  ninfo( 3 ) = ninfo( 3 ) + 1
414  go to 200
415  END IF
416  IF( s.NE.stmp )
417  $ vmax = one / eps
418  IF( -one.NE.septmp )
419  $ vmax = one / eps
420  DO 130 i = 1, n
421  DO 120 j = 1, n
422  IF( ttmp( i, j ).NE.t( i, j ) )
423  $ vmax = one / eps
424  IF( qtmp( i, j ).NE.q( i, j ) )
425  $ vmax = one / eps
426  120 continue
427  130 continue
428 *
429 * Compute invariant subspace condition number only and compare
430 * Update Q
431 *
432  CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
433  CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
434  septmp = -one
435  stmp = -one
436  CALL ctrsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
437  $ m, stmp, septmp, work, lwork, info )
438  IF( info.NE.0 ) THEN
439  lmax( 3 ) = knt
440  ninfo( 3 ) = ninfo( 3 ) + 1
441  go to 200
442  END IF
443  IF( -one.NE.stmp )
444  $ vmax = one / eps
445  IF( sep.NE.septmp )
446  $ vmax = one / eps
447  DO 150 i = 1, n
448  DO 140 j = 1, n
449  IF( ttmp( i, j ).NE.t( i, j ) )
450  $ vmax = one / eps
451  IF( qtmp( i, j ).NE.q( i, j ) )
452  $ vmax = one / eps
453  140 continue
454  150 continue
455 *
456 * Compute eigenvalue condition number only and compare
457 * Do not update Q
458 *
459  CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
460  CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
461  septmp = -one
462  stmp = -one
463  CALL ctrsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
464  $ m, stmp, septmp, work, lwork, info )
465  IF( info.NE.0 ) THEN
466  lmax( 3 ) = knt
467  ninfo( 3 ) = ninfo( 3 ) + 1
468  go to 200
469  END IF
470  IF( s.NE.stmp )
471  $ vmax = one / eps
472  IF( -one.NE.septmp )
473  $ vmax = one / eps
474  DO 170 i = 1, n
475  DO 160 j = 1, n
476  IF( ttmp( i, j ).NE.t( i, j ) )
477  $ vmax = one / eps
478  IF( qtmp( i, j ).NE.qsav( i, j ) )
479  $ vmax = one / eps
480  160 continue
481  170 continue
482 *
483 * Compute invariant subspace condition number only and compare
484 * Do not update Q
485 *
486  CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
487  CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
488  septmp = -one
489  stmp = -one
490  CALL ctrsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
491  $ m, stmp, septmp, work, lwork, info )
492  IF( info.NE.0 ) THEN
493  lmax( 3 ) = knt
494  ninfo( 3 ) = ninfo( 3 ) + 1
495  go to 200
496  END IF
497  IF( -one.NE.stmp )
498  $ vmax = one / eps
499  IF( sep.NE.septmp )
500  $ vmax = one / eps
501  DO 190 i = 1, n
502  DO 180 j = 1, n
503  IF( ttmp( i, j ).NE.t( i, j ) )
504  $ vmax = one / eps
505  IF( qtmp( i, j ).NE.qsav( i, j ) )
506  $ vmax = one / eps
507  180 continue
508  190 continue
509  IF( vmax.GT.rmax( 1 ) ) THEN
510  rmax( 1 ) = vmax
511  IF( ninfo( 1 ).EQ.0 )
512  $ lmax( 1 ) = knt
513  END IF
514  200 continue
515  go to 10
516 *
517 * End of CGET38
518 *
519  END