92 SUBROUTINE cget38( RMAX, LMAX, NINFO, KNT, NIN )
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 )
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 ctrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
CTRSEN
subroutine cget38(RMAX, LMAX, NINFO, KNT, NIN)
CGET38
subroutine chst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
CHST01
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine chseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
CHSEQR
real function slamch(CMACH)
SLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine csscal(N, SA, CX, INCX)
CSSCAL
subroutine cunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CUNGHR