103 INTEGER lmax( 3 ), ninfo( 3 )
111 parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
113 parameter ( epsin = 5.9605e-8 )
115 parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
117 parameter ( liwork = ldt*ldt )
120 INTEGER i, info, iscl, itmp, j, kmin, m, n, ndim
121 REAL bignum, eps, s, sep, sepin, septmp, sin,
122 $ smlnum, stmp, tnrm, tol, tolin, v, vimin, vmax,
126 LOGICAL select( ldt )
127 INTEGER ipnt( ldt ), iselec( ldt ), iwork( liwork )
128 REAL q( ldt, ldt ), qsav( ldt, ldt ),
129 $ qtmp( ldt, ldt ), result( 2 ), t( ldt, ldt ),
130 $ tmp( ldt, ldt ), tsav( ldt, ldt ),
131 $ tsav1( ldt, ldt ), ttmp( ldt, ldt ), val( 3 ),
132 $ wi( ldt ), witmp( ldt ), work( lwork ),
133 $ wr( ldt ), wrtmp( ldt )
144 INTRINSIC max,
REAL, sqrt
149 smlnum =
slamch(
'S' ) / eps
150 bignum = one / smlnum
151 CALL slabad( smlnum, bignum )
155 eps = max( eps, epsin )
167 val( 1 ) = sqrt( smlnum )
169 val( 3 ) = sqrt( sqrt( bignum ) )
176 READ( nin, fmt = * )n, ndim
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 =
slange(
'M', n, n, tmp, ldt, work )
191 CALL slacpy(
'F', n, n, tmp, ldt, t, ldt )
194 CALL sscal( n, vmul, t( 1, i ), 1 )
198 CALL slacpy(
'F', n, n, t, ldt, tsav, ldt )
202 CALL sgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
206 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL slacpy(
'L', n, n, t, ldt, q, ldt )
213 CALL sorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
218 CALL shseqr(
'S',
'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
222 ninfo( 2 ) = ninfo( 2 ) + 1
230 SELECT( i ) = .false.
232 CALL scopy( n, wr, 1, wrtmp, 1 )
233 CALL scopy( n, wi, 1, witmp, 1 )
239 IF( wrtmp( j ).LT.vrmin )
THEN
245 wrtmp( kmin ) = wrtmp( i )
246 witmp( kmin ) = witmp( i )
250 ipnt( i ) = ipnt( kmin )
254 SELECT( ipnt( iselec( i ) ) ) = .true.
259 CALL slacpy(
'F', n, n, q, ldt, qsav, ldt )
260 CALL slacpy(
'F', n, n, t, ldt, tsav1, ldt )
261 CALL strsen(
'B',
'V',
SELECT, n, t, ldt, q, ldt, wrtmp, witmp,
262 $ m, s, sep, work, lwork, iwork, liwork, info )
265 ninfo( 3 ) = ninfo( 3 ) + 1
273 CALL shst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
275 vmax = max( result( 1 ), result( 2 ) )
276 IF( vmax.GT.rmax( 1 ) )
THEN
278 IF( ninfo( 1 ).EQ.0 )
285 v = max( two*
REAL( n )*eps*tnrm, smlnum )
288 IF( v.GT.septmp )
THEN
293 IF( v.GT.sepin )
THEN
298 tol = max( tol, smlnum / eps )
299 tolin = max( tolin, smlnum / eps )
300 IF( eps*( sin-tolin ).GT.stmp+tol )
THEN
302 ELSE IF( sin-tolin.GT.stmp+tol )
THEN
303 vmax = ( sin-tolin ) / ( stmp+tol )
304 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) )
THEN
306 ELSE IF( sin+tolin.LT.stmp-tol )
THEN
307 vmax = ( stmp-tol ) / ( sin+tolin )
311 IF( vmax.GT.rmax( 2 ) )
THEN
313 IF( ninfo( 2 ).EQ.0 )
320 IF( v.GT.septmp*stmp )
THEN
325 IF( v.GT.sepin*sin )
THEN
330 tol = max( tol, smlnum / eps )
331 tolin = max( tolin, smlnum / eps )
332 IF( eps*( sepin-tolin ).GT.septmp+tol )
THEN
334 ELSE IF( sepin-tolin.GT.septmp+tol )
THEN
335 vmax = ( sepin-tolin ) / ( septmp+tol )
336 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) )
THEN
338 ELSE IF( sepin+tolin.LT.septmp-tol )
THEN
339 vmax = ( septmp-tol ) / ( sepin+tolin )
343 IF( vmax.GT.rmax( 2 ) )
THEN
345 IF( ninfo( 2 ).EQ.0 )
352 IF( sin.LE.
REAL( 2*n )*eps .AND. stmp.LE.
REAL( 2*n )*eps ) then
354 ELSE IF( eps*sin.GT.stmp )
THEN
356 ELSE IF( sin.GT.stmp )
THEN
358 ELSE IF( sin.LT.eps*stmp )
THEN
360 ELSE IF( sin.LT.stmp )
THEN
365 IF( vmax.GT.rmax( 3 ) )
THEN
367 IF( ninfo( 3 ).EQ.0 )
374 IF( sepin.LE.v .AND. septmp.LE.v )
THEN
376 ELSE IF( eps*sepin.GT.septmp )
THEN
378 ELSE IF( sepin.GT.septmp )
THEN
379 vmax = sepin / septmp
380 ELSE IF( sepin.LT.eps*septmp )
THEN
382 ELSE IF( sepin.LT.septmp )
THEN
383 vmax = septmp / sepin
387 IF( vmax.GT.rmax( 3 ) )
THEN
389 IF( ninfo( 3 ).EQ.0 )
397 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
398 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
401 CALL strsen(
'E',
'V',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
402 $ witmp, m, stmp, septmp, work, lwork, iwork,
406 ninfo( 3 ) = ninfo( 3 ) + 1
415 IF( ttmp( i, j ).NE.t( i, j ) )
417 IF( qtmp( i, j ).NE.q( i, j ) )
425 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
426 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
429 CALL strsen(
'V',
'V',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
430 $ witmp, m, stmp, septmp, work, lwork, iwork,
434 ninfo( 3 ) = ninfo( 3 ) + 1
443 IF( ttmp( i, j ).NE.t( i, j ) )
445 IF( qtmp( i, j ).NE.q( i, j ) )
453 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
454 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
457 CALL strsen(
'E',
'N',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
458 $ witmp, m, stmp, septmp, work, lwork, iwork,
462 ninfo( 3 ) = ninfo( 3 ) + 1
471 IF( ttmp( i, j ).NE.t( i, j ) )
473 IF( qtmp( i, j ).NE.qsav( i, j ) )
481 CALL slacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
482 CALL slacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
485 CALL strsen(
'V',
'N',
SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
486 $ witmp, m, stmp, septmp, work, lwork, iwork,
490 ninfo( 3 ) = ninfo( 3 ) + 1
499 IF( ttmp( i, j ).NE.t( i, j ) )
501 IF( qtmp( i, j ).NE.qsav( i, j ) )
505 IF( vmax.GT.rmax( 1 ) )
THEN
507 IF( ninfo( 1 ).EQ.0 )
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine strsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
STRSEN
real function slamch(CMACH)
SLAMCH
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR