103 INTEGER lmax( 3 ), ninfo( 3 )
104 DOUBLE PRECISION rmax( 3 )
110 DOUBLE PRECISION zero, one, two
111 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
112 DOUBLE PRECISION epsin
113 parameter ( epsin = 5.9605d-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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dble, max, sqrt
149 smlnum =
dlamch(
'S' ) / eps
150 bignum = one / smlnum
151 CALL dlabad( 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 =
dlange(
'M', n, n, tmp, ldt, work )
191 CALL dlacpy(
'F', n, n, tmp, ldt, t, ldt )
194 CALL dscal( n, vmul, t( 1, i ), 1 )
198 CALL dlacpy(
'F', n, n, t, ldt, tsav, ldt )
202 CALL dgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
206 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL dlacpy(
'L', n, n, t, ldt, q, ldt )
213 CALL dorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
218 CALL dhseqr(
'S',
'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
222 ninfo( 2 ) = ninfo( 2 ) + 1
230 SELECT( i ) = .false.
232 CALL dcopy( n, wr, 1, wrtmp, 1 )
233 CALL dcopy( 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 dlacpy(
'F', n, n, q, ldt, qsav, ldt )
260 CALL dlacpy(
'F', n, n, t, ldt, tsav1, ldt )
261 CALL dtrsen(
'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 dhst01( 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*dble( 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.dble( 2*n )*eps .AND. stmp.LE.dble( 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 dlacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
398 CALL dlacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
401 CALL dtrsen(
'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 dlacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
426 CALL dlacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
429 CALL dtrsen(
'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 dlacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
454 CALL dlacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
457 CALL dtrsen(
'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 dlacpy(
'F', n, n, tsav1, ldt, ttmp, ldt )
482 CALL dlacpy(
'F', n, n, qsav, ldt, qtmp, ldt )
485 CALL dtrsen(
'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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(CMACH)
DLAMCH
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dtrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
DTRSEN
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
subroutine dlabad(SMALL, LARGE)
DLABAD