91 SUBROUTINE zget37( 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, INFO, ISCL, ISRT, J, KMIN, M, N
118 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
119 $ vcmin, vmax, vmin, vmul
122 LOGICAL SELECT( ldt )
124 DOUBLE PRECISION DUM( 1 ), RWORK( 2*ldt ), S( ldt ), SEP( ldt ),
125 $ sepin( ldt ), septmp( ldt ), sin( ldt ),
126 $ stmp( ldt ), val( 3 ), wiin( ldt ),
127 $ wrin( ldt ), wsrt( ldt )
128 COMPLEX*16 CDUM( 1 ), LE( ldt, ldt ), RE( ldt, ldt ),
129 $ t( ldt, ldt ), tmp( ldt, ldt ), w( ldt ),
130 $ work( lwork ), wtmp( ldt )
133 DOUBLE PRECISION DLAMCH, ZLANGE
134 EXTERNAL dlamch, zlange
141 INTRINSIC dble, dimag, max, sqrt
146 smlnum = dlamch(
'S' ) / eps
147 bignum = one / smlnum
148 CALL dlabad( smlnum, bignum )
152 eps = max( eps, epsin )
163 val( 1 ) = sqrt( smlnum )
165 val( 3 ) = sqrt( bignum )
172 READ( nin, fmt = * )n, isrt
176 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
179 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
181 tnrm = zlange(
'M', n, n, tmp, ldt, rwork )
187 CALL zlacpy(
'F', n, n, tmp, ldt, t, ldt )
190 CALL zdscal( n, vmul, t( 1, i ), 1 )
197 CALL zgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
201 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL zhseqr(
'S',
'N', n, 1, n, t, ldt, w, cdum, 1, work,
216 ninfo( 2 ) = ninfo( 2 ) + 1
225 CALL ztrevc(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, n,
226 $ m, work, rwork, info )
230 CALL ztrsna(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, s,
231 $ sep, n, m, work, n, rwork, info )
234 ninfo( 3 ) = ninfo( 3 ) + 1
241 CALL zcopy( n, w, 1, wtmp, 1 )
247 wsrt( i ) = dble( w( i ) )
254 wsrt( i ) = dimag( w( i ) )
257 CALL dcopy( n, s, 1, stmp, 1 )
258 CALL dcopy( n, sep, 1, septmp, 1 )
259 CALL dscal( n, one / vmul, septmp, 1 )
264 IF( wsrt( j ).LT.vmin )
THEN
269 wsrt( kmin ) = wsrt( i )
272 wtmp( i ) = w( kmin )
275 stmp( kmin ) = stmp( i )
277 vmin = septmp( kmin )
278 septmp( kmin ) = septmp( i )
285 v = max( two*dble( n )*eps*tnrm, smlnum )
289 IF( v.GT.septmp( i ) )
THEN
292 tol = v / septmp( i )
294 IF( v.GT.sepin( i ) )
THEN
297 tolin = v / sepin( i )
299 tol = max( tol, smlnum / eps )
300 tolin = max( tolin, smlnum / eps )
301 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol )
THEN
303 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol )
THEN
304 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
305 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) )
THEN
307 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol )
THEN
308 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
312 IF( vmax.GT.rmax( 2 ) )
THEN
314 IF( ninfo( 2 ).EQ.0 )
323 IF( v.GT.septmp( i )*stmp( i ) )
THEN
328 IF( v.GT.sepin( i )*sin( i ) )
THEN
333 tol = max( tol, smlnum / eps )
334 tolin = max( tolin, smlnum / eps )
335 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol )
THEN
337 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol )
THEN
338 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
339 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) )
THEN
341 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol )
THEN
342 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
346 IF( vmax.GT.rmax( 2 ) )
THEN
348 IF( ninfo( 2 ).EQ.0 )
357 IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
358 $ dble( 2*n )*eps )
THEN
360 ELSE IF( eps*sin( i ).GT.stmp( i ) )
THEN
362 ELSE IF( sin( i ).GT.stmp( i ) )
THEN
363 vmax = sin( i ) / stmp( i )
364 ELSE IF( sin( i ).LT.eps*stmp( i ) )
THEN
366 ELSE IF( sin( i ).LT.stmp( i ) )
THEN
367 vmax = stmp( i ) / sin( i )
371 IF( vmax.GT.rmax( 3 ) )
THEN
373 IF( ninfo( 3 ).EQ.0 )
382 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v )
THEN
384 ELSE IF( eps*sepin( i ).GT.septmp( i ) )
THEN
386 ELSE IF( sepin( i ).GT.septmp( i ) )
THEN
387 vmax = sepin( i ) / septmp( i )
388 ELSE IF( sepin( i ).LT.eps*septmp( i ) )
THEN
390 ELSE IF( sepin( i ).LT.septmp( i ) )
THEN
391 vmax = septmp( i ) / sepin( i )
395 IF( vmax.GT.rmax( 3 ) )
THEN
397 IF( ninfo( 3 ).EQ.0 )
406 CALL dcopy( n, dum, 0, stmp, 1 )
407 CALL dcopy( n, dum, 0, septmp, 1 )
408 CALL ztrsna(
'E',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
409 $ stmp, septmp, n, m, work, n, rwork, info )
412 ninfo( 3 ) = ninfo( 3 ) + 1
416 IF( stmp( i ).NE.s( i ) )
418 IF( septmp( i ).NE.dum( 1 ) )
424 CALL dcopy( n, dum, 0, stmp, 1 )
425 CALL dcopy( n, dum, 0, septmp, 1 )
426 CALL ztrsna(
'V',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
427 $ stmp, septmp, n, m, work, n, rwork, info )
430 ninfo( 3 ) = ninfo( 3 ) + 1
434 IF( stmp( i ).NE.dum( 1 ) )
436 IF( septmp( i ).NE.sep( i ) )
445 CALL dcopy( n, dum, 0, stmp, 1 )
446 CALL dcopy( n, dum, 0, septmp, 1 )
447 CALL ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
448 $ stmp, septmp, n, m, work, n, rwork, info )
451 ninfo( 3 ) = ninfo( 3 ) + 1
455 IF( septmp( i ).NE.sep( i ) )
457 IF( stmp( i ).NE.s( i ) )
463 CALL dcopy( n, dum, 0, stmp, 1 )
464 CALL dcopy( n, dum, 0, septmp, 1 )
465 CALL ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
466 $ stmp, septmp, n, m, work, n, rwork, info )
469 ninfo( 3 ) = ninfo( 3 ) + 1
473 IF( stmp( i ).NE.s( i ) )
475 IF( septmp( i ).NE.dum( 1 ) )
481 CALL dcopy( n, dum, 0, stmp, 1 )
482 CALL dcopy( n, dum, 0, septmp, 1 )
483 CALL ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
484 $ stmp, septmp, n, m, work, n, rwork, info )
487 ninfo( 3 ) = ninfo( 3 ) + 1
491 IF( stmp( i ).NE.dum( 1 ) )
493 IF( septmp( i ).NE.sep( i ) )
496 IF( vmax.GT.rmax( 1 ) )
THEN
498 IF( ninfo( 1 ).EQ.0 )
505 SELECT( i ) = .false.
512 CALL zcopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
513 CALL zcopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
518 SELECT( n-1 ) = .true.
519 CALL zcopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
520 CALL zcopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
525 CALL dcopy( icmp, dum, 0, stmp, 1 )
526 CALL dcopy( icmp, dum, 0, septmp, 1 )
527 CALL ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
528 $ stmp, septmp, n, m, work, n, rwork, info )
531 ninfo( 3 ) = ninfo( 3 ) + 1
536 IF( septmp( i ).NE.sep( j ) )
538 IF( stmp( i ).NE.s( j ) )
544 CALL dcopy( icmp, dum, 0, stmp, 1 )
545 CALL dcopy( icmp, dum, 0, septmp, 1 )
546 CALL ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
547 $ stmp, septmp, n, m, work, n, rwork, info )
550 ninfo( 3 ) = ninfo( 3 ) + 1
555 IF( stmp( i ).NE.s( j ) )
557 IF( septmp( i ).NE.dum( 1 ) )
563 CALL dcopy( icmp, dum, 0, stmp, 1 )
564 CALL dcopy( icmp, dum, 0, septmp, 1 )
565 CALL ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
566 $ stmp, septmp, n, m, work, n, rwork, info )
569 ninfo( 3 ) = ninfo( 3 ) + 1
574 IF( stmp( i ).NE.dum( 1 ) )
576 IF( septmp( i ).NE.sep( j ) )
579 IF( vmax.GT.rmax( 1 ) )
THEN
581 IF( ninfo( 1 ).EQ.0 )
subroutine zget37(RMAX, LMAX, NINFO, KNT, NIN)
ZGET37
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine ztrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, INFO)
ZTRSNA
subroutine ztrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTREVC
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
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 dlabad(SMALL, LARGE)
DLABAD
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL