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

SGET38

Purpose: 
 SGET38 tests STRSEN, 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 SHST01 or comparing
                    different calls to STRSEN
          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 STRSEN 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 STRSEN 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 sget38.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  REAL zero, one, two
111  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
112  REAL epsin
113  parameter ( epsin = 5.9605e-8 )
114  INTEGER ldt, lwork
115  parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
116  INTEGER liwork
117  parameter ( liwork = ldt*ldt )
118 * ..
119 * .. Local Scalars ..
120  INTEGER i, info, iscl, itmp, j, kmin, m, n, ndim
121  REAL bignum, eps, s, sep, sepin, septmp, sin,
122  $ smlnum, stmp, tnrm, tol, tolin, v, vimin, vmax,
123  $ vmul, vrmin
124 * ..
125 * .. Local Arrays ..
126  LOGICAL select( ldt )
127  INTEGER ipnt( ldt ), iselec( ldt ), iwork( liwork )
128  REAL q( ldt, ldt ), qsav( ldt, ldt ),
129  $ qtmp( ldt, ldt ), result( 2 ), t( ldt, ldt ),
130  $ tmp( ldt, ldt ), tsav( ldt, ldt ),
131  $ tsav1( ldt, ldt ), ttmp( ldt, ldt ), val( 3 ),
132  $ wi( ldt ), witmp( ldt ), work( lwork ),
133  $ wr( ldt ), wrtmp( ldt )
134 * ..
135 * .. External Functions ..
136  REAL slamch, slange
137  EXTERNAL slamch, slange
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL scopy, sgehrd, shseqr, shst01, slabad, slacpy,
141  $ sorghr, sscal, strsen
142 * ..
143 * .. Intrinsic Functions ..
144  INTRINSIC max, REAL, sqrt
145 * ..
146 * .. Executable Statements ..
147 *
148  eps = slamch( 'P' )
149  smlnum = slamch( 'S' ) / eps
150  bignum = one / smlnum
151  CALL slabad( smlnum, bignum )
152 *
153 * EPSIN = 2**(-24) = precision to which input data computed
154 *
155  eps = max( eps, epsin )
156  rmax( 1 ) = zero
157  rmax( 2 ) = zero
158  rmax( 3 ) = zero
159  lmax( 1 ) = 0
160  lmax( 2 ) = 0
161  lmax( 3 ) = 0
162  knt = 0
163  ninfo( 1 ) = 0
164  ninfo( 2 ) = 0
165  ninfo( 3 ) = 0
166 *
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
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 = slange( 'M', n, n, tmp, ldt, work )
186  DO 160 iscl = 1, 3
187 *
188 * Scale input matrix
189 *
190  knt = knt + 1
191  CALL slacpy( 'F', n, n, tmp, ldt, t, ldt )
192  vmul = val( iscl )
193  DO 30 i = 1, n
194  CALL sscal( n, vmul, t( 1, i ), 1 )
195  30 CONTINUE
196  IF( tnrm.EQ.zero )
197  $ vmul = one
198  CALL slacpy( 'F', n, n, t, ldt, tsav, ldt )
199 *
200 * Compute Schur form
201 *
202  CALL sgehrd( 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 160
208  END IF
209 *
210 * Generate orthogonal matrix
211 *
212  CALL slacpy( 'L', n, n, t, ldt, q, ldt )
213  CALL sorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
214  $ info )
215 *
216 * Compute Schur form
217 *
218  CALL shseqr( 'S', 'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
219  $ lwork, info )
220  IF( info.NE.0 ) THEN
221  lmax( 2 ) = knt
222  ninfo( 2 ) = ninfo( 2 ) + 1
223  GO TO 160
224  END IF
225 *
226 * Sort, select eigenvalues
227 *
228  DO 40 i = 1, n
229  ipnt( i ) = i
230  SELECT( i ) = .false.
231  40 CONTINUE
232  CALL scopy( n, wr, 1, wrtmp, 1 )
233  CALL scopy( n, wi, 1, witmp, 1 )
234  DO 60 i = 1, n - 1
235  kmin = i
236  vrmin = wrtmp( i )
237  vimin = witmp( i )
238  DO 50 j = i + 1, n
239  IF( wrtmp( j ).LT.vrmin ) THEN
240  kmin = j
241  vrmin = wrtmp( j )
242  vimin = witmp( j )
243  END IF
244  50 CONTINUE
245  wrtmp( kmin ) = wrtmp( i )
246  witmp( kmin ) = witmp( i )
247  wrtmp( i ) = vrmin
248  witmp( i ) = vimin
249  itmp = ipnt( i )
250  ipnt( i ) = ipnt( kmin )
251  ipnt( kmin ) = itmp
252  60 CONTINUE
253  DO 70 i = 1, ndim
254  SELECT( ipnt( iselec( i ) ) ) = .true.
255  70 CONTINUE
256 *
257 * Compute condition numbers
258 *
259  CALL slacpy( 'F', n, n, q, ldt, qsav, ldt )
260  CALL slacpy( 'F', n, n, t, ldt, tsav1, ldt )
261  CALL strsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wrtmp, witmp,
262  $ m, s, sep, work, lwork, iwork, liwork, info )
263  IF( info.NE.0 ) THEN
264  lmax( 3 ) = knt
265  ninfo( 3 ) = ninfo( 3 ) + 1
266  GO TO 160
267  END IF
268  septmp = sep / vmul
269  stmp = s
270 *
271 * Compute residuals
272 *
273  CALL shst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
274  $ result )
275  vmax = max( result( 1 ), result( 2 ) )
276  IF( vmax.GT.rmax( 1 ) ) THEN
277  rmax( 1 ) = vmax
278  IF( ninfo( 1 ).EQ.0 )
279  $ lmax( 1 ) = knt
280  END IF
281 *
282 * Compare condition number for eigenvalue cluster
283 * taking its condition number into account
284 *
285  v = max( two*REAL( n )*eps*tnrm, smlnum )
286  IF( tnrm.EQ.zero )
287  $ v = one
288  IF( v.GT.septmp ) THEN
289  tol = one
290  ELSE
291  tol = v / septmp
292  END IF
293  IF( v.GT.sepin ) THEN
294  tolin = one
295  ELSE
296  tolin = v / sepin
297  END IF
298  tol = max( tol, smlnum / eps )
299  tolin = max( tolin, smlnum / eps )
300  IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
301  vmax = one / eps
302  ELSE IF( sin-tolin.GT.stmp+tol ) THEN
303  vmax = ( sin-tolin ) / ( stmp+tol )
304  ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
305  vmax = one / eps
306  ELSE IF( sin+tolin.LT.stmp-tol ) THEN
307  vmax = ( stmp-tol ) / ( sin+tolin )
308  ELSE
309  vmax = one
310  END IF
311  IF( vmax.GT.rmax( 2 ) ) THEN
312  rmax( 2 ) = vmax
313  IF( ninfo( 2 ).EQ.0 )
314  $ lmax( 2 ) = knt
315  END IF
316 *
317 * Compare condition numbers for invariant subspace
318 * taking its condition number into account
319 *
320  IF( v.GT.septmp*stmp ) THEN
321  tol = septmp
322  ELSE
323  tol = v / stmp
324  END IF
325  IF( v.GT.sepin*sin ) THEN
326  tolin = sepin
327  ELSE
328  tolin = v / sin
329  END IF
330  tol = max( tol, smlnum / eps )
331  tolin = max( tolin, smlnum / eps )
332  IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
333  vmax = one / eps
334  ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
335  vmax = ( sepin-tolin ) / ( septmp+tol )
336  ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
337  vmax = one / eps
338  ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
339  vmax = ( septmp-tol ) / ( sepin+tolin )
340  ELSE
341  vmax = one
342  END IF
343  IF( vmax.GT.rmax( 2 ) ) THEN
344  rmax( 2 ) = vmax
345  IF( ninfo( 2 ).EQ.0 )
346  $ lmax( 2 ) = knt
347  END IF
348 *
349 * Compare condition number for eigenvalue cluster
350 * without taking its condition number into account
351 *
352  IF( sin.LE.REAL( 2*n )*eps .AND. stmp.LE.REAL( 2*n )*eps ) then
353  vmax = one
354  ELSE IF( eps*sin.GT.stmp ) THEN
355  vmax = one / eps
356  ELSE IF( sin.GT.stmp ) THEN
357  vmax = sin / stmp
358  ELSE IF( sin.LT.eps*stmp ) THEN
359  vmax = one / eps
360  ELSE IF( sin.LT.stmp ) THEN
361  vmax = stmp / sin
362  ELSE
363  vmax = one
364  END IF
365  IF( vmax.GT.rmax( 3 ) ) THEN
366  rmax( 3 ) = vmax
367  IF( ninfo( 3 ).EQ.0 )
368  $ lmax( 3 ) = knt
369  END IF
370 *
371 * Compare condition numbers for invariant subspace
372 * without taking its condition number into account
373 *
374  IF( sepin.LE.v .AND. septmp.LE.v ) THEN
375  vmax = one
376  ELSE IF( eps*sepin.GT.septmp ) THEN
377  vmax = one / eps
378  ELSE IF( sepin.GT.septmp ) THEN
379  vmax = sepin / septmp
380  ELSE IF( sepin.LT.eps*septmp ) THEN
381  vmax = one / eps
382  ELSE IF( sepin.LT.septmp ) THEN
383  vmax = septmp / sepin
384  ELSE
385  vmax = one
386  END IF
387  IF( vmax.GT.rmax( 3 ) ) THEN
388  rmax( 3 ) = vmax
389  IF( ninfo( 3 ).EQ.0 )
390  $ lmax( 3 ) = knt
391  END IF
392 *
393 * Compute eigenvalue condition number only and compare
394 * Update Q
395 *
396  vmax = zero
397  CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
398  CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
399  septmp = -one
400  stmp = -one
401  CALL strsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
402  $ witmp, m, stmp, septmp, work, lwork, iwork,
403  $ liwork, info )
404  IF( info.NE.0 ) THEN
405  lmax( 3 ) = knt
406  ninfo( 3 ) = ninfo( 3 ) + 1
407  GO TO 160
408  END IF
409  IF( s.NE.stmp )
410  $ vmax = one / eps
411  IF( -one.NE.septmp )
412  $ vmax = one / eps
413  DO 90 i = 1, n
414  DO 80 j = 1, n
415  IF( ttmp( i, j ).NE.t( i, j ) )
416  $ vmax = one / eps
417  IF( qtmp( i, j ).NE.q( i, j ) )
418  $ vmax = one / eps
419  80 CONTINUE
420  90 CONTINUE
421 *
422 * Compute invariant subspace condition number only and compare
423 * Update Q
424 *
425  CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
426  CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
427  septmp = -one
428  stmp = -one
429  CALL strsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
430  $ witmp, m, stmp, septmp, work, lwork, iwork,
431  $ liwork, info )
432  IF( info.NE.0 ) THEN
433  lmax( 3 ) = knt
434  ninfo( 3 ) = ninfo( 3 ) + 1
435  GO TO 160
436  END IF
437  IF( -one.NE.stmp )
438  $ vmax = one / eps
439  IF( sep.NE.septmp )
440  $ vmax = one / eps
441  DO 110 i = 1, n
442  DO 100 j = 1, n
443  IF( ttmp( i, j ).NE.t( i, j ) )
444  $ vmax = one / eps
445  IF( qtmp( i, j ).NE.q( i, j ) )
446  $ vmax = one / eps
447  100 CONTINUE
448  110 CONTINUE
449 *
450 * Compute eigenvalue condition number only and compare
451 * Do not update Q
452 *
453  CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
454  CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
455  septmp = -one
456  stmp = -one
457  CALL strsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
458  $ witmp, m, stmp, septmp, work, lwork, iwork,
459  $ liwork, info )
460  IF( info.NE.0 ) THEN
461  lmax( 3 ) = knt
462  ninfo( 3 ) = ninfo( 3 ) + 1
463  GO TO 160
464  END IF
465  IF( s.NE.stmp )
466  $ vmax = one / eps
467  IF( -one.NE.septmp )
468  $ vmax = one / eps
469  DO 130 i = 1, n
470  DO 120 j = 1, n
471  IF( ttmp( i, j ).NE.t( i, j ) )
472  $ vmax = one / eps
473  IF( qtmp( i, j ).NE.qsav( i, j ) )
474  $ vmax = one / eps
475  120 CONTINUE
476  130 CONTINUE
477 *
478 * Compute invariant subspace condition number only and compare
479 * Do not update Q
480 *
481  CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
482  CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
483  septmp = -one
484  stmp = -one
485  CALL strsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
486  $ witmp, m, stmp, septmp, work, lwork, iwork,
487  $ liwork, info )
488  IF( info.NE.0 ) THEN
489  lmax( 3 ) = knt
490  ninfo( 3 ) = ninfo( 3 ) + 1
491  GO TO 160
492  END IF
493  IF( -one.NE.stmp )
494  $ vmax = one / eps
495  IF( sep.NE.septmp )
496  $ vmax = one / eps
497  DO 150 i = 1, n
498  DO 140 j = 1, n
499  IF( ttmp( i, j ).NE.t( i, j ) )
500  $ vmax = one / eps
501  IF( qtmp( i, j ).NE.qsav( i, j ) )
502  $ vmax = one / eps
503  140 CONTINUE
504  150 CONTINUE
505  IF( vmax.GT.rmax( 1 ) ) THEN
506  rmax( 1 ) = vmax
507  IF( ninfo( 1 ).EQ.0 )
508  $ lmax( 1 ) = knt
509  END IF
510  160 CONTINUE
511  GO TO 10
512 *
513 * End of SGET38
514 *
sgehrd
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
shst01
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01
Definition: shst01.f:136
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
sorghr
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128
slange
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:116
shseqr
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
strsen
subroutine strsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
STRSEN
Definition: strsen.f:316
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: