89 SUBROUTINE dget37( RMAX, LMAX, NINFO, KNT, NIN )
99 INTEGER LMAX( 3 ), NINFO( 3 )
100 DOUBLE PRECISION RMAX( 3 )
106 DOUBLE PRECISION ZERO, ONE, TWO
107 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
108 DOUBLE PRECISION EPSIN
109 parameter( epsin = 5.9605d-8 )
111 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
114 INTEGER I, ICMP, IFND, INFO, ISCL, J, KMIN, M, N
115 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116 $ VIMIN, VMAX, VMUL, VRMIN
119 LOGICAL SELECT( LDT )
120 INTEGER IWORK( 2*LDT ), LCMP( 3 )
121 DOUBLE PRECISION DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
122 $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
123 $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
124 $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
125 $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
126 $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
130 DOUBLE PRECISION DLAMCH, DLANGE
131 EXTERNAL dlamch, dlange
138 INTRINSIC dble, max, sqrt
143 smlnum = dlamch(
'S' ) / eps
144 bignum = one / smlnum
148 eps = max( eps, epsin )
160 val( 1 ) = sqrt( smlnum )
162 val( 3 ) = sqrt( bignum )
169 READ( nin, fmt = * )n
173 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
176 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
178 tnrm = dlange(
'M', n, n, tmp, ldt, work )
187 CALL dlacpy(
'F', n, n, tmp, ldt, t, ldt )
190 CALL dscal( n, vmul, t( 1, i ), 1 )
197 CALL dgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
201 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL dhseqr(
'S',
'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
216 ninfo( 2 ) = ninfo( 2 ) + 1
222 CALL dtrevc(
'Both',
'All',
SELECT, n, t, ldt, le, ldt, re,
223 $ ldt, n, m, work, info )
227 CALL dtrsna(
'Both',
'All',
SELECT, n, t, ldt, le, ldt, re,
228 $ ldt, s, sep, n, m, work, n, iwork, info )
231 ninfo( 3 ) = ninfo( 3 ) + 1
238 CALL dcopy( n, wr, 1, wrtmp, 1 )
239 CALL dcopy( n, wi, 1, witmp, 1 )
240 CALL dcopy( n, s, 1, stmp, 1 )
241 CALL dcopy( n, sep, 1, septmp, 1 )
242 CALL dscal( n, one / vmul, septmp, 1 )
248 IF( wrtmp( j ).LT.vrmin )
THEN
254 wrtmp( kmin ) = wrtmp( i )
255 witmp( kmin ) = witmp( i )
259 stmp( kmin ) = stmp( i )
261 vrmin = septmp( kmin )
262 septmp( kmin ) = septmp( i )
269 v = max( two*dble( n )*eps*tnrm, smlnum )
273 IF( v.GT.septmp( i ) )
THEN
276 tol = v / septmp( i )
278 IF( v.GT.sepin( i ) )
THEN
281 tolin = v / sepin( i )
283 tol = max( tol, smlnum / eps )
284 tolin = max( tolin, smlnum / eps )
285 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol )
THEN
287 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol )
THEN
288 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
289 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) )
THEN
291 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol )
THEN
292 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
296 IF( vmax.GT.rmax( 2 ) )
THEN
298 IF( ninfo( 2 ).EQ.0 )
307 IF( v.GT.septmp( i )*stmp( i ) )
THEN
312 IF( v.GT.sepin( i )*sin( i ) )
THEN
317 tol = max( tol, smlnum / eps )
318 tolin = max( tolin, smlnum / eps )
319 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol )
THEN
321 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol )
THEN
322 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
323 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) )
THEN
325 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol )
THEN
326 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
330 IF( vmax.GT.rmax( 2 ) )
THEN
332 IF( ninfo( 2 ).EQ.0 )
341 IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
342 $ dble( 2*n )*eps )
THEN
344 ELSE IF( eps*sin( i ).GT.stmp( i ) )
THEN
346 ELSE IF( sin( i ).GT.stmp( i ) )
THEN
347 vmax = sin( i ) / stmp( i )
348 ELSE IF( sin( i ).LT.eps*stmp( i ) )
THEN
350 ELSE IF( sin( i ).LT.stmp( i ) )
THEN
351 vmax = stmp( i ) / sin( i )
355 IF( vmax.GT.rmax( 3 ) )
THEN
357 IF( ninfo( 3 ).EQ.0 )
366 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v )
THEN
368 ELSE IF( eps*sepin( i ).GT.septmp( i ) )
THEN
370 ELSE IF( sepin( i ).GT.septmp( i ) )
THEN
371 vmax = sepin( i ) / septmp( i )
372 ELSE IF( sepin( i ).LT.eps*septmp( i ) )
THEN
374 ELSE IF( sepin( i ).LT.septmp( i ) )
THEN
375 vmax = septmp( i ) / sepin( i )
379 IF( vmax.GT.rmax( 3 ) )
THEN
381 IF( ninfo( 3 ).EQ.0 )
390 CALL dcopy( n, dum, 0, stmp, 1 )
391 CALL dcopy( n, dum, 0, septmp, 1 )
392 CALL dtrsna(
'Eigcond',
'All',
SELECT, n, t, ldt, le, ldt, re,
393 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
396 ninfo( 3 ) = ninfo( 3 ) + 1
400 IF( stmp( i ).NE.s( i ) )
402 IF( septmp( i ).NE.dum( 1 ) )
408 CALL dcopy( n, dum, 0, stmp, 1 )
409 CALL dcopy( n, dum, 0, septmp, 1 )
410 CALL dtrsna(
'Veccond',
'All',
SELECT, n, t, ldt, le, ldt, re,
411 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
414 ninfo( 3 ) = ninfo( 3 ) + 1
418 IF( stmp( i ).NE.dum( 1 ) )
420 IF( septmp( i ).NE.sep( i ) )
429 CALL dcopy( n, dum, 0, stmp, 1 )
430 CALL dcopy( n, dum, 0, septmp, 1 )
431 CALL dtrsna(
'Bothcond',
'Some',
SELECT, n, t, ldt, le, ldt,
432 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
436 ninfo( 3 ) = ninfo( 3 ) + 1
440 IF( septmp( i ).NE.sep( i ) )
442 IF( stmp( i ).NE.s( i ) )
448 CALL dcopy( n, dum, 0, stmp, 1 )
449 CALL dcopy( n, dum, 0, septmp, 1 )
450 CALL dtrsna(
'Eigcond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
451 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
454 ninfo( 3 ) = ninfo( 3 ) + 1
458 IF( stmp( i ).NE.s( i ) )
460 IF( septmp( i ).NE.dum( 1 ) )
466 CALL dcopy( n, dum, 0, stmp, 1 )
467 CALL dcopy( n, dum, 0, septmp, 1 )
468 CALL dtrsna(
'Veccond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
469 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
472 ninfo( 3 ) = ninfo( 3 ) + 1
476 IF( stmp( i ).NE.dum( 1 ) )
478 IF( septmp( i ).NE.sep( i ) )
481 IF( vmax.GT.rmax( 1 ) )
THEN
483 IF( ninfo( 1 ).EQ.0 )
489 IF( wi( 1 ).EQ.zero )
THEN
493 IF( ifnd.EQ.1 .OR. wi( i ).EQ.zero )
THEN
494 SELECT( i ) = .false.
499 CALL dcopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
500 CALL dcopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
501 CALL dcopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
502 CALL dcopy( n, le( 1, i+1 ), 1, le( 1, 3 ), 1 )
515 IF( ifnd.EQ.1 .OR. wi( i ).NE.zero )
THEN
516 SELECT( i ) = .false.
520 CALL dcopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
521 CALL dcopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
533 CALL dcopy( icmp, dum, 0, stmp, 1 )
534 CALL dcopy( icmp, dum, 0, septmp, 1 )
535 CALL dtrsna(
'Bothcond',
'Some',
SELECT, n, t, ldt, le, ldt,
536 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
540 ninfo( 3 ) = ninfo( 3 ) + 1
545 IF( septmp( i ).NE.sep( j ) )
547 IF( stmp( i ).NE.s( j ) )
553 CALL dcopy( icmp, dum, 0, stmp, 1 )
554 CALL dcopy( icmp, dum, 0, septmp, 1 )
555 CALL dtrsna(
'Eigcond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
556 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
559 ninfo( 3 ) = ninfo( 3 ) + 1
564 IF( stmp( i ).NE.s( j ) )
566 IF( septmp( i ).NE.dum( 1 ) )
572 CALL dcopy( icmp, dum, 0, stmp, 1 )
573 CALL dcopy( icmp, dum, 0, septmp, 1 )
574 CALL dtrsna(
'Veccond',
'Some',
SELECT, n, t, ldt, le, ldt, re,
575 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
578 ninfo( 3 ) = ninfo( 3 ) + 1
583 IF( stmp( i ).NE.dum( 1 ) )
585 IF( septmp( i ).NE.sep( j ) )
588 IF( vmax.GT.rmax( 1 ) )
THEN
590 IF( ninfo( 1 ).EQ.0 )