91 SUBROUTINE sget37( RMAX, LMAX, NINFO, KNT, NIN )
102 INTEGER LMAX( 3 ), NINFO( 3 )
110 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
112 parameter( epsin = 5.9605e-8 )
114 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
117 INTEGER I, ICMP, IFND, INFO, ISCL, J, KMIN, M, N
118 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
119 $ vimin, vmax, vmul, vrmin
122 LOGICAL SELECT( ldt )
123 INTEGER IWORK( 2*ldt ), LCMP( 3 )
124 REAL DUM( 1 ), LE( ldt, ldt ), RE( ldt, ldt ),
125 $ s( ldt ), sep( ldt ), sepin( ldt ),
126 $ septmp( ldt ), sin( ldt ), stmp( ldt ),
127 $ t( ldt, ldt ), tmp( ldt, ldt ), val( 3 ),
128 $ wi( ldt ), wiin( ldt ), witmp( ldt ),
129 $ work( lwork ), wr( ldt ), wrin( 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( bignum )
173 READ( nin, fmt = * )n
177 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
180 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
182 tnrm = slange(
'M', n, n, tmp, ldt, work )
191 CALL
slacpy(
'F', n, n, tmp, ldt, t, ldt )
194 CALL
sscal( n, vmul, t( 1, i ), 1 )
201 CALL
sgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
205 ninfo( 1 ) = ninfo( 1 ) + 1
216 CALL
shseqr(
'S',
'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
220 ninfo( 2 ) = ninfo( 2 ) + 1
226 CALL
strevc(
'Both',
'All',
SELECT, n, t, ldt, le, ldt, re,
227 $ ldt, n, m, work, info )
231 CALL
strsna(
'Both',
'All',
SELECT, n, t, ldt, le, ldt, re,
232 $ ldt, s, sep, n, m, work, n, iwork, info )
235 ninfo( 3 ) = ninfo( 3 ) + 1
242 CALL
scopy( n, wr, 1, wrtmp, 1 )
243 CALL
scopy( n, wi, 1, witmp, 1 )
244 CALL
scopy( n, s, 1, stmp, 1 )
245 CALL
scopy( n, sep, 1, septmp, 1 )
246 CALL
sscal( n, one / vmul, septmp, 1 )
252 IF( wrtmp( j ).LT.vrmin )
THEN
258 wrtmp( kmin ) = wrtmp( i )
259 witmp( kmin ) = witmp( i )
263 stmp( kmin ) = stmp( i )
265 vrmin = septmp( kmin )
266 septmp( kmin ) = septmp( i )
273 v = max( two*
REAL( n )*EPS*TNRM, SMLNUM )
277 IF( v.GT.septmp( i ) )
THEN
280 tol = v / septmp( i )
282 IF( v.GT.sepin( i ) )
THEN
285 tolin = v / sepin( i )
287 tol = max( tol, smlnum / eps )
288 tolin = max( tolin, smlnum / eps )
289 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol )
THEN
291 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol )
THEN
292 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
293 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) )
THEN
295 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol )
THEN
296 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
300 IF( vmax.GT.rmax( 2 ) )
THEN
302 IF( ninfo( 2 ).EQ.0 )
311 IF( v.GT.septmp( i )*stmp( i ) )
THEN
316 IF( v.GT.sepin( i )*sin( i ) )
THEN
321 tol = max( tol, smlnum / eps )
322 tolin = max( tolin, smlnum / eps )
323 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol )
THEN
325 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol )
THEN
326 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
327 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) )
THEN
329 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol )
THEN
330 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
334 IF( vmax.GT.rmax( 2 ) )
THEN
336 IF( ninfo( 2 ).EQ.0 )
345 IF( sin( i ).LE.
REAL( 2*n )*EPS .AND. STMP( i ).LE.
346 $
REAL( 2*n )*EPS ) then
348 ELSE IF( eps*sin( i ).GT.stmp( i ) )
THEN
350 ELSE IF( sin( i ).GT.stmp( i ) )
THEN
351 vmax = sin( i ) / stmp( i )
352 ELSE IF( sin( i ).LT.eps*stmp( i ) )
THEN
354 ELSE IF( sin( i ).LT.stmp( i ) )
THEN
355 vmax = stmp( i ) / sin( i )
359 IF( vmax.GT.rmax( 3 ) )
THEN
361 IF( ninfo( 3 ).EQ.0 )
370 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v )
THEN
372 ELSE IF( eps*sepin( i ).GT.septmp( i ) )
THEN
374 ELSE IF( sepin( i ).GT.septmp( i ) )
THEN
375 vmax = sepin( i ) / septmp( i )
376 ELSE IF( sepin( i ).LT.eps*septmp( i ) )
THEN
378 ELSE IF( sepin( i ).LT.septmp( i ) )
THEN
379 vmax = septmp( i ) / sepin( i )
383 IF( vmax.GT.rmax( 3 ) )
THEN
385 IF( ninfo( 3 ).EQ.0 )
394 CALL
scopy( n, dum, 0, stmp, 1 )
395 CALL
scopy( n, dum, 0, septmp, 1 )
396 CALL
strsna(
'Eigcond',
'All',
SELECT, n, t, ldt, le, ldt, re,
397 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
400 ninfo( 3 ) = ninfo( 3 ) + 1
404 IF( stmp( i ).NE.s( i ) )
406 IF( septmp( i ).NE.dum( 1 ) )
412 CALL
scopy( n, dum, 0, stmp, 1 )
413 CALL
scopy( n, dum, 0, septmp, 1 )
414 CALL
strsna(
'Veccond',
'All',
SELECT, n, t, ldt, le, ldt, re,
415 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
418 ninfo( 3 ) = ninfo( 3 ) + 1
422 IF( stmp( i ).NE.dum( 1 ) )
424 IF( septmp( i ).NE.sep( i ) )
433 CALL
scopy( n, dum, 0, stmp, 1 )
434 CALL
scopy( n, dum, 0, septmp, 1 )
435 CALL
strsna(
'Bothcond',
'Some',
SELECT, n, t, ldt, le, ldt,
436 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
440 ninfo( 3 ) = ninfo( 3 ) + 1
444 IF( septmp( i ).NE.sep( i ) )
446 IF( stmp( i ).NE.s( i ) )
452 CALL
scopy( n, dum, 0, stmp, 1 )
453 CALL
scopy( n, dum, 0, septmp, 1 )
454 CALL
strsna(
'Eigcond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
455 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
458 ninfo( 3 ) = ninfo( 3 ) + 1
462 IF( stmp( i ).NE.s( i ) )
464 IF( septmp( i ).NE.dum( 1 ) )
470 CALL
scopy( n, dum, 0, stmp, 1 )
471 CALL
scopy( n, dum, 0, septmp, 1 )
472 CALL
strsna(
'Veccond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
473 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
476 ninfo( 3 ) = ninfo( 3 ) + 1
480 IF( stmp( i ).NE.dum( 1 ) )
482 IF( septmp( i ).NE.sep( i ) )
485 IF( vmax.GT.rmax( 1 ) )
THEN
487 IF( ninfo( 1 ).EQ.0 )
493 IF( wi( 1 ).EQ.zero )
THEN
497 IF( ifnd.EQ.1 .OR. wi( i ).EQ.zero )
THEN
498 SELECT( i ) = .false.
503 CALL
scopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
504 CALL
scopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
505 CALL
scopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
506 CALL
scopy( n, le( 1, i+1 ), 1, le( 1, 3 ), 1 )
519 IF( ifnd.EQ.1 .OR. wi( i ).NE.zero )
THEN
520 SELECT( i ) = .false.
524 CALL
scopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
525 CALL
scopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
537 CALL
scopy( icmp, dum, 0, stmp, 1 )
538 CALL
scopy( icmp, dum, 0, septmp, 1 )
539 CALL
strsna(
'Bothcond',
'Some',
SELECT, n, t, ldt, le, ldt,
540 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
544 ninfo( 3 ) = ninfo( 3 ) + 1
549 IF( septmp( i ).NE.sep( j ) )
551 IF( stmp( i ).NE.s( j ) )
557 CALL
scopy( icmp, dum, 0, stmp, 1 )
558 CALL
scopy( icmp, dum, 0, septmp, 1 )
559 CALL
strsna(
'Eigcond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
560 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
563 ninfo( 3 ) = ninfo( 3 ) + 1
568 IF( stmp( i ).NE.s( j ) )
570 IF( septmp( i ).NE.dum( 1 ) )
576 CALL
scopy( icmp, dum, 0, stmp, 1 )
577 CALL
scopy( icmp, dum, 0, septmp, 1 )
578 CALL
strsna(
'Veccond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
579 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
582 ninfo( 3 ) = ninfo( 3 ) + 1
587 IF( stmp( i ).NE.dum( 1 ) )
589 IF( septmp( i ).NE.sep( j ) )
592 IF( vmax.GT.rmax( 1 ) )
THEN
594 IF( ninfo( 1 ).EQ.0 )
subroutine strsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
STRSNA
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine strevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STREVC
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
subroutine sget37(RMAX, LMAX, NINFO, KNT, NIN)
SGET37
subroutine sscal(N, SA, SX, INCX)
SSCAL