91 SUBROUTINE dget37( RMAX, LMAX, NINFO, KNT, NIN )
102 INTEGER LMAX( 3 ), NINFO( 3 )
103 DOUBLE PRECISION RMAX( 3 )
109 DOUBLE PRECISION ZERO, ONE, TWO
110 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
111 DOUBLE PRECISION EPSIN
112 parameter ( epsin = 5.9605d-8 )
114 parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
117 INTEGER I, ICMP, IFND, INFO, ISCL, J, KMIN, M, N
118 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
119 $ vimin, vmax, vmul, vrmin
122 LOGICAL SELECT( ldt )
123 INTEGER IWORK( 2*ldt ), LCMP( 3 )
124 DOUBLE PRECISION 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 ),
133 DOUBLE PRECISION DLAMCH, DLANGE
134 EXTERNAL dlamch, dlange
141 INTRINSIC dble, max, sqrt
146 smlnum = dlamch(
'S' ) / eps
147 bignum = one / smlnum
148 CALL dlabad( 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 = dlange(
'M', n, n, tmp, ldt, work )
191 CALL dlacpy(
'F', n, n, tmp, ldt, t, ldt )
194 CALL dscal( n, vmul, t( 1, i ), 1 )
201 CALL dgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
205 ninfo( 1 ) = ninfo( 1 ) + 1
216 CALL dhseqr(
'S',
'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
220 ninfo( 2 ) = ninfo( 2 ) + 1
226 CALL dtrevc(
'Both',
'All',
SELECT, n, t, ldt, le, ldt, re,
227 $ ldt, n, m, work, info )
231 CALL dtrsna(
'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 dcopy( n, wr, 1, wrtmp, 1 )
243 CALL dcopy( n, wi, 1, witmp, 1 )
244 CALL dcopy( n, s, 1, stmp, 1 )
245 CALL dcopy( n, sep, 1, septmp, 1 )
246 CALL dscal( 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*dble( 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.dble( 2*n )*eps .AND. stmp( i ).LE.
346 $ dble( 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 dcopy( n, dum, 0, stmp, 1 )
395 CALL dcopy( n, dum, 0, septmp, 1 )
396 CALL dtrsna(
'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 dcopy( n, dum, 0, stmp, 1 )
413 CALL dcopy( n, dum, 0, septmp, 1 )
414 CALL dtrsna(
'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 dcopy( n, dum, 0, stmp, 1 )
434 CALL dcopy( n, dum, 0, septmp, 1 )
435 CALL dtrsna(
'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 dcopy( n, dum, 0, stmp, 1 )
453 CALL dcopy( n, dum, 0, septmp, 1 )
454 CALL dtrsna(
'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 dcopy( n, dum, 0, stmp, 1 )
471 CALL dcopy( n, dum, 0, septmp, 1 )
472 CALL dtrsna(
'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 dcopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
504 CALL dcopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
505 CALL dcopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
506 CALL dcopy( 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 dcopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
525 CALL dcopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
537 CALL dcopy( icmp, dum, 0, stmp, 1 )
538 CALL dcopy( icmp, dum, 0, septmp, 1 )
539 CALL dtrsna(
'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 dcopy( icmp, dum, 0, stmp, 1 )
558 CALL dcopy( icmp, dum, 0, septmp, 1 )
559 CALL dtrsna(
'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 dcopy( icmp, dum, 0, stmp, 1 )
577 CALL dcopy( icmp, dum, 0, septmp, 1 )
578 CALL dtrsna(
'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 dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
subroutine dlabad(SMALL, LARGE)
DLABAD
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
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
subroutine dget37(RMAX, LMAX, NINFO, KNT, NIN)
DGET37