LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sget38()

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.

Definition at line 90 of file sget38.f.

91*
92* -- LAPACK test routine --
93* -- LAPACK is a software package provided by Univ. of Tennessee, --
94* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95*
96* .. Scalar Arguments ..
97 INTEGER KNT, NIN
98* ..
99* .. Array Arguments ..
100 INTEGER LMAX( 3 ), NINFO( 3 )
101 REAL RMAX( 3 )
102* ..
103*
104* =====================================================================
105*
106* .. Parameters ..
107 REAL ZERO, ONE, TWO
108 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
109 REAL EPSIN
110 parameter( epsin = 5.9605e-8 )
111 INTEGER LDT, LWORK
112 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
113 INTEGER LIWORK
114 parameter( liwork = ldt*ldt )
115* ..
116* .. Local Scalars ..
117 INTEGER I, INFO, ISCL, ITMP, J, KMIN, M, N, NDIM
118 REAL BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
119 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VIMIN, VMAX,
120 $ VMUL, VRMIN
121* ..
122* .. Local Arrays ..
123 LOGICAL SELECT( LDT )
124 INTEGER IPNT( LDT ), ISELEC( LDT ), IWORK( LIWORK )
125 REAL Q( LDT, LDT ), QSAV( LDT, LDT ),
126 $ QTMP( LDT, LDT ), RESULT( 2 ), T( LDT, LDT ),
127 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
128 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), VAL( 3 ),
129 $ WI( LDT ), WITMP( LDT ), WORK( LWORK ),
130 $ WR( LDT ), WRTMP( LDT )
131* ..
132* .. External Functions ..
133 REAL SLAMCH, SLANGE
134 EXTERNAL slamch, slange
135* ..
136* .. External Subroutines ..
137 EXTERNAL scopy, sgehrd, shseqr, shst01, slacpy, sorghr,
138 $ sscal, strsen
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC max, real, sqrt
142* ..
143* .. Executable Statements ..
144*
145 eps = slamch( 'P' )
146 smlnum = slamch( 'S' ) / eps
147 bignum = one / smlnum
148*
149* EPSIN = 2**(-24) = precision to which input data computed
150*
151 eps = max( eps, epsin )
152 rmax( 1 ) = zero
153 rmax( 2 ) = zero
154 rmax( 3 ) = zero
155 lmax( 1 ) = 0
156 lmax( 2 ) = 0
157 lmax( 3 ) = 0
158 knt = 0
159 ninfo( 1 ) = 0
160 ninfo( 2 ) = 0
161 ninfo( 3 ) = 0
162*
163 val( 1 ) = sqrt( smlnum )
164 val( 2 ) = one
165 val( 3 ) = sqrt( sqrt( bignum ) )
166*
167* Read input data until N=0. Assume input eigenvalues are sorted
168* lexicographically (increasing by real part, then decreasing by
169* imaginary part)
170*
171 10 CONTINUE
172 READ( nin, fmt = * )n, ndim
173 IF( n.EQ.0 )
174 $ RETURN
175 READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
176 DO 20 i = 1, n
177 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
178 20 CONTINUE
179 READ( nin, fmt = * )sin, sepin
180*
181 tnrm = slange( 'M', n, n, tmp, ldt, work )
182 DO 160 iscl = 1, 3
183*
184* Scale input matrix
185*
186 knt = knt + 1
187 CALL slacpy( 'F', n, n, tmp, ldt, t, ldt )
188 vmul = val( iscl )
189 DO 30 i = 1, n
190 CALL sscal( n, vmul, t( 1, i ), 1 )
191 30 CONTINUE
192 IF( tnrm.EQ.zero )
193 $ vmul = one
194 CALL slacpy( 'F', n, n, t, ldt, tsav, ldt )
195*
196* Compute Schur form
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 160
204 END IF
205*
206* Generate orthogonal matrix
207*
208 CALL slacpy( 'L', n, n, t, ldt, q, ldt )
209 CALL sorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
210 $ info )
211*
212* Compute Schur form
213*
214 CALL shseqr( 'S', 'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
215 $ lwork, info )
216 IF( info.NE.0 ) THEN
217 lmax( 2 ) = knt
218 ninfo( 2 ) = ninfo( 2 ) + 1
219 GO TO 160
220 END IF
221*
222* Sort, select eigenvalues
223*
224 DO 40 i = 1, n
225 ipnt( i ) = i
226 SELECT( i ) = .false.
227 40 CONTINUE
228 CALL scopy( n, wr, 1, wrtmp, 1 )
229 CALL scopy( n, wi, 1, witmp, 1 )
230 DO 60 i = 1, n - 1
231 kmin = i
232 vrmin = wrtmp( i )
233 vimin = witmp( i )
234 DO 50 j = i + 1, n
235 IF( wrtmp( j ).LT.vrmin ) THEN
236 kmin = j
237 vrmin = wrtmp( j )
238 vimin = witmp( j )
239 END IF
240 50 CONTINUE
241 wrtmp( kmin ) = wrtmp( i )
242 witmp( kmin ) = witmp( i )
243 wrtmp( i ) = vrmin
244 witmp( i ) = vimin
245 itmp = ipnt( i )
246 ipnt( i ) = ipnt( kmin )
247 ipnt( kmin ) = itmp
248 60 CONTINUE
249 DO 70 i = 1, ndim
250 SELECT( ipnt( iselec( i ) ) ) = .true.
251 70 CONTINUE
252*
253* Compute condition numbers
254*
255 CALL slacpy( 'F', n, n, q, ldt, qsav, ldt )
256 CALL slacpy( 'F', n, n, t, ldt, tsav1, ldt )
257 CALL strsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wrtmp, witmp,
258 $ m, s, sep, work, lwork, iwork, liwork, info )
259 IF( info.NE.0 ) THEN
260 lmax( 3 ) = knt
261 ninfo( 3 ) = ninfo( 3 ) + 1
262 GO TO 160
263 END IF
264 septmp = sep / vmul
265 stmp = s
266*
267* Compute residuals
268*
269 CALL shst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
270 $ result )
271 vmax = max( result( 1 ), result( 2 ) )
272 IF( vmax.GT.rmax( 1 ) ) THEN
273 rmax( 1 ) = vmax
274 IF( ninfo( 1 ).EQ.0 )
275 $ lmax( 1 ) = knt
276 END IF
277*
278* Compare condition number for eigenvalue cluster
279* taking its condition number into account
280*
281 v = max( two*real( n )*eps*tnrm, smlnum )
282 IF( tnrm.EQ.zero )
283 $ v = one
284 IF( v.GT.septmp ) THEN
285 tol = one
286 ELSE
287 tol = v / septmp
288 END IF
289 IF( v.GT.sepin ) THEN
290 tolin = one
291 ELSE
292 tolin = v / sepin
293 END IF
294 tol = max( tol, smlnum / eps )
295 tolin = max( tolin, smlnum / eps )
296 IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
297 vmax = one / eps
298 ELSE IF( sin-tolin.GT.stmp+tol ) THEN
299 vmax = ( sin-tolin ) / ( stmp+tol )
300 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
301 vmax = one / eps
302 ELSE IF( sin+tolin.LT.stmp-tol ) THEN
303 vmax = ( stmp-tol ) / ( sin+tolin )
304 ELSE
305 vmax = one
306 END IF
307 IF( vmax.GT.rmax( 2 ) ) THEN
308 rmax( 2 ) = vmax
309 IF( ninfo( 2 ).EQ.0 )
310 $ lmax( 2 ) = knt
311 END IF
312*
313* Compare condition numbers for invariant subspace
314* taking its condition number into account
315*
316 IF( v.GT.septmp*stmp ) THEN
317 tol = septmp
318 ELSE
319 tol = v / stmp
320 END IF
321 IF( v.GT.sepin*sin ) THEN
322 tolin = sepin
323 ELSE
324 tolin = v / sin
325 END IF
326 tol = max( tol, smlnum / eps )
327 tolin = max( tolin, smlnum / eps )
328 IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
329 vmax = one / eps
330 ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
331 vmax = ( sepin-tolin ) / ( septmp+tol )
332 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
333 vmax = one / eps
334 ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
335 vmax = ( septmp-tol ) / ( sepin+tolin )
336 ELSE
337 vmax = one
338 END IF
339 IF( vmax.GT.rmax( 2 ) ) THEN
340 rmax( 2 ) = vmax
341 IF( ninfo( 2 ).EQ.0 )
342 $ lmax( 2 ) = knt
343 END IF
344*
345* Compare condition number for eigenvalue cluster
346* without taking its condition number into account
347*
348 IF( sin.LE.real( 2*n )*eps .AND. stmp.LE.real( 2*n )*eps ) THEN
349 vmax = one
350 ELSE IF( eps*sin.GT.stmp ) THEN
351 vmax = one / eps
352 ELSE IF( sin.GT.stmp ) THEN
353 vmax = sin / stmp
354 ELSE IF( sin.LT.eps*stmp ) THEN
355 vmax = one / eps
356 ELSE IF( sin.LT.stmp ) THEN
357 vmax = stmp / sin
358 ELSE
359 vmax = one
360 END IF
361 IF( vmax.GT.rmax( 3 ) ) THEN
362 rmax( 3 ) = vmax
363 IF( ninfo( 3 ).EQ.0 )
364 $ lmax( 3 ) = knt
365 END IF
366*
367* Compare condition numbers for invariant subspace
368* without taking its condition number into account
369*
370 IF( sepin.LE.v .AND. septmp.LE.v ) THEN
371 vmax = one
372 ELSE IF( eps*sepin.GT.septmp ) THEN
373 vmax = one / eps
374 ELSE IF( sepin.GT.septmp ) THEN
375 vmax = sepin / septmp
376 ELSE IF( sepin.LT.eps*septmp ) THEN
377 vmax = one / eps
378 ELSE IF( sepin.LT.septmp ) THEN
379 vmax = septmp / sepin
380 ELSE
381 vmax = one
382 END IF
383 IF( vmax.GT.rmax( 3 ) ) THEN
384 rmax( 3 ) = vmax
385 IF( ninfo( 3 ).EQ.0 )
386 $ lmax( 3 ) = knt
387 END IF
388*
389* Compute eigenvalue condition number only and compare
390* Update Q
391*
392 vmax = zero
393 CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
394 CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
395 septmp = -one
396 stmp = -one
397 CALL strsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
398 $ witmp, m, stmp, septmp, work, lwork, iwork,
399 $ liwork, info )
400 IF( info.NE.0 ) THEN
401 lmax( 3 ) = knt
402 ninfo( 3 ) = ninfo( 3 ) + 1
403 GO TO 160
404 END IF
405 IF( s.NE.stmp )
406 $ vmax = one / eps
407 IF( -one.NE.septmp )
408 $ vmax = one / eps
409 DO 90 i = 1, n
410 DO 80 j = 1, n
411 IF( ttmp( i, j ).NE.t( i, j ) )
412 $ vmax = one / eps
413 IF( qtmp( i, j ).NE.q( i, j ) )
414 $ vmax = one / eps
415 80 CONTINUE
416 90 CONTINUE
417*
418* Compute invariant subspace condition number only and compare
419* Update Q
420*
421 CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
422 CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
423 septmp = -one
424 stmp = -one
425 CALL strsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
426 $ witmp, m, stmp, septmp, work, lwork, iwork,
427 $ liwork, info )
428 IF( info.NE.0 ) THEN
429 lmax( 3 ) = knt
430 ninfo( 3 ) = ninfo( 3 ) + 1
431 GO TO 160
432 END IF
433 IF( -one.NE.stmp )
434 $ vmax = one / eps
435 IF( sep.NE.septmp )
436 $ vmax = one / eps
437 DO 110 i = 1, n
438 DO 100 j = 1, n
439 IF( ttmp( i, j ).NE.t( i, j ) )
440 $ vmax = one / eps
441 IF( qtmp( i, j ).NE.q( i, j ) )
442 $ vmax = one / eps
443 100 CONTINUE
444 110 CONTINUE
445*
446* Compute eigenvalue condition number only and compare
447* Do not update Q
448*
449 CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
450 CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
451 septmp = -one
452 stmp = -one
453 CALL strsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
454 $ witmp, m, stmp, septmp, work, lwork, iwork,
455 $ liwork, info )
456 IF( info.NE.0 ) THEN
457 lmax( 3 ) = knt
458 ninfo( 3 ) = ninfo( 3 ) + 1
459 GO TO 160
460 END IF
461 IF( s.NE.stmp )
462 $ vmax = one / eps
463 IF( -one.NE.septmp )
464 $ vmax = one / eps
465 DO 130 i = 1, n
466 DO 120 j = 1, n
467 IF( ttmp( i, j ).NE.t( i, j ) )
468 $ vmax = one / eps
469 IF( qtmp( i, j ).NE.qsav( i, j ) )
470 $ vmax = one / eps
471 120 CONTINUE
472 130 CONTINUE
473*
474* Compute invariant subspace condition number only and compare
475* Do not update Q
476*
477 CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
478 CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
479 septmp = -one
480 stmp = -one
481 CALL strsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
482 $ witmp, m, stmp, septmp, work, lwork, iwork,
483 $ liwork, info )
484 IF( info.NE.0 ) THEN
485 lmax( 3 ) = knt
486 ninfo( 3 ) = ninfo( 3 ) + 1
487 GO TO 160
488 END IF
489 IF( -one.NE.stmp )
490 $ vmax = one / eps
491 IF( sep.NE.septmp )
492 $ vmax = one / eps
493 DO 150 i = 1, n
494 DO 140 j = 1, n
495 IF( ttmp( i, j ).NE.t( i, j ) )
496 $ vmax = one / eps
497 IF( qtmp( i, j ).NE.qsav( i, j ) )
498 $ vmax = one / eps
499 140 CONTINUE
500 150 CONTINUE
501 IF( vmax.GT.rmax( 1 ) ) THEN
502 rmax( 1 ) = vmax
503 IF( ninfo( 1 ).EQ.0 )
504 $ lmax( 1 ) = knt
505 END IF
506 160 CONTINUE
507 GO TO 10
508*
509* End of SGET38
510*
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgehrd
subroutine sgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
SGEHRD
Definition sgehrd.f:167
shseqr
subroutine shseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
SHSEQR
Definition shseqr.f:316
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:103
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:114
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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:314
sorghr
subroutine sorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
SORGHR
Definition sorghr.f:126
shst01
subroutine shst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
SHST01
Definition shst01.f:134
Here is the call graph for this function:
Here is the caller graph for this function: