89 SUBROUTINE zget37( 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, INFO, ISCL, ISRT, J, KMIN, M, N
115 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116 $ VCMIN, VMAX, VMIN, VMUL
119 LOGICAL SELECT( LDT )
121 DOUBLE PRECISION DUM( 1 ), RWORK( 2*LDT ), S( LDT ), SEP( LDT ),
122 $ SEPIN( LDT ), SEPTMP( LDT ), SIN( LDT ),
123 $ STMP( LDT ), VAL( 3 ), WIIN( LDT ),
124 $ WRIN( LDT ), WSRT( LDT )
125 COMPLEX*16 CDUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
126 $ T( LDT, LDT ), TMP( LDT, LDT ), W( LDT ),
127 $ WORK( LWORK ), WTMP( LDT )
130 DOUBLE PRECISION DLAMCH, ZLANGE
131 EXTERNAL dlamch, zlange
138 INTRINSIC dble, dimag, max, sqrt
143 smlnum = dlamch(
'S' ) / eps
144 bignum = one / smlnum
148 eps = max( eps, epsin )
159 val( 1 ) = sqrt( smlnum )
161 val( 3 ) = sqrt( bignum )
168 READ( nin, fmt = * )n, isrt
172 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
175 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
177 tnrm = zlange(
'M', n, n, tmp, ldt, rwork )
183 CALL zlacpy(
'F', n, n, tmp, ldt, t, ldt )
186 CALL zdscal( n, vmul, t( 1, i ), 1 )
193 CALL zgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
197 ninfo( 1 ) = ninfo( 1 ) + 1
208 CALL zhseqr(
'S',
'N', n, 1, n, t, ldt, w, cdum, 1, work,
212 ninfo( 2 ) = ninfo( 2 ) + 1
221 CALL ztrevc(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, n,
222 $ m, work, rwork, info )
226 CALL ztrsna(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, s,
227 $ sep, n, m, work, n, rwork, info )
230 ninfo( 3 ) = ninfo( 3 ) + 1
237 CALL zcopy( n, w, 1, wtmp, 1 )
243 wsrt( i ) = dble( w( i ) )
250 wsrt( i ) = dimag( w( i ) )
253 CALL dcopy( n, s, 1, stmp, 1 )
254 CALL dcopy( n, sep, 1, septmp, 1 )
255 CALL dscal( n, one / vmul, septmp, 1 )
260 IF( wsrt( j ).LT.vmin )
THEN
265 wsrt( kmin ) = wsrt( i )
267 vcmin = dble( wtmp( i ) )
268 wtmp( i ) = w( kmin )
271 stmp( kmin ) = stmp( i )
273 vmin = septmp( kmin )
274 septmp( kmin ) = septmp( i )
281 v = max( two*dble( n )*eps*tnrm, smlnum )
285 IF( v.GT.septmp( i ) )
THEN
288 tol = v / septmp( i )
290 IF( v.GT.sepin( i ) )
THEN
293 tolin = v / sepin( i )
295 tol = max( tol, smlnum / eps )
296 tolin = max( tolin, smlnum / eps )
297 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol )
THEN
299 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol )
THEN
300 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
301 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) )
THEN
303 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol )
THEN
304 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
308 IF( vmax.GT.rmax( 2 ) )
THEN
310 IF( ninfo( 2 ).EQ.0 )
319 IF( v.GT.septmp( i )*stmp( i ) )
THEN
324 IF( v.GT.sepin( i )*sin( i ) )
THEN
329 tol = max( tol, smlnum / eps )
330 tolin = max( tolin, smlnum / eps )
331 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol )
THEN
333 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol )
THEN
334 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
335 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) )
THEN
337 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol )
THEN
338 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
342 IF( vmax.GT.rmax( 2 ) )
THEN
344 IF( ninfo( 2 ).EQ.0 )
353 IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
354 $ dble( 2*n )*eps )
THEN
356 ELSE IF( eps*sin( i ).GT.stmp( i ) )
THEN
358 ELSE IF( sin( i ).GT.stmp( i ) )
THEN
359 vmax = sin( i ) / stmp( i )
360 ELSE IF( sin( i ).LT.eps*stmp( i ) )
THEN
362 ELSE IF( sin( i ).LT.stmp( i ) )
THEN
363 vmax = stmp( i ) / sin( i )
367 IF( vmax.GT.rmax( 3 ) )
THEN
369 IF( ninfo( 3 ).EQ.0 )
378 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v )
THEN
380 ELSE IF( eps*sepin( i ).GT.septmp( i ) )
THEN
382 ELSE IF( sepin( i ).GT.septmp( i ) )
THEN
383 vmax = sepin( i ) / septmp( i )
384 ELSE IF( sepin( i ).LT.eps*septmp( i ) )
THEN
386 ELSE IF( sepin( i ).LT.septmp( i ) )
THEN
387 vmax = septmp( i ) / sepin( i )
391 IF( vmax.GT.rmax( 3 ) )
THEN
393 IF( ninfo( 3 ).EQ.0 )
402 CALL dcopy( n, dum, 0, stmp, 1 )
403 CALL dcopy( n, dum, 0, septmp, 1 )
404 CALL ztrsna(
'E',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
405 $ stmp, septmp, n, m, work, n, rwork, info )
408 ninfo( 3 ) = ninfo( 3 ) + 1
412 IF( stmp( i ).NE.s( i ) )
414 IF( septmp( i ).NE.dum( 1 ) )
420 CALL dcopy( n, dum, 0, stmp, 1 )
421 CALL dcopy( n, dum, 0, septmp, 1 )
422 CALL ztrsna(
'V',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
423 $ stmp, septmp, n, m, work, n, rwork, info )
426 ninfo( 3 ) = ninfo( 3 ) + 1
430 IF( stmp( i ).NE.dum( 1 ) )
432 IF( septmp( i ).NE.sep( i ) )
441 CALL dcopy( n, dum, 0, stmp, 1 )
442 CALL dcopy( n, dum, 0, septmp, 1 )
443 CALL ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
444 $ stmp, septmp, n, m, work, n, rwork, info )
447 ninfo( 3 ) = ninfo( 3 ) + 1
451 IF( septmp( i ).NE.sep( i ) )
453 IF( stmp( i ).NE.s( i ) )
459 CALL dcopy( n, dum, 0, stmp, 1 )
460 CALL dcopy( n, dum, 0, septmp, 1 )
461 CALL ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
462 $ stmp, septmp, n, m, work, n, rwork, info )
465 ninfo( 3 ) = ninfo( 3 ) + 1
469 IF( stmp( i ).NE.s( i ) )
471 IF( septmp( i ).NE.dum( 1 ) )
477 CALL dcopy( n, dum, 0, stmp, 1 )
478 CALL dcopy( n, dum, 0, septmp, 1 )
479 CALL ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
480 $ stmp, septmp, n, m, work, n, rwork, info )
483 ninfo( 3 ) = ninfo( 3 ) + 1
487 IF( stmp( i ).NE.dum( 1 ) )
489 IF( septmp( i ).NE.sep( i ) )
492 IF( vmax.GT.rmax( 1 ) )
THEN
494 IF( ninfo( 1 ).EQ.0 )
501 SELECT( i ) = .false.
508 CALL zcopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
509 CALL zcopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
514 SELECT( n-1 ) = .true.
515 CALL zcopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
516 CALL zcopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
521 CALL dcopy( icmp, dum, 0, stmp, 1 )
522 CALL dcopy( icmp, dum, 0, septmp, 1 )
523 CALL ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
524 $ stmp, septmp, n, m, work, n, rwork, info )
527 ninfo( 3 ) = ninfo( 3 ) + 1
532 IF( septmp( i ).NE.sep( j ) )
534 IF( stmp( i ).NE.s( j ) )
540 CALL dcopy( icmp, dum, 0, stmp, 1 )
541 CALL dcopy( icmp, dum, 0, septmp, 1 )
542 CALL ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
543 $ stmp, septmp, n, m, work, n, rwork, info )
546 ninfo( 3 ) = ninfo( 3 ) + 1
551 IF( stmp( i ).NE.s( j ) )
553 IF( septmp( i ).NE.dum( 1 ) )
559 CALL dcopy( icmp, dum, 0, stmp, 1 )
560 CALL dcopy( icmp, dum, 0, septmp, 1 )
561 CALL ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
562 $ stmp, septmp, n, m, work, n, rwork, info )
565 ninfo( 3 ) = ninfo( 3 ) + 1
570 IF( stmp( i ).NE.dum( 1 ) )
572 IF( septmp( i ).NE.sep( j ) )
575 IF( vmax.GT.rmax( 1 ) )
THEN
577 IF( ninfo( 1 ).EQ.0 )