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 )
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
subroutine zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZGEHRD
subroutine zhseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
ZHSEQR
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine ztrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
ZTREVC
subroutine ztrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, rwork, info)
ZTRSNA
subroutine zget37(rmax, lmax, ninfo, knt, nin)
ZGET37