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.```
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 *
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01
Definition: shst01.f:136
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
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128
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
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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: