103 INTEGER lmax( 3 ), ninfo( 3 )
111 parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
113 parameter ( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
115 parameter ( epsin = 5.9605e-8 )
117 parameter ( czero = ( 0.0e+0, 0.0e+0 ) )
120 INTEGER i, info, iscl, isrt, itmp, j, kmin, m, n, ndim
121 REAL bignum, eps, s, sep, sepin, septmp, sin,
122 $ smlnum, stmp, tnrm, tol, tolin, v, vmax, vmin,
126 LOGICAL select( ldt )
127 INTEGER ipnt( ldt ), iselec( ldt )
128 REAL result( 2 ), rwork( ldt ), val( 3 ),
130 COMPLEX q( ldt, ldt ), qsav( ldt, ldt ),
131 $ qtmp( ldt, ldt ), t( ldt, ldt ),
132 $ tmp( ldt, ldt ), tsav( ldt, ldt ),
133 $ tsav1( ldt, ldt ), ttmp( ldt, ldt ), w( ldt ),
134 $ work( lwork ), wtmp( ldt )
145 INTRINSIC aimag, max,
REAL, sqrt
150 smlnum =
slamch(
'S' ) / eps
151 bignum = one / smlnum
152 CALL slabad( smlnum, bignum )
156 eps = max( eps, epsin )
167 val( 1 ) = sqrt( smlnum )
169 val( 3 ) = sqrt( sqrt( bignum ) )
176 READ( nin, fmt = * )n, ndim, isrt
179 READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
181 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
183 READ( nin, fmt = * )sin, sepin
185 tnrm =
clange(
'M', n, n, tmp, ldt, rwork )
191 CALL clacpy(
'F', n, n, tmp, ldt, t, ldt )
194 CALL csscal( n, vmul, t( 1, i ), 1 )
198 CALL clacpy(
'F', n, n, t, ldt, tsav, ldt )
202 CALL cgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
206 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL clacpy(
'L', n, n, t, ldt, q, ldt )
213 CALL cunghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
223 CALL chseqr(
'S',
'V', n, 1, n, t, ldt, w, q, ldt, work, lwork,
227 ninfo( 2 ) = ninfo( 2 ) + 1
235 SELECT( i ) = .false.
239 wsrt( i ) =
REAL( W( I ) )
243 wsrt( i ) = aimag( w( i ) )
250 IF( wsrt( j ).LT.vmin )
THEN
255 wsrt( kmin ) = wsrt( i )
258 ipnt( i ) = ipnt( kmin )
262 SELECT( ipnt( iselec( i ) ) ) = .true.
267 CALL clacpy(
'F', n, n, q, ldt, qsav, ldt )
268 CALL clacpy(
'F', n, n, t, ldt, tsav1, ldt )
269 CALL ctrsen(
'B',
'V',
SELECT, n, t, ldt, q, ldt, wtmp, m, s,
270 $ sep, work, lwork, info )
273 ninfo( 3 ) = ninfo( 3 ) + 1
281 CALL chst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
283 vmax = max( result( 1 ), result( 2 ) )
284 IF( vmax.GT.rmax( 1 ) )
THEN
286 IF( ninfo( 1 ).EQ.0 )
293 v = max( two*
REAL( n )*eps*tnrm, smlnum )
296 IF( v.GT.septmp )
THEN
301 IF( v.GT.sepin )
THEN
306 tol = max( tol, smlnum / eps )
307 tolin = max( tolin, smlnum / eps )
308 IF( eps*( sin-tolin ).GT.stmp+tol )
THEN
310 ELSE IF( sin-tolin.GT.stmp+tol )
THEN
311 vmax = ( sin-tolin ) / ( stmp+tol )
312 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) )
THEN
314 ELSE IF( sin+tolin.LT.stmp-tol )
THEN
315 vmax = ( stmp-tol ) / ( sin+tolin )
319 IF( vmax.GT.rmax( 2 ) )
THEN
321 IF( ninfo( 2 ).EQ.0 )
328 IF( v.GT.septmp*stmp )
THEN
333 IF( v.GT.sepin*sin )
THEN
338 tol = max( tol, smlnum / eps )
339 tolin = max( tolin, smlnum / eps )
340 IF( eps*( sepin-tolin ).GT.septmp+tol )
THEN
342 ELSE IF( sepin-tolin.GT.septmp+tol )
THEN
343 vmax = ( sepin-tolin ) / ( septmp+tol )
344 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) )
THEN
346 ELSE IF( sepin+tolin.LT.septmp-tol )
THEN
347 vmax = ( septmp-tol ) / ( sepin+tolin )
351 IF( vmax.GT.rmax( 2 ) )
THEN
353 IF( ninfo( 2 ).EQ.0 )
360 IF( sin.LE.
REAL( 2*n )*eps .AND. stmp.LE.
REAL( 2*n )*eps ) then
362 ELSE IF( eps*sin.GT.stmp )
THEN
364 ELSE IF( sin.GT.stmp )
THEN
366 ELSE IF( sin.LT.eps*stmp )
THEN
368 ELSE IF( sin.LT.stmp )
THEN
373 IF( vmax.GT.rmax( 3 ) )
THEN
375 IF( ninfo( 3 ).EQ.0 )
382 IF( sepin.LE.v .AND. septmp.LE.v )
THEN
384 ELSE IF( eps*sepin.GT.septmp )
THEN
386 ELSE IF( sepin.GT.septmp )
THEN
387 vmax = sepin / septmp
388 ELSE IF( sepin.LT.eps*septmp )
THEN
390 ELSE IF( sepin.LT.septmp )
THEN
391 vmax = septmp / sepin
395 IF( vmax.GT.rmax( 3 ) )
THEN
397 IF( ninfo( 3 ).EQ.0 )
405 CALL clacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
406 CALL clacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
409 CALL ctrsen(
'E',
'V',
SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
410 $ m, stmp, septmp, work, lwork, info )
413 ninfo( 3 ) = ninfo( 3 ) + 1
422 IF( ttmp( i, j ).NE.t( i, j ) )
424 IF( qtmp( i, j ).NE.q( i, j ) )
432 CALL clacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
433 CALL clacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
436 CALL ctrsen(
'V',
'V',
SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
437 $ m, stmp, septmp, work, lwork, info )
440 ninfo( 3 ) = ninfo( 3 ) + 1
449 IF( ttmp( i, j ).NE.t( i, j ) )
451 IF( qtmp( i, j ).NE.q( i, j ) )
459 CALL clacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
460 CALL clacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
463 CALL ctrsen(
'E',
'N',
SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
464 $ m, stmp, septmp, work, lwork, info )
467 ninfo( 3 ) = ninfo( 3 ) + 1
476 IF( ttmp( i, j ).NE.t( i, j ) )
478 IF( qtmp( i, j ).NE.qsav( i, j ) )
486 CALL clacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
487 CALL clacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
490 CALL ctrsen(
'V',
'N',
SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
491 $ m, stmp, septmp, work, lwork, info )
494 ninfo( 3 ) = ninfo( 3 ) + 1
503 IF( ttmp( i, j ).NE.t( i, j ) )
505 IF( qtmp( i, j ).NE.qsav( i, j ) )
509 IF( vmax.GT.rmax( 1 ) )
THEN
511 IF( ninfo( 1 ).EQ.0 )
subroutine chseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
CHSEQR
subroutine ctrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
CTRSEN
subroutine chst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
CHST01
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 ...
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 cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
subroutine cunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CUNGHR
subroutine csscal(N, SA, CX, INCX)
CSSCAL