102 INTEGER lmax( 3 ), ninfo( 3 )
110 parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
112 parameter ( epsin = 5.9605e-8 )
114 parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
117 INTEGER i, icmp, info, iscl, isrt, j, kmin, m, n
118 REAL bignum, eps, smlnum, tnrm, tol, tolin, v,
119 $ vcmin, vmax, vmin, vmul
122 LOGICAL select( ldt )
124 REAL 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 cdum( 1 ), le( ldt, ldt ), re( ldt, ldt ),
129 $ t( ldt, ldt ), tmp( ldt, ldt ), w( ldt ),
130 $ work( lwork ), wtmp( ldt )
141 INTRINSIC aimag, max,
REAL, sqrt
146 smlnum =
slamch(
'S' ) / eps
147 bignum = one / smlnum
148 CALL slabad( 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 =
clange(
'M', n, n, tmp, ldt, rwork )
187 CALL clacpy(
'F', n, n, tmp, ldt, t, ldt )
190 CALL csscal( n, vmul, t( 1, i ), 1 )
197 CALL cgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
201 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL chseqr(
'S',
'N', n, 1, n, t, ldt, w, cdum, 1, work,
216 ninfo( 2 ) = ninfo( 2 ) + 1
225 CALL ctrevc(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, n,
226 $ m, work, rwork, info )
230 CALL ctrsna(
'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 ccopy( n, w, 1, wtmp, 1 )
247 wsrt( i ) =
REAL( W( I ) )
254 wsrt( i ) = aimag( w( i ) )
257 CALL scopy( n, s, 1, stmp, 1 )
258 CALL scopy( n, sep, 1, septmp, 1 )
259 CALL sscal( 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*
REAL( 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.
REAL( 2*n )*eps .AND. stmp( i ).LE.
358 $
REAL( 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 scopy( n, dum, 0, stmp, 1 )
407 CALL scopy( n, dum, 0, septmp, 1 )
408 CALL ctrsna(
'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 scopy( n, dum, 0, stmp, 1 )
425 CALL scopy( n, dum, 0, septmp, 1 )
426 CALL ctrsna(
'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 scopy( n, dum, 0, stmp, 1 )
446 CALL scopy( n, dum, 0, septmp, 1 )
447 CALL ctrsna(
'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 scopy( n, dum, 0, stmp, 1 )
464 CALL scopy( n, dum, 0, septmp, 1 )
465 CALL ctrsna(
'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 scopy( n, dum, 0, stmp, 1 )
482 CALL scopy( n, dum, 0, septmp, 1 )
483 CALL ctrsna(
'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 ccopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
513 CALL ccopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
518 SELECT( n-1 ) = .true.
519 CALL ccopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
520 CALL ccopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
525 CALL scopy( icmp, dum, 0, stmp, 1 )
526 CALL scopy( icmp, dum, 0, septmp, 1 )
527 CALL ctrsna(
'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 scopy( icmp, dum, 0, stmp, 1 )
545 CALL scopy( icmp, dum, 0, septmp, 1 )
546 CALL ctrsna(
'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 scopy( icmp, dum, 0, stmp, 1 )
564 CALL scopy( icmp, dum, 0, septmp, 1 )
565 CALL ctrsna(
'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 chseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
CHSEQR
subroutine ctrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTREVC
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine ctrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, INFO)
CTRSNA
real function slamch(CMACH)
SLAMCH
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine csscal(N, SA, CX, INCX)
CSSCAL