LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cget38 ( real, dimension( 3 )  RMAX,
integer, dimension( 3 )  LMAX,
integer, dimension( 3 )  NINFO,
integer  KNT,
integer  NIN 
)

CGET38

Purpose: 
 CGET38 tests CTRSEN, a routine for estimating condition numbers of a
 cluster of eigenvalues and/or its associated right invariant subspace

 The test matrices are read from a file with logical unit number NIN.
Parameters
[out]RMAX
          RMAX is REAL array, dimension (3)
          Values of the largest test ratios.
          RMAX(1) = largest residuals from CHST01 or comparing
                    different calls to CTRSEN
          RMAX(2) = largest error in reciprocal condition
                    numbers taking their conditioning into account
          RMAX(3) = largest error in reciprocal condition
                    numbers not taking their conditioning into
                    account (may be larger than RMAX(2))
[out]LMAX
          LMAX is INTEGER array, dimension (3)
          LMAX(i) is example number where largest test ratio
          RMAX(i) is achieved. Also:
          If CGEHRD returns INFO nonzero on example i, LMAX(1)=i
          If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
          If CTRSEN returns INFO nonzero on example i, LMAX(3)=i
[out]NINFO
          NINFO is INTEGER array, dimension (3)
          NINFO(1) = No. of times CGEHRD returned INFO nonzero
          NINFO(2) = No. of times CHSEQR returned INFO nonzero
          NINFO(3) = No. of times CTRSEN returned INFO nonzero
[out]KNT
          KNT is INTEGER
          Total number of examples tested.
[in]NIN
          NIN is INTEGER
          Input logical unit number.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 93 of file cget38.f.

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 *
cunghr
subroutine cunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CUNGHR
Definition: cunghr.f:128
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
chst01
subroutine chst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
CHST01
Definition: chst01.f:142
clange
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
ctrsen
subroutine ctrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
CTRSEN
Definition: ctrsen.f:266
chseqr
subroutine chseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
CHSEQR
Definition: chseqr.f:301
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
cgehrd
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
Definition: cgehrd.f:169
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: