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 )
subroutine dget37(rmax, lmax, ninfo, knt, nin)
DGET37
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine dtrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, info)
DTREVC
subroutine dtrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, iwork, info)
DTRSNA