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
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 )