LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sget37()

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

SGET37

Purpose:
 SGET37 tests STRSNA, 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 STRSNA
          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 SGEHRD returns INFO nonzero on example i, LMAX(1)=i
          If SHSEQR returns INFO nonzero on example i, LMAX(2)=i
          If STRSNA returns INFO nonzero on example i, LMAX(3)=i
[out]NINFO
          NINFO is INTEGER array, dimension (3)
          NINFO(1) = No. of times SGEHRD returned INFO nonzero
          NINFO(2) = No. of times SHSEQR returned INFO nonzero
          NINFO(3) = No. of times STRSNA 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.

Definition at line 89 of file sget37.f.

90 *
91 * -- LAPACK test routine --
92 * -- LAPACK is a software package provided by Univ. of Tennessee, --
93 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
94 *
95 * .. Scalar Arguments ..
96  INTEGER KNT, NIN
97 * ..
98 * .. Array Arguments ..
99  INTEGER LMAX( 3 ), NINFO( 3 )
100  REAL RMAX( 3 )
101 * ..
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106  REAL ZERO, ONE, TWO
107  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
108  REAL EPSIN
109  parameter( epsin = 5.9605e-8 )
110  INTEGER LDT, LWORK
111  parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
112 * ..
113 * .. Local Scalars ..
114  INTEGER I, ICMP, IFND, INFO, ISCL, J, KMIN, M, N
115  REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116  $ VIMIN, VMAX, VMUL, VRMIN
117 * ..
118 * .. Local Arrays ..
119  LOGICAL SELECT( LDT )
120  INTEGER IWORK( 2*LDT ), LCMP( 3 )
121  REAL DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
122  $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
123  $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
124  $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
125  $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
126  $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
127  $ WRTMP( LDT )
128 * ..
129 * .. External Functions ..
130  REAL SLAMCH, SLANGE
131  EXTERNAL slamch, slange
132 * ..
133 * .. External Subroutines ..
134  EXTERNAL scopy, sgehrd, shseqr, slabad, slacpy, sscal,
135  $ strevc, strsna
136 * ..
137 * .. Intrinsic Functions ..
138  INTRINSIC max, real, sqrt
139 * ..
140 * .. Executable Statements ..
141 *
142  eps = slamch( 'P' )
143  smlnum = slamch( 'S' ) / eps
144  bignum = one / smlnum
145  CALL slabad( smlnum, bignum )
146 *
147 * EPSIN = 2**(-24) = precision to which input data computed
148 *
149  eps = max( eps, epsin )
150  rmax( 1 ) = zero
151  rmax( 2 ) = zero
152  rmax( 3 ) = zero
153  lmax( 1 ) = 0
154  lmax( 2 ) = 0
155  lmax( 3 ) = 0
156  knt = 0
157  ninfo( 1 ) = 0
158  ninfo( 2 ) = 0
159  ninfo( 3 ) = 0
160 *
161  val( 1 ) = sqrt( smlnum )
162  val( 2 ) = one
163  val( 3 ) = sqrt( bignum )
164 *
165 * Read input data until N=0. Assume input eigenvalues are sorted
166 * lexicographically (increasing by real part, then decreasing by
167 * imaginary part)
168 *
169  10 CONTINUE
170  READ( nin, fmt = * )n
171  IF( n.EQ.0 )
172  $ RETURN
173  DO 20 i = 1, n
174  READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
175  20 CONTINUE
176  DO 30 i = 1, n
177  READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
178  30 CONTINUE
179  tnrm = slange( 'M', n, n, tmp, ldt, work )
180 *
181 * Begin test
182 *
183  DO 240 iscl = 1, 3
184 *
185 * Scale input matrix
186 *
187  knt = knt + 1
188  CALL slacpy( 'F', n, n, tmp, ldt, t, ldt )
189  vmul = val( iscl )
190  DO 40 i = 1, n
191  CALL sscal( n, vmul, t( 1, i ), 1 )
192  40 CONTINUE
193  IF( tnrm.EQ.zero )
194  $ vmul = one
195 *
196 * Compute eigenvalues and eigenvectors
197 *
198  CALL sgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
199  $ info )
200  IF( info.NE.0 ) THEN
201  lmax( 1 ) = knt
202  ninfo( 1 ) = ninfo( 1 ) + 1
203  GO TO 240
204  END IF
205  DO 60 j = 1, n - 2
206  DO 50 i = j + 2, n
207  t( i, j ) = zero
208  50 CONTINUE
209  60 CONTINUE
210 *
211 * Compute Schur form
212 *
213  CALL shseqr( 'S', 'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
214  $ lwork, info )
215  IF( info.NE.0 ) THEN
216  lmax( 2 ) = knt
217  ninfo( 2 ) = ninfo( 2 ) + 1
218  GO TO 240
219  END IF
220 *
221 * Compute eigenvectors
222 *
223  CALL strevc( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
224  $ ldt, n, m, work, info )
225 *
226 * Compute condition numbers
227 *
228  CALL strsna( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
229  $ ldt, s, sep, n, m, work, n, iwork, info )
230  IF( info.NE.0 ) THEN
231  lmax( 3 ) = knt
232  ninfo( 3 ) = ninfo( 3 ) + 1
233  GO TO 240
234  END IF
235 *
236 * Sort eigenvalues and condition numbers lexicographically
237 * to compare with inputs
238 *
239  CALL scopy( n, wr, 1, wrtmp, 1 )
240  CALL scopy( n, wi, 1, witmp, 1 )
241  CALL scopy( n, s, 1, stmp, 1 )
242  CALL scopy( n, sep, 1, septmp, 1 )
243  CALL sscal( n, one / vmul, septmp, 1 )
244  DO 80 i = 1, n - 1
245  kmin = i
246  vrmin = wrtmp( i )
247  vimin = witmp( i )
248  DO 70 j = i + 1, n
249  IF( wrtmp( j ).LT.vrmin ) THEN
250  kmin = j
251  vrmin = wrtmp( j )
252  vimin = witmp( j )
253  END IF
254  70 CONTINUE
255  wrtmp( kmin ) = wrtmp( i )
256  witmp( kmin ) = witmp( i )
257  wrtmp( i ) = vrmin
258  witmp( i ) = vimin
259  vrmin = stmp( kmin )
260  stmp( kmin ) = stmp( i )
261  stmp( i ) = vrmin
262  vrmin = septmp( kmin )
263  septmp( kmin ) = septmp( i )
264  septmp( i ) = vrmin
265  80 CONTINUE
266 *
267 * Compare condition numbers for eigenvalues
268 * taking their condition numbers into account
269 *
270  v = max( two*real( n )*eps*tnrm, smlnum )
271  IF( tnrm.EQ.zero )
272  $ v = one
273  DO 90 i = 1, n
274  IF( v.GT.septmp( i ) ) THEN
275  tol = one
276  ELSE
277  tol = v / septmp( i )
278  END IF
279  IF( v.GT.sepin( i ) ) THEN
280  tolin = one
281  ELSE
282  tolin = v / sepin( i )
283  END IF
284  tol = max( tol, smlnum / eps )
285  tolin = max( tolin, smlnum / eps )
286  IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
287  vmax = one / eps
288  ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
289  vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
290  ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
291  vmax = one / eps
292  ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
293  vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
294  ELSE
295  vmax = one
296  END IF
297  IF( vmax.GT.rmax( 2 ) ) THEN
298  rmax( 2 ) = vmax
299  IF( ninfo( 2 ).EQ.0 )
300  $ lmax( 2 ) = knt
301  END IF
302  90 CONTINUE
303 *
304 * Compare condition numbers for eigenvectors
305 * taking their condition numbers into account
306 *
307  DO 100 i = 1, n
308  IF( v.GT.septmp( i )*stmp( i ) ) THEN
309  tol = septmp( i )
310  ELSE
311  tol = v / stmp( i )
312  END IF
313  IF( v.GT.sepin( i )*sin( i ) ) THEN
314  tolin = sepin( i )
315  ELSE
316  tolin = v / sin( i )
317  END IF
318  tol = max( tol, smlnum / eps )
319  tolin = max( tolin, smlnum / eps )
320  IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
321  vmax = one / eps
322  ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
323  vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
324  ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
325  vmax = one / eps
326  ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
327  vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
328  ELSE
329  vmax = one
330  END IF
331  IF( vmax.GT.rmax( 2 ) ) THEN
332  rmax( 2 ) = vmax
333  IF( ninfo( 2 ).EQ.0 )
334  $ lmax( 2 ) = knt
335  END IF
336  100 CONTINUE
337 *
338 * Compare condition numbers for eigenvalues
339 * without taking their condition numbers into account
340 *
341  DO 110 i = 1, n
342  IF( sin( i ).LE.real( 2*n )*eps .AND. stmp( i ).LE.
343  $ real( 2*n )*eps ) THEN
344  vmax = one
345  ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
346  vmax = one / eps
347  ELSE IF( sin( i ).GT.stmp( i ) ) THEN
348  vmax = sin( i ) / stmp( i )
349  ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
350  vmax = one / eps
351  ELSE IF( sin( i ).LT.stmp( i ) ) THEN
352  vmax = stmp( i ) / sin( i )
353  ELSE
354  vmax = one
355  END IF
356  IF( vmax.GT.rmax( 3 ) ) THEN
357  rmax( 3 ) = vmax
358  IF( ninfo( 3 ).EQ.0 )
359  $ lmax( 3 ) = knt
360  END IF
361  110 CONTINUE
362 *
363 * Compare condition numbers for eigenvectors
364 * without taking their condition numbers into account
365 *
366  DO 120 i = 1, n
367  IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
368  vmax = one
369  ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
370  vmax = one / eps
371  ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
372  vmax = sepin( i ) / septmp( i )
373  ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
374  vmax = one / eps
375  ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
376  vmax = septmp( i ) / sepin( i )
377  ELSE
378  vmax = one
379  END IF
380  IF( vmax.GT.rmax( 3 ) ) THEN
381  rmax( 3 ) = vmax
382  IF( ninfo( 3 ).EQ.0 )
383  $ lmax( 3 ) = knt
384  END IF
385  120 CONTINUE
386 *
387 * Compute eigenvalue condition numbers only and compare
388 *
389  vmax = zero
390  dum( 1 ) = -one
391  CALL scopy( n, dum, 0, stmp, 1 )
392  CALL scopy( n, dum, 0, septmp, 1 )
393  CALL strsna( 'Eigcond', 'All', SELECT, n, t, ldt, le, ldt, re,
394  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
395  IF( info.NE.0 ) THEN
396  lmax( 3 ) = knt
397  ninfo( 3 ) = ninfo( 3 ) + 1
398  GO TO 240
399  END IF
400  DO 130 i = 1, n
401  IF( stmp( i ).NE.s( i ) )
402  $ vmax = one / eps
403  IF( septmp( i ).NE.dum( 1 ) )
404  $ vmax = one / eps
405  130 CONTINUE
406 *
407 * Compute eigenvector condition numbers only and compare
408 *
409  CALL scopy( n, dum, 0, stmp, 1 )
410  CALL scopy( n, dum, 0, septmp, 1 )
411  CALL strsna( 'Veccond', 'All', SELECT, n, t, ldt, le, ldt, re,
412  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
413  IF( info.NE.0 ) THEN
414  lmax( 3 ) = knt
415  ninfo( 3 ) = ninfo( 3 ) + 1
416  GO TO 240
417  END IF
418  DO 140 i = 1, n
419  IF( stmp( i ).NE.dum( 1 ) )
420  $ vmax = one / eps
421  IF( septmp( i ).NE.sep( i ) )
422  $ vmax = one / eps
423  140 CONTINUE
424 *
425 * Compute all condition numbers using SELECT and compare
426 *
427  DO 150 i = 1, n
428  SELECT( i ) = .true.
429  150 CONTINUE
430  CALL scopy( n, dum, 0, stmp, 1 )
431  CALL scopy( n, dum, 0, septmp, 1 )
432  CALL strsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
433  $ re, ldt, stmp, septmp, n, m, work, n, iwork,
434  $ info )
435  IF( info.NE.0 ) THEN
436  lmax( 3 ) = knt
437  ninfo( 3 ) = ninfo( 3 ) + 1
438  GO TO 240
439  END IF
440  DO 160 i = 1, n
441  IF( septmp( i ).NE.sep( i ) )
442  $ vmax = one / eps
443  IF( stmp( i ).NE.s( i ) )
444  $ vmax = one / eps
445  160 CONTINUE
446 *
447 * Compute eigenvalue condition numbers using SELECT and compare
448 *
449  CALL scopy( n, dum, 0, stmp, 1 )
450  CALL scopy( n, dum, 0, septmp, 1 )
451  CALL strsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
452  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
453  IF( info.NE.0 ) THEN
454  lmax( 3 ) = knt
455  ninfo( 3 ) = ninfo( 3 ) + 1
456  GO TO 240
457  END IF
458  DO 170 i = 1, n
459  IF( stmp( i ).NE.s( i ) )
460  $ vmax = one / eps
461  IF( septmp( i ).NE.dum( 1 ) )
462  $ vmax = one / eps
463  170 CONTINUE
464 *
465 * Compute eigenvector condition numbers using SELECT and compare
466 *
467  CALL scopy( n, dum, 0, stmp, 1 )
468  CALL scopy( n, dum, 0, septmp, 1 )
469  CALL strsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
470  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
471  IF( info.NE.0 ) THEN
472  lmax( 3 ) = knt
473  ninfo( 3 ) = ninfo( 3 ) + 1
474  GO TO 240
475  END IF
476  DO 180 i = 1, n
477  IF( stmp( i ).NE.dum( 1 ) )
478  $ vmax = one / eps
479  IF( septmp( i ).NE.sep( i ) )
480  $ vmax = one / eps
481  180 CONTINUE
482  IF( vmax.GT.rmax( 1 ) ) THEN
483  rmax( 1 ) = vmax
484  IF( ninfo( 1 ).EQ.0 )
485  $ lmax( 1 ) = knt
486  END IF
487 *
488 * Select first real and first complex eigenvalue
489 *
490  IF( wi( 1 ).EQ.zero ) THEN
491  lcmp( 1 ) = 1
492  ifnd = 0
493  DO 190 i = 2, n
494  IF( ifnd.EQ.1 .OR. wi( i ).EQ.zero ) THEN
495  SELECT( i ) = .false.
496  ELSE
497  ifnd = 1
498  lcmp( 2 ) = i
499  lcmp( 3 ) = i + 1
500  CALL scopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
501  CALL scopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
502  CALL scopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
503  CALL scopy( n, le( 1, i+1 ), 1, le( 1, 3 ), 1 )
504  END IF
505  190 CONTINUE
506  IF( ifnd.EQ.0 ) THEN
507  icmp = 1
508  ELSE
509  icmp = 3
510  END IF
511  ELSE
512  lcmp( 1 ) = 1
513  lcmp( 2 ) = 2
514  ifnd = 0
515  DO 200 i = 3, n
516  IF( ifnd.EQ.1 .OR. wi( i ).NE.zero ) THEN
517  SELECT( i ) = .false.
518  ELSE
519  lcmp( 3 ) = i
520  ifnd = 1
521  CALL scopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
522  CALL scopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
523  END IF
524  200 CONTINUE
525  IF( ifnd.EQ.0 ) THEN
526  icmp = 2
527  ELSE
528  icmp = 3
529  END IF
530  END IF
531 *
532 * Compute all selected condition numbers
533 *
534  CALL scopy( icmp, dum, 0, stmp, 1 )
535  CALL scopy( icmp, dum, 0, septmp, 1 )
536  CALL strsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
537  $ re, ldt, stmp, septmp, n, m, work, n, iwork,
538  $ info )
539  IF( info.NE.0 ) THEN
540  lmax( 3 ) = knt
541  ninfo( 3 ) = ninfo( 3 ) + 1
542  GO TO 240
543  END IF
544  DO 210 i = 1, icmp
545  j = lcmp( i )
546  IF( septmp( i ).NE.sep( j ) )
547  $ vmax = one / eps
548  IF( stmp( i ).NE.s( j ) )
549  $ vmax = one / eps
550  210 CONTINUE
551 *
552 * Compute selected eigenvalue condition numbers
553 *
554  CALL scopy( icmp, dum, 0, stmp, 1 )
555  CALL scopy( icmp, dum, 0, septmp, 1 )
556  CALL strsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
557  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
558  IF( info.NE.0 ) THEN
559  lmax( 3 ) = knt
560  ninfo( 3 ) = ninfo( 3 ) + 1
561  GO TO 240
562  END IF
563  DO 220 i = 1, icmp
564  j = lcmp( i )
565  IF( stmp( i ).NE.s( j ) )
566  $ vmax = one / eps
567  IF( septmp( i ).NE.dum( 1 ) )
568  $ vmax = one / eps
569  220 CONTINUE
570 *
571 * Compute selected eigenvector condition numbers
572 *
573  CALL scopy( icmp, dum, 0, stmp, 1 )
574  CALL scopy( icmp, dum, 0, septmp, 1 )
575  CALL strsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
576  $ ldt, stmp, septmp, n, m, work, n, iwork, info )
577  IF( info.NE.0 ) THEN
578  lmax( 3 ) = knt
579  ninfo( 3 ) = ninfo( 3 ) + 1
580  GO TO 240
581  END IF
582  DO 230 i = 1, icmp
583  j = lcmp( i )
584  IF( stmp( i ).NE.dum( 1 ) )
585  $ vmax = one / eps
586  IF( septmp( i ).NE.sep( j ) )
587  $ vmax = one / eps
588  230 CONTINUE
589  IF( vmax.GT.rmax( 1 ) ) THEN
590  rmax( 1 ) = vmax
591  IF( ninfo( 1 ).EQ.0 )
592  $ lmax( 1 ) = knt
593  END IF
594  240 CONTINUE
595  GO TO 10
596 *
597 * End of SGET37
598 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:167
subroutine strevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STREVC
Definition: strevc.f:222
subroutine strsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
STRSNA
Definition: strsna.f:265
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:316
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: