90 SUBROUTINE sget38( RMAX, LMAX, NINFO, KNT, NIN )
100 INTEGER LMAX( 3 ), NINFO( 3 )
108 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
110 parameter( epsin = 5.9605e-8 )
112 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
114 parameter( liwork = ldt*ldt )
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,
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 )
134 EXTERNAL slamch, slange
141 INTRINSIC max, real, sqrt
146 smlnum = slamch(
'S' ) / eps
147 bignum = one / smlnum
148 CALL slabad( smlnum, bignum )
152 eps = max( eps, epsin )
164 val( 1 ) = sqrt( smlnum )
166 val( 3 ) = sqrt( sqrt( bignum ) )
173 READ( nin, fmt = * )n, ndim
176 READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
178 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
180 READ( nin, fmt = * )sin, sepin
182 tnrm = slange(
'M', n, n, tmp, ldt, work )
188 CALL slacpy(
'F', n, n, tmp, ldt, t, ldt )
191 CALL sscal( n, vmul, t( 1, i ), 1 )
195 CALL slacpy(
'F', n, n, t, ldt, tsav, ldt )
199 CALL sgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
203 ninfo( 1 ) = ninfo( 1 ) + 1
209 CALL slacpy(
'L', n, n, t, ldt, q, ldt )
210 CALL sorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
215 CALL shseqr(
'S',
'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
219 ninfo( 2 ) = ninfo( 2 ) + 1
227 SELECT( i ) = .false.
229 CALL scopy( n, wr, 1, wrtmp, 1 )
230 CALL scopy( n, wi, 1, witmp, 1 )
236 IF( wrtmp( j ).LT.vrmin )
THEN
242 wrtmp( kmin ) = wrtmp( i )
243 witmp( kmin ) = witmp( i )
247 ipnt( i ) = ipnt( kmin )
251 SELECT( ipnt( iselec( i ) ) ) = .true.
256 CALL slacpy(
'F', n, n, q, ldt, qsav, ldt )
257 CALL slacpy(
'F', n, n, t, ldt, tsav1, ldt )
258 CALL strsen(
'B',
'V',
SELECT, n, t, ldt, q, ldt, wrtmp, witmp,
259 $ m, s, sep, work, lwork, iwork, liwork, info )
262 ninfo( 3 ) = ninfo( 3 ) + 1
270 CALL shst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
272 vmax = max( result( 1 ), result( 2 ) )
273 IF( vmax.GT.rmax( 1 ) )
THEN
275 IF( ninfo( 1 ).EQ.0 )
282 v = max( two*real( n )*eps*tnrm, smlnum )
285 IF( v.GT.septmp )
THEN
290 IF( v.GT.sepin )
THEN
295 tol = max( tol, smlnum / eps )
296 tolin = max( tolin, smlnum / eps )
297 IF( eps*( sin-tolin ).GT.stmp+tol )
THEN
299 ELSE IF( sin-tolin.GT.stmp+tol )
THEN
300 vmax = ( sin-tolin ) / ( stmp+tol )
301 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) )
THEN
303 ELSE IF( sin+tolin.LT.stmp-tol )
THEN
304 vmax = ( stmp-tol ) / ( sin+tolin )
308 IF( vmax.GT.rmax( 2 ) )
THEN
310 IF( ninfo( 2 ).EQ.0 )
317 IF( v.GT.septmp*stmp )
THEN
322 IF( v.GT.sepin*sin )
THEN
327 tol = max( tol, smlnum / eps )
328 tolin = max( tolin, smlnum / eps )
329 IF( eps*( sepin-tolin ).GT.septmp+tol )
THEN
331 ELSE IF( sepin-tolin.GT.septmp+tol )
THEN
332 vmax = ( sepin-tolin ) / ( septmp+tol )
333 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) )
THEN
335 ELSE IF( sepin+tolin.LT.septmp-tol )
THEN
336 vmax = ( septmp-tol ) / ( sepin+tolin )
340 IF( vmax.GT.rmax( 2 ) )
THEN
342 IF( ninfo( 2 ).EQ.0 )
349 IF( sin.LE.real( 2*n )*eps .AND. stmp.LE.real( 2*n )*eps )
THEN
351 ELSE IF( eps*sin.GT.stmp )
THEN
353 ELSE IF( sin.GT.stmp )
THEN
355 ELSE IF( sin.LT.eps*stmp )
THEN
357 ELSE IF( sin.LT.stmp )
THEN
362 IF( vmax.GT.rmax( 3 ) )
THEN
364 IF( ninfo( 3 ).EQ.0 )
371 IF( sepin.LE.v .AND. septmp.LE.v )
THEN
373 ELSE IF( eps*sepin.GT.septmp )
THEN
375 ELSE IF( sepin.GT.septmp )
THEN
376 vmax = sepin / septmp
377 ELSE IF( sepin.LT.eps*septmp )
THEN
379 ELSE IF( sepin.LT.septmp )
THEN
380 vmax = septmp / sepin
384 IF( vmax.GT.rmax( 3 ) )
THEN
386 IF( ninfo( 3 ).EQ.0 )
394 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
395 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
398 CALL strsen(
'E',
'V',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
399 $ witmp, m, stmp, septmp, work, lwork, iwork,
403 ninfo( 3 ) = ninfo( 3 ) + 1
412 IF( ttmp( i, j ).NE.t( i, j ) )
414 IF( qtmp( i, j ).NE.q( i, j ) )
422 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
423 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
426 CALL strsen(
'V',
'V',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
427 $ witmp, m, stmp, septmp, work, lwork, iwork,
431 ninfo( 3 ) = ninfo( 3 ) + 1
440 IF( ttmp( i, j ).NE.t( i, j ) )
442 IF( qtmp( i, j ).NE.q( i, j ) )
450 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
451 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
454 CALL strsen(
'E',
'N',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
455 $ witmp, m, stmp, septmp, work, lwork, iwork,
459 ninfo( 3 ) = ninfo( 3 ) + 1
468 IF( ttmp( i, j ).NE.t( i, j ) )
470 IF( qtmp( i, j ).NE.qsav( i, j ) )
478 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
479 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
482 CALL strsen(
'V',
'N',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
483 $ witmp, m, stmp, septmp, work, lwork, iwork,
487 ninfo( 3 ) = ninfo( 3 ) + 1
496 IF( ttmp( i, j ).NE.t( i, j ) )
498 IF( qtmp( i, j ).NE.qsav( i, j ) )
502 IF( vmax.GT.rmax( 1 ) )
THEN
504 IF( ninfo( 1 ).EQ.0 )
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
subroutine strsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
STRSEN
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine sget38(RMAX, LMAX, NINFO, KNT, NIN)
SGET38
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01