LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ cget37()

subroutine cget37 ( real, dimension( 3 )  RMAX,
integer, dimension( 3 )  LMAX,
integer, dimension( 3 )  NINFO,
integer  KNT,
integer  NIN 
)

CGET37

Purpose:
 CGET37 tests CTRSNA, a routine for estimating condition numbers of
 eigenvalues and/or right eigenvectors of a matrix.

 The test matrices are read from a file with logical unit number NIN.
Parameters
[out]RMAX
          RMAX is REAL array, dimension (3)
          Value of the largest test ratio.
          RMAX(1) = largest ratio comparing different calls to CTRSNA
          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 CTRSNA 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 CTRSNA 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
December 2016

Definition at line 92 of file cget37.f.

92 *
93 * -- LAPACK test routine (version 3.7.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 * December 2016
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 *
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
subroutine ctrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, INFO)
CTRSNA
Definition: ctrsna.f:251
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
Definition: cgehrd.f:169
subroutine chseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
CHSEQR
Definition: chseqr.f:301
subroutine ctrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTREVC
Definition: ctrevc.f:220
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:80
Here is the call graph for this function:
Here is the caller graph for this function: