91 SUBROUTINE zget37( RMAX, LMAX, NINFO, KNT, NIN )
102 INTEGER lmax( 3 ), ninfo( 3 )
103 DOUBLE PRECISION rmax( 3 )
109 DOUBLE PRECISION zero, one, two
110 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
111 DOUBLE PRECISION epsin
112 parameter( epsin = 5.9605d-8 )
114 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
117 INTEGER i, icmp, info, iscl, isrt,
j, kmin, m, n
118 DOUBLE PRECISION bignum, eps, smlnum, tnrm, tol, tolin, v,
119 $ vcmin, vmax, vmin, vmul
122 LOGICAL select( ldt )
124 DOUBLE PRECISION dum( 1 ), rwork( 2*ldt ), s( ldt ), sep( ldt ),
125 $ sepin( ldt ), septmp( ldt ), sin( ldt ),
126 $ stmp( ldt ), val( 3 ), wiin( ldt ),
127 $ wrin( ldt ), wsrt( ldt )
128 COMPLEX*16 cdum( 1 ), le( ldt, ldt ), re( ldt, ldt ),
129 $ t( ldt, ldt ), tmp( ldt, ldt ), w( ldt ),
130 $ work( lwork ), wtmp( ldt )
141 INTRINSIC dble, dimag, max, sqrt
146 smlnum =
dlamch(
'S' ) / eps
147 bignum = one / smlnum
148 CALL
dlabad( smlnum, bignum )
152 eps = max( eps, epsin )
163 val( 1 ) = sqrt( smlnum )
165 val( 3 ) = sqrt( bignum )
172 READ( nin, fmt = * )n, isrt
176 READ( nin, fmt = * )( tmp( i,
j ),
j = 1, n )
179 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
181 tnrm =
zlange(
'M', n, n, tmp, ldt, rwork )
187 CALL
zlacpy(
'F', n, n, tmp, ldt, t, ldt )
190 CALL
zdscal( n, vmul, t( 1, i ), 1 )
197 CALL
zgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
201 ninfo( 1 ) = ninfo( 1 ) + 1
212 CALL
zhseqr(
'S',
'N', n, 1, n, t, ldt, w, cdum, 1, work,
216 ninfo( 2 ) = ninfo( 2 ) + 1
225 CALL
ztrevc(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, n,
226 $ m, work, rwork, info )
230 CALL
ztrsna(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, s,
231 $ sep, n, m, work, n, rwork, info )
234 ninfo( 3 ) = ninfo( 3 ) + 1
241 CALL
zcopy( n, w, 1, wtmp, 1 )
247 wsrt( i ) = dble( w( i ) )
254 wsrt( i ) = dimag( w( i ) )
257 CALL
dcopy( n, s, 1, stmp, 1 )
258 CALL
dcopy( n, sep, 1, septmp, 1 )
259 CALL
dscal( n, one / vmul, septmp, 1 )
264 IF( wsrt(
j ).LT.vmin )
THEN
269 wsrt( kmin ) = wsrt( i )
272 wtmp( i ) = w( kmin )
275 stmp( kmin ) = stmp( i )
277 vmin = septmp( kmin )
278 septmp( kmin ) = septmp( i )
285 v = max( two*dble( n )*eps*tnrm, smlnum )
289 IF( v.GT.septmp( i ) )
THEN
292 tol = v / septmp( i )
294 IF( v.GT.sepin( i ) )
THEN
297 tolin = v / sepin( i )
299 tol = max( tol, smlnum / eps )
300 tolin = max( tolin, smlnum / eps )
301 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol )
THEN
303 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol )
THEN
304 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
305 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) )
THEN
307 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol )
THEN
308 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
312 IF( vmax.GT.rmax( 2 ) )
THEN
314 IF( ninfo( 2 ).EQ.0 )
323 IF( v.GT.septmp( i )*stmp( i ) )
THEN
328 IF( v.GT.sepin( i )*sin( i ) )
THEN
333 tol = max( tol, smlnum / eps )
334 tolin = max( tolin, smlnum / eps )
335 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol )
THEN
337 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol )
THEN
338 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
339 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) )
THEN
341 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol )
THEN
342 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
346 IF( vmax.GT.rmax( 2 ) )
THEN
348 IF( ninfo( 2 ).EQ.0 )
357 IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
358 $ dble( 2*n )*eps )
THEN
360 ELSE IF( eps*sin( i ).GT.stmp( i ) )
THEN
362 ELSE IF( sin( i ).GT.stmp( i ) )
THEN
363 vmax = sin( i ) / stmp( i )
364 ELSE IF( sin( i ).LT.eps*stmp( i ) )
THEN
366 ELSE IF( sin( i ).LT.stmp( i ) )
THEN
367 vmax = stmp( i ) / sin( i )
371 IF( vmax.GT.rmax( 3 ) )
THEN
373 IF( ninfo( 3 ).EQ.0 )
382 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v )
THEN
384 ELSE IF( eps*sepin( i ).GT.septmp( i ) )
THEN
386 ELSE IF( sepin( i ).GT.septmp( i ) )
THEN
387 vmax = sepin( i ) / septmp( i )
388 ELSE IF( sepin( i ).LT.eps*septmp( i ) )
THEN
390 ELSE IF( sepin( i ).LT.septmp( i ) )
THEN
391 vmax = septmp( i ) / sepin( i )
395 IF( vmax.GT.rmax( 3 ) )
THEN
397 IF( ninfo( 3 ).EQ.0 )
406 CALL
dcopy( n, dum, 0, stmp, 1 )
407 CALL
dcopy( n, dum, 0, septmp, 1 )
408 CALL
ztrsna(
'E',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
409 $ stmp, septmp, n, m, work, n, rwork, info )
412 ninfo( 3 ) = ninfo( 3 ) + 1
416 IF( stmp( i ).NE.s( i ) )
418 IF( septmp( i ).NE.dum( 1 ) )
424 CALL
dcopy( n, dum, 0, stmp, 1 )
425 CALL
dcopy( n, dum, 0, septmp, 1 )
426 CALL
ztrsna(
'V',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
427 $ stmp, septmp, n, m, work, n, rwork, info )
430 ninfo( 3 ) = ninfo( 3 ) + 1
434 IF( stmp( i ).NE.dum( 1 ) )
436 IF( septmp( i ).NE.sep( i ) )
445 CALL
dcopy( n, dum, 0, stmp, 1 )
446 CALL
dcopy( n, dum, 0, septmp, 1 )
447 CALL
ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
448 $ stmp, septmp, n, m, work, n, rwork, info )
451 ninfo( 3 ) = ninfo( 3 ) + 1
455 IF( septmp( i ).NE.sep( i ) )
457 IF( stmp( i ).NE.s( i ) )
463 CALL
dcopy( n, dum, 0, stmp, 1 )
464 CALL
dcopy( n, dum, 0, septmp, 1 )
465 CALL
ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
466 $ stmp, septmp, n, m, work, n, rwork, info )
469 ninfo( 3 ) = ninfo( 3 ) + 1
473 IF( stmp( i ).NE.s( i ) )
475 IF( septmp( i ).NE.dum( 1 ) )
481 CALL
dcopy( n, dum, 0, stmp, 1 )
482 CALL
dcopy( n, dum, 0, septmp, 1 )
483 CALL
ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
484 $ stmp, septmp, n, m, work, n, rwork, info )
487 ninfo( 3 ) = ninfo( 3 ) + 1
491 IF( stmp( i ).NE.dum( 1 ) )
493 IF( septmp( i ).NE.sep( i ) )
496 IF( vmax.GT.rmax( 1 ) )
THEN
498 IF( ninfo( 1 ).EQ.0 )
505 SELECT( i ) = .false.
512 CALL
zcopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
513 CALL
zcopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
518 SELECT( n-1 ) = .true.
519 CALL
zcopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
520 CALL
zcopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
525 CALL
dcopy( icmp, dum, 0, stmp, 1 )
526 CALL
dcopy( icmp, dum, 0, septmp, 1 )
527 CALL
ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
528 $ stmp, septmp, n, m, work, n, rwork, info )
531 ninfo( 3 ) = ninfo( 3 ) + 1
536 IF( septmp( i ).NE.sep(
j ) )
538 IF( stmp( i ).NE.s(
j ) )
544 CALL
dcopy( icmp, dum, 0, stmp, 1 )
545 CALL
dcopy( icmp, dum, 0, septmp, 1 )
546 CALL
ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
547 $ stmp, septmp, n, m, work, n, rwork, info )
550 ninfo( 3 ) = ninfo( 3 ) + 1
555 IF( stmp( i ).NE.s(
j ) )
557 IF( septmp( i ).NE.dum( 1 ) )
563 CALL
dcopy( icmp, dum, 0, stmp, 1 )
564 CALL
dcopy( icmp, dum, 0, septmp, 1 )
565 CALL
ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
566 $ stmp, septmp, n, m, work, n, rwork, info )
569 ninfo( 3 ) = ninfo( 3 ) + 1
574 IF( stmp( i ).NE.dum( 1 ) )
576 IF( septmp( i ).NE.sep(
j ) )
579 IF( vmax.GT.rmax( 1 ) )
THEN
581 IF( ninfo( 1 ).EQ.0 )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine ztrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTREVC
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine zget37(RMAX, LMAX, NINFO, KNT, NIN)
ZGET37
subroutine ztrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, INFO)
ZTRSNA
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE 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