92 SUBROUTINE dget38( RMAX, LMAX, NINFO, KNT, NIN )
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 dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
subroutine dget38(RMAX, LMAX, NINFO, KNT, NIN)
DGET38
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dtrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
DTRSEN
subroutine dscal(N, DA, DX, INCX)
DSCAL
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 ...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR