00001 SUBROUTINE ZGET37( RMAX, LMAX, NINFO, KNT, NIN )
00002
00003
00004
00005
00006
00007
00008 INTEGER KNT, NIN
00009
00010
00011 INTEGER LMAX( 3 ), NINFO( 3 )
00012 DOUBLE PRECISION RMAX( 3 )
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 DOUBLE PRECISION ZERO, ONE, TWO
00057 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00058 DOUBLE PRECISION EPSIN
00059 PARAMETER ( EPSIN = 5.9605D-8 )
00060 INTEGER LDT, LWORK
00061 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
00062
00063
00064 INTEGER I, ICMP, INFO, ISCL, ISRT, J, KMIN, M, N
00065 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
00066 $ VCMIN, VMAX, VMIN, VMUL
00067
00068
00069 LOGICAL SELECT( LDT )
00070 INTEGER LCMP( 3 )
00071 DOUBLE PRECISION DUM( 1 ), RWORK( 2*LDT ), S( LDT ), SEP( LDT ),
00072 $ SEPIN( LDT ), SEPTMP( LDT ), SIN( LDT ),
00073 $ STMP( LDT ), VAL( 3 ), WIIN( LDT ),
00074 $ WRIN( LDT ), WSRT( LDT )
00075 COMPLEX*16 CDUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
00076 $ T( LDT, LDT ), TMP( LDT, LDT ), W( LDT ),
00077 $ WORK( LWORK ), WTMP( LDT )
00078
00079
00080 DOUBLE PRECISION DLAMCH, ZLANGE
00081 EXTERNAL DLAMCH, ZLANGE
00082
00083
00084 EXTERNAL DCOPY, DLABAD, DSCAL, ZCOPY, ZDSCAL, ZGEHRD,
00085 $ ZHSEQR, ZLACPY, ZTREVC, ZTRSNA
00086
00087
00088 INTRINSIC DBLE, DIMAG, MAX, SQRT
00089
00090
00091
00092 EPS = DLAMCH( 'P' )
00093 SMLNUM = DLAMCH( 'S' ) / EPS
00094 BIGNUM = ONE / SMLNUM
00095 CALL DLABAD( SMLNUM, BIGNUM )
00096
00097
00098
00099 EPS = MAX( EPS, EPSIN )
00100 RMAX( 1 ) = ZERO
00101 RMAX( 2 ) = ZERO
00102 RMAX( 3 ) = ZERO
00103 LMAX( 1 ) = 0
00104 LMAX( 2 ) = 0
00105 LMAX( 3 ) = 0
00106 KNT = 0
00107 NINFO( 1 ) = 0
00108 NINFO( 2 ) = 0
00109 NINFO( 3 ) = 0
00110 VAL( 1 ) = SQRT( SMLNUM )
00111 VAL( 2 ) = ONE
00112 VAL( 3 ) = SQRT( BIGNUM )
00113
00114
00115
00116
00117
00118 10 CONTINUE
00119 READ( NIN, FMT = * )N, ISRT
00120 IF( N.EQ.0 )
00121 $ RETURN
00122 DO 20 I = 1, N
00123 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
00124 20 CONTINUE
00125 DO 30 I = 1, N
00126 READ( NIN, FMT = * )WRIN( I ), WIIN( I ), SIN( I ), SEPIN( I )
00127 30 CONTINUE
00128 TNRM = ZLANGE( 'M', N, N, TMP, LDT, RWORK )
00129 DO 260 ISCL = 1, 3
00130
00131
00132
00133 KNT = KNT + 1
00134 CALL ZLACPY( 'F', N, N, TMP, LDT, T, LDT )
00135 VMUL = VAL( ISCL )
00136 DO 40 I = 1, N
00137 CALL ZDSCAL( N, VMUL, T( 1, I ), 1 )
00138 40 CONTINUE
00139 IF( TNRM.EQ.ZERO )
00140 $ VMUL = ONE
00141
00142
00143
00144 CALL ZGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
00145 $ INFO )
00146 IF( INFO.NE.0 ) THEN
00147 LMAX( 1 ) = KNT
00148 NINFO( 1 ) = NINFO( 1 ) + 1
00149 GO TO 260
00150 END IF
00151 DO 60 J = 1, N - 2
00152 DO 50 I = J + 2, N
00153 T( I, J ) = ZERO
00154 50 CONTINUE
00155 60 CONTINUE
00156
00157
00158
00159 CALL ZHSEQR( 'S', 'N', N, 1, N, T, LDT, W, CDUM, 1, WORK,
00160 $ LWORK, INFO )
00161 IF( INFO.NE.0 ) THEN
00162 LMAX( 2 ) = KNT
00163 NINFO( 2 ) = NINFO( 2 ) + 1
00164 GO TO 260
00165 END IF
00166
00167
00168
00169 DO 70 I = 1, N
00170 SELECT( I ) = .TRUE.
00171 70 CONTINUE
00172 CALL ZTREVC( 'B', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT, N,
00173 $ M, WORK, RWORK, INFO )
00174
00175
00176
00177 CALL ZTRSNA( 'B', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT, S,
00178 $ SEP, N, M, WORK, N, RWORK, INFO )
00179 IF( INFO.NE.0 ) THEN
00180 LMAX( 3 ) = KNT
00181 NINFO( 3 ) = NINFO( 3 ) + 1
00182 GO TO 260
00183 END IF
00184
00185
00186
00187
00188 CALL ZCOPY( N, W, 1, WTMP, 1 )
00189 IF( ISRT.EQ.0 ) THEN
00190
00191
00192
00193 DO 80 I = 1, N
00194 WSRT( I ) = DBLE( W( I ) )
00195 80 CONTINUE
00196 ELSE
00197
00198
00199
00200 DO 90 I = 1, N
00201 WSRT( I ) = DIMAG( W( I ) )
00202 90 CONTINUE
00203 END IF
00204 CALL DCOPY( N, S, 1, STMP, 1 )
00205 CALL DCOPY( N, SEP, 1, SEPTMP, 1 )
00206 CALL DSCAL( N, ONE / VMUL, SEPTMP, 1 )
00207 DO 110 I = 1, N - 1
00208 KMIN = I
00209 VMIN = WSRT( I )
00210 DO 100 J = I + 1, N
00211 IF( WSRT( J ).LT.VMIN ) THEN
00212 KMIN = J
00213 VMIN = WSRT( J )
00214 END IF
00215 100 CONTINUE
00216 WSRT( KMIN ) = WSRT( I )
00217 WSRT( I ) = VMIN
00218 VCMIN = WTMP( I )
00219 WTMP( I ) = W( KMIN )
00220 WTMP( KMIN ) = VCMIN
00221 VMIN = STMP( KMIN )
00222 STMP( KMIN ) = STMP( I )
00223 STMP( I ) = VMIN
00224 VMIN = SEPTMP( KMIN )
00225 SEPTMP( KMIN ) = SEPTMP( I )
00226 SEPTMP( I ) = VMIN
00227 110 CONTINUE
00228
00229
00230
00231
00232 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM )
00233 IF( TNRM.EQ.ZERO )
00234 $ V = ONE
00235 DO 120 I = 1, N
00236 IF( V.GT.SEPTMP( I ) ) THEN
00237 TOL = ONE
00238 ELSE
00239 TOL = V / SEPTMP( I )
00240 END IF
00241 IF( V.GT.SEPIN( I ) ) THEN
00242 TOLIN = ONE
00243 ELSE
00244 TOLIN = V / SEPIN( I )
00245 END IF
00246 TOL = MAX( TOL, SMLNUM / EPS )
00247 TOLIN = MAX( TOLIN, SMLNUM / EPS )
00248 IF( EPS*( SIN( I )-TOLIN ).GT.STMP( I )+TOL ) THEN
00249 VMAX = ONE / EPS
00250 ELSE IF( SIN( I )-TOLIN.GT.STMP( I )+TOL ) THEN
00251 VMAX = ( SIN( I )-TOLIN ) / ( STMP( I )+TOL )
00252 ELSE IF( SIN( I )+TOLIN.LT.EPS*( STMP( I )-TOL ) ) THEN
00253 VMAX = ONE / EPS
00254 ELSE IF( SIN( I )+TOLIN.LT.STMP( I )-TOL ) THEN
00255 VMAX = ( STMP( I )-TOL ) / ( SIN( I )+TOLIN )
00256 ELSE
00257 VMAX = ONE
00258 END IF
00259 IF( VMAX.GT.RMAX( 2 ) ) THEN
00260 RMAX( 2 ) = VMAX
00261 IF( NINFO( 2 ).EQ.0 )
00262 $ LMAX( 2 ) = KNT
00263 END IF
00264 120 CONTINUE
00265
00266
00267
00268
00269 DO 130 I = 1, N
00270 IF( V.GT.SEPTMP( I )*STMP( I ) ) THEN
00271 TOL = SEPTMP( I )
00272 ELSE
00273 TOL = V / STMP( I )
00274 END IF
00275 IF( V.GT.SEPIN( I )*SIN( I ) ) THEN
00276 TOLIN = SEPIN( I )
00277 ELSE
00278 TOLIN = V / SIN( I )
00279 END IF
00280 TOL = MAX( TOL, SMLNUM / EPS )
00281 TOLIN = MAX( TOLIN, SMLNUM / EPS )
00282 IF( EPS*( SEPIN( I )-TOLIN ).GT.SEPTMP( I )+TOL ) THEN
00283 VMAX = ONE / EPS
00284 ELSE IF( SEPIN( I )-TOLIN.GT.SEPTMP( I )+TOL ) THEN
00285 VMAX = ( SEPIN( I )-TOLIN ) / ( SEPTMP( I )+TOL )
00286 ELSE IF( SEPIN( I )+TOLIN.LT.EPS*( SEPTMP( I )-TOL ) ) THEN
00287 VMAX = ONE / EPS
00288 ELSE IF( SEPIN( I )+TOLIN.LT.SEPTMP( I )-TOL ) THEN
00289 VMAX = ( SEPTMP( I )-TOL ) / ( SEPIN( I )+TOLIN )
00290 ELSE
00291 VMAX = ONE
00292 END IF
00293 IF( VMAX.GT.RMAX( 2 ) ) THEN
00294 RMAX( 2 ) = VMAX
00295 IF( NINFO( 2 ).EQ.0 )
00296 $ LMAX( 2 ) = KNT
00297 END IF
00298 130 CONTINUE
00299
00300
00301
00302
00303 DO 140 I = 1, N
00304 IF( SIN( I ).LE.DBLE( 2*N )*EPS .AND. STMP( I ).LE.
00305 $ DBLE( 2*N )*EPS ) THEN
00306 VMAX = ONE
00307 ELSE IF( EPS*SIN( I ).GT.STMP( I ) ) THEN
00308 VMAX = ONE / EPS
00309 ELSE IF( SIN( I ).GT.STMP( I ) ) THEN
00310 VMAX = SIN( I ) / STMP( I )
00311 ELSE IF( SIN( I ).LT.EPS*STMP( I ) ) THEN
00312 VMAX = ONE / EPS
00313 ELSE IF( SIN( I ).LT.STMP( I ) ) THEN
00314 VMAX = STMP( I ) / SIN( I )
00315 ELSE
00316 VMAX = ONE
00317 END IF
00318 IF( VMAX.GT.RMAX( 3 ) ) THEN
00319 RMAX( 3 ) = VMAX
00320 IF( NINFO( 3 ).EQ.0 )
00321 $ LMAX( 3 ) = KNT
00322 END IF
00323 140 CONTINUE
00324
00325
00326
00327
00328 DO 150 I = 1, N
00329 IF( SEPIN( I ).LE.V .AND. SEPTMP( I ).LE.V ) THEN
00330 VMAX = ONE
00331 ELSE IF( EPS*SEPIN( I ).GT.SEPTMP( I ) ) THEN
00332 VMAX = ONE / EPS
00333 ELSE IF( SEPIN( I ).GT.SEPTMP( I ) ) THEN
00334 VMAX = SEPIN( I ) / SEPTMP( I )
00335 ELSE IF( SEPIN( I ).LT.EPS*SEPTMP( I ) ) THEN
00336 VMAX = ONE / EPS
00337 ELSE IF( SEPIN( I ).LT.SEPTMP( I ) ) THEN
00338 VMAX = SEPTMP( I ) / SEPIN( I )
00339 ELSE
00340 VMAX = ONE
00341 END IF
00342 IF( VMAX.GT.RMAX( 3 ) ) THEN
00343 RMAX( 3 ) = VMAX
00344 IF( NINFO( 3 ).EQ.0 )
00345 $ LMAX( 3 ) = KNT
00346 END IF
00347 150 CONTINUE
00348
00349
00350
00351 VMAX = ZERO
00352 DUM( 1 ) = -ONE
00353 CALL DCOPY( N, DUM, 0, STMP, 1 )
00354 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00355 CALL ZTRSNA( 'E', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00356 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00357 IF( INFO.NE.0 ) THEN
00358 LMAX( 3 ) = KNT
00359 NINFO( 3 ) = NINFO( 3 ) + 1
00360 GO TO 260
00361 END IF
00362 DO 160 I = 1, N
00363 IF( STMP( I ).NE.S( I ) )
00364 $ VMAX = ONE / EPS
00365 IF( SEPTMP( I ).NE.DUM( 1 ) )
00366 $ VMAX = ONE / EPS
00367 160 CONTINUE
00368
00369
00370
00371 CALL DCOPY( N, DUM, 0, STMP, 1 )
00372 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00373 CALL ZTRSNA( 'V', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00374 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00375 IF( INFO.NE.0 ) THEN
00376 LMAX( 3 ) = KNT
00377 NINFO( 3 ) = NINFO( 3 ) + 1
00378 GO TO 260
00379 END IF
00380 DO 170 I = 1, N
00381 IF( STMP( I ).NE.DUM( 1 ) )
00382 $ VMAX = ONE / EPS
00383 IF( SEPTMP( I ).NE.SEP( I ) )
00384 $ VMAX = ONE / EPS
00385 170 CONTINUE
00386
00387
00388
00389 DO 180 I = 1, N
00390 SELECT( I ) = .TRUE.
00391 180 CONTINUE
00392 CALL DCOPY( N, DUM, 0, STMP, 1 )
00393 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00394 CALL ZTRSNA( 'B', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00395 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00396 IF( INFO.NE.0 ) THEN
00397 LMAX( 3 ) = KNT
00398 NINFO( 3 ) = NINFO( 3 ) + 1
00399 GO TO 260
00400 END IF
00401 DO 190 I = 1, N
00402 IF( SEPTMP( I ).NE.SEP( I ) )
00403 $ VMAX = ONE / EPS
00404 IF( STMP( I ).NE.S( I ) )
00405 $ VMAX = ONE / EPS
00406 190 CONTINUE
00407
00408
00409
00410 CALL DCOPY( N, DUM, 0, STMP, 1 )
00411 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00412 CALL ZTRSNA( 'E', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00413 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00414 IF( INFO.NE.0 ) THEN
00415 LMAX( 3 ) = KNT
00416 NINFO( 3 ) = NINFO( 3 ) + 1
00417 GO TO 260
00418 END IF
00419 DO 200 I = 1, N
00420 IF( STMP( I ).NE.S( I ) )
00421 $ VMAX = ONE / EPS
00422 IF( SEPTMP( I ).NE.DUM( 1 ) )
00423 $ VMAX = ONE / EPS
00424 200 CONTINUE
00425
00426
00427
00428 CALL DCOPY( N, DUM, 0, STMP, 1 )
00429 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00430 CALL ZTRSNA( 'V', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00431 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00432 IF( INFO.NE.0 ) THEN
00433 LMAX( 3 ) = KNT
00434 NINFO( 3 ) = NINFO( 3 ) + 1
00435 GO TO 260
00436 END IF
00437 DO 210 I = 1, N
00438 IF( STMP( I ).NE.DUM( 1 ) )
00439 $ VMAX = ONE / EPS
00440 IF( SEPTMP( I ).NE.SEP( I ) )
00441 $ VMAX = ONE / EPS
00442 210 CONTINUE
00443 IF( VMAX.GT.RMAX( 1 ) ) THEN
00444 RMAX( 1 ) = VMAX
00445 IF( NINFO( 1 ).EQ.0 )
00446 $ LMAX( 1 ) = KNT
00447 END IF
00448
00449
00450
00451 DO 220 I = 1, N
00452 SELECT( I ) = .FALSE.
00453 220 CONTINUE
00454 ICMP = 0
00455 IF( N.GT.1 ) THEN
00456 ICMP = 1
00457 LCMP( 1 ) = 2
00458 SELECT( 2 ) = .TRUE.
00459 CALL ZCOPY( N, RE( 1, 2 ), 1, RE( 1, 1 ), 1 )
00460 CALL ZCOPY( N, LE( 1, 2 ), 1, LE( 1, 1 ), 1 )
00461 END IF
00462 IF( N.GT.3 ) THEN
00463 ICMP = 2
00464 LCMP( 2 ) = N - 1
00465 SELECT( N-1 ) = .TRUE.
00466 CALL ZCOPY( N, RE( 1, N-1 ), 1, RE( 1, 2 ), 1 )
00467 CALL ZCOPY( N, LE( 1, N-1 ), 1, LE( 1, 2 ), 1 )
00468 END IF
00469
00470
00471
00472 CALL DCOPY( ICMP, DUM, 0, STMP, 1 )
00473 CALL DCOPY( ICMP, DUM, 0, SEPTMP, 1 )
00474 CALL ZTRSNA( 'B', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00475 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00476 IF( INFO.NE.0 ) THEN
00477 LMAX( 3 ) = KNT
00478 NINFO( 3 ) = NINFO( 3 ) + 1
00479 GO TO 260
00480 END IF
00481 DO 230 I = 1, ICMP
00482 J = LCMP( I )
00483 IF( SEPTMP( I ).NE.SEP( J ) )
00484 $ VMAX = ONE / EPS
00485 IF( STMP( I ).NE.S( J ) )
00486 $ VMAX = ONE / EPS
00487 230 CONTINUE
00488
00489
00490
00491 CALL DCOPY( ICMP, DUM, 0, STMP, 1 )
00492 CALL DCOPY( ICMP, DUM, 0, SEPTMP, 1 )
00493 CALL ZTRSNA( 'E', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00494 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00495 IF( INFO.NE.0 ) THEN
00496 LMAX( 3 ) = KNT
00497 NINFO( 3 ) = NINFO( 3 ) + 1
00498 GO TO 260
00499 END IF
00500 DO 240 I = 1, ICMP
00501 J = LCMP( I )
00502 IF( STMP( I ).NE.S( J ) )
00503 $ VMAX = ONE / EPS
00504 IF( SEPTMP( I ).NE.DUM( 1 ) )
00505 $ VMAX = ONE / EPS
00506 240 CONTINUE
00507
00508
00509
00510 CALL DCOPY( ICMP, DUM, 0, STMP, 1 )
00511 CALL DCOPY( ICMP, DUM, 0, SEPTMP, 1 )
00512 CALL ZTRSNA( 'V', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
00513 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
00514 IF( INFO.NE.0 ) THEN
00515 LMAX( 3 ) = KNT
00516 NINFO( 3 ) = NINFO( 3 ) + 1
00517 GO TO 260
00518 END IF
00519 DO 250 I = 1, ICMP
00520 J = LCMP( I )
00521 IF( STMP( I ).NE.DUM( 1 ) )
00522 $ VMAX = ONE / EPS
00523 IF( SEPTMP( I ).NE.SEP( J ) )
00524 $ VMAX = ONE / EPS
00525 250 CONTINUE
00526 IF( VMAX.GT.RMAX( 1 ) ) THEN
00527 RMAX( 1 ) = VMAX
00528 IF( NINFO( 1 ).EQ.0 )
00529 $ LMAX( 1 ) = KNT
00530 END IF
00531 260 CONTINUE
00532 GO TO 10
00533
00534
00535
00536 END