00001 SUBROUTINE DGET37( 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, IFND, INFO, ISCL, J, KMIN, M, N
00065 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
00066 $ VIMIN, VMAX, VMUL, VRMIN
00067
00068
00069 LOGICAL SELECT( LDT )
00070 INTEGER IWORK( 2*LDT ), LCMP( 3 )
00071 DOUBLE PRECISION DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
00072 $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
00073 $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
00074 $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
00075 $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
00076 $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
00077 $ WRTMP( LDT )
00078
00079
00080 DOUBLE PRECISION DLAMCH, DLANGE
00081 EXTERNAL DLAMCH, DLANGE
00082
00083
00084 EXTERNAL DCOPY, DGEHRD, DHSEQR, DLABAD, DLACPY, DSCAL,
00085 $ DTREVC, DTRSNA
00086
00087
00088 INTRINSIC DBLE, 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
00111 VAL( 1 ) = SQRT( SMLNUM )
00112 VAL( 2 ) = ONE
00113 VAL( 3 ) = SQRT( BIGNUM )
00114
00115
00116
00117
00118
00119 10 CONTINUE
00120 READ( NIN, FMT = * )N
00121 IF( N.EQ.0 )
00122 $ RETURN
00123 DO 20 I = 1, N
00124 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
00125 20 CONTINUE
00126 DO 30 I = 1, N
00127 READ( NIN, FMT = * )WRIN( I ), WIIN( I ), SIN( I ), SEPIN( I )
00128 30 CONTINUE
00129 TNRM = DLANGE( 'M', N, N, TMP, LDT, WORK )
00130
00131
00132
00133 DO 240 ISCL = 1, 3
00134
00135
00136
00137 KNT = KNT + 1
00138 CALL DLACPY( 'F', N, N, TMP, LDT, T, LDT )
00139 VMUL = VAL( ISCL )
00140 DO 40 I = 1, N
00141 CALL DSCAL( N, VMUL, T( 1, I ), 1 )
00142 40 CONTINUE
00143 IF( TNRM.EQ.ZERO )
00144 $ VMUL = ONE
00145
00146
00147
00148 CALL DGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
00149 $ INFO )
00150 IF( INFO.NE.0 ) THEN
00151 LMAX( 1 ) = KNT
00152 NINFO( 1 ) = NINFO( 1 ) + 1
00153 GO TO 240
00154 END IF
00155 DO 60 J = 1, N - 2
00156 DO 50 I = J + 2, N
00157 T( I, J ) = ZERO
00158 50 CONTINUE
00159 60 CONTINUE
00160
00161
00162
00163 CALL DHSEQR( 'S', 'N', N, 1, N, T, LDT, WR, WI, DUM, 1, WORK,
00164 $ LWORK, INFO )
00165 IF( INFO.NE.0 ) THEN
00166 LMAX( 2 ) = KNT
00167 NINFO( 2 ) = NINFO( 2 ) + 1
00168 GO TO 240
00169 END IF
00170
00171
00172
00173 CALL DTREVC( 'Both', 'All', SELECT, N, T, LDT, LE, LDT, RE,
00174 $ LDT, N, M, WORK, INFO )
00175
00176
00177
00178 CALL DTRSNA( 'Both', 'All', SELECT, N, T, LDT, LE, LDT, RE,
00179 $ LDT, S, SEP, N, M, WORK, N, IWORK, INFO )
00180 IF( INFO.NE.0 ) THEN
00181 LMAX( 3 ) = KNT
00182 NINFO( 3 ) = NINFO( 3 ) + 1
00183 GO TO 240
00184 END IF
00185
00186
00187
00188
00189 CALL DCOPY( N, WR, 1, WRTMP, 1 )
00190 CALL DCOPY( N, WI, 1, WITMP, 1 )
00191 CALL DCOPY( N, S, 1, STMP, 1 )
00192 CALL DCOPY( N, SEP, 1, SEPTMP, 1 )
00193 CALL DSCAL( N, ONE / VMUL, SEPTMP, 1 )
00194 DO 80 I = 1, N - 1
00195 KMIN = I
00196 VRMIN = WRTMP( I )
00197 VIMIN = WITMP( I )
00198 DO 70 J = I + 1, N
00199 IF( WRTMP( J ).LT.VRMIN ) THEN
00200 KMIN = J
00201 VRMIN = WRTMP( J )
00202 VIMIN = WITMP( J )
00203 END IF
00204 70 CONTINUE
00205 WRTMP( KMIN ) = WRTMP( I )
00206 WITMP( KMIN ) = WITMP( I )
00207 WRTMP( I ) = VRMIN
00208 WITMP( I ) = VIMIN
00209 VRMIN = STMP( KMIN )
00210 STMP( KMIN ) = STMP( I )
00211 STMP( I ) = VRMIN
00212 VRMIN = SEPTMP( KMIN )
00213 SEPTMP( KMIN ) = SEPTMP( I )
00214 SEPTMP( I ) = VRMIN
00215 80 CONTINUE
00216
00217
00218
00219
00220 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM )
00221 IF( TNRM.EQ.ZERO )
00222 $ V = ONE
00223 DO 90 I = 1, N
00224 IF( V.GT.SEPTMP( I ) ) THEN
00225 TOL = ONE
00226 ELSE
00227 TOL = V / SEPTMP( I )
00228 END IF
00229 IF( V.GT.SEPIN( I ) ) THEN
00230 TOLIN = ONE
00231 ELSE
00232 TOLIN = V / SEPIN( I )
00233 END IF
00234 TOL = MAX( TOL, SMLNUM / EPS )
00235 TOLIN = MAX( TOLIN, SMLNUM / EPS )
00236 IF( EPS*( SIN( I )-TOLIN ).GT.STMP( I )+TOL ) THEN
00237 VMAX = ONE / EPS
00238 ELSE IF( SIN( I )-TOLIN.GT.STMP( I )+TOL ) THEN
00239 VMAX = ( SIN( I )-TOLIN ) / ( STMP( I )+TOL )
00240 ELSE IF( SIN( I )+TOLIN.LT.EPS*( STMP( I )-TOL ) ) THEN
00241 VMAX = ONE / EPS
00242 ELSE IF( SIN( I )+TOLIN.LT.STMP( I )-TOL ) THEN
00243 VMAX = ( STMP( I )-TOL ) / ( SIN( I )+TOLIN )
00244 ELSE
00245 VMAX = ONE
00246 END IF
00247 IF( VMAX.GT.RMAX( 2 ) ) THEN
00248 RMAX( 2 ) = VMAX
00249 IF( NINFO( 2 ).EQ.0 )
00250 $ LMAX( 2 ) = KNT
00251 END IF
00252 90 CONTINUE
00253
00254
00255
00256
00257 DO 100 I = 1, N
00258 IF( V.GT.SEPTMP( I )*STMP( I ) ) THEN
00259 TOL = SEPTMP( I )
00260 ELSE
00261 TOL = V / STMP( I )
00262 END IF
00263 IF( V.GT.SEPIN( I )*SIN( I ) ) THEN
00264 TOLIN = SEPIN( I )
00265 ELSE
00266 TOLIN = V / SIN( I )
00267 END IF
00268 TOL = MAX( TOL, SMLNUM / EPS )
00269 TOLIN = MAX( TOLIN, SMLNUM / EPS )
00270 IF( EPS*( SEPIN( I )-TOLIN ).GT.SEPTMP( I )+TOL ) THEN
00271 VMAX = ONE / EPS
00272 ELSE IF( SEPIN( I )-TOLIN.GT.SEPTMP( I )+TOL ) THEN
00273 VMAX = ( SEPIN( I )-TOLIN ) / ( SEPTMP( I )+TOL )
00274 ELSE IF( SEPIN( I )+TOLIN.LT.EPS*( SEPTMP( I )-TOL ) ) THEN
00275 VMAX = ONE / EPS
00276 ELSE IF( SEPIN( I )+TOLIN.LT.SEPTMP( I )-TOL ) THEN
00277 VMAX = ( SEPTMP( I )-TOL ) / ( SEPIN( I )+TOLIN )
00278 ELSE
00279 VMAX = ONE
00280 END IF
00281 IF( VMAX.GT.RMAX( 2 ) ) THEN
00282 RMAX( 2 ) = VMAX
00283 IF( NINFO( 2 ).EQ.0 )
00284 $ LMAX( 2 ) = KNT
00285 END IF
00286 100 CONTINUE
00287
00288
00289
00290
00291 DO 110 I = 1, N
00292 IF( SIN( I ).LE.DBLE( 2*N )*EPS .AND. STMP( I ).LE.
00293 $ DBLE( 2*N )*EPS ) THEN
00294 VMAX = ONE
00295 ELSE IF( EPS*SIN( I ).GT.STMP( I ) ) THEN
00296 VMAX = ONE / EPS
00297 ELSE IF( SIN( I ).GT.STMP( I ) ) THEN
00298 VMAX = SIN( I ) / STMP( I )
00299 ELSE IF( SIN( I ).LT.EPS*STMP( I ) ) THEN
00300 VMAX = ONE / EPS
00301 ELSE IF( SIN( I ).LT.STMP( I ) ) THEN
00302 VMAX = STMP( I ) / SIN( I )
00303 ELSE
00304 VMAX = ONE
00305 END IF
00306 IF( VMAX.GT.RMAX( 3 ) ) THEN
00307 RMAX( 3 ) = VMAX
00308 IF( NINFO( 3 ).EQ.0 )
00309 $ LMAX( 3 ) = KNT
00310 END IF
00311 110 CONTINUE
00312
00313
00314
00315
00316 DO 120 I = 1, N
00317 IF( SEPIN( I ).LE.V .AND. SEPTMP( I ).LE.V ) THEN
00318 VMAX = ONE
00319 ELSE IF( EPS*SEPIN( I ).GT.SEPTMP( I ) ) THEN
00320 VMAX = ONE / EPS
00321 ELSE IF( SEPIN( I ).GT.SEPTMP( I ) ) THEN
00322 VMAX = SEPIN( I ) / SEPTMP( I )
00323 ELSE IF( SEPIN( I ).LT.EPS*SEPTMP( I ) ) THEN
00324 VMAX = ONE / EPS
00325 ELSE IF( SEPIN( I ).LT.SEPTMP( I ) ) THEN
00326 VMAX = SEPTMP( I ) / SEPIN( I )
00327 ELSE
00328 VMAX = ONE
00329 END IF
00330 IF( VMAX.GT.RMAX( 3 ) ) THEN
00331 RMAX( 3 ) = VMAX
00332 IF( NINFO( 3 ).EQ.0 )
00333 $ LMAX( 3 ) = KNT
00334 END IF
00335 120 CONTINUE
00336
00337
00338
00339 VMAX = ZERO
00340 DUM( 1 ) = -ONE
00341 CALL DCOPY( N, DUM, 0, STMP, 1 )
00342 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00343 CALL DTRSNA( 'Eigcond', 'All', SELECT, N, T, LDT, LE, LDT, RE,
00344 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
00345 IF( INFO.NE.0 ) THEN
00346 LMAX( 3 ) = KNT
00347 NINFO( 3 ) = NINFO( 3 ) + 1
00348 GO TO 240
00349 END IF
00350 DO 130 I = 1, N
00351 IF( STMP( I ).NE.S( I ) )
00352 $ VMAX = ONE / EPS
00353 IF( SEPTMP( I ).NE.DUM( 1 ) )
00354 $ VMAX = ONE / EPS
00355 130 CONTINUE
00356
00357
00358
00359 CALL DCOPY( N, DUM, 0, STMP, 1 )
00360 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00361 CALL DTRSNA( 'Veccond', 'All', SELECT, N, T, LDT, LE, LDT, RE,
00362 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
00363 IF( INFO.NE.0 ) THEN
00364 LMAX( 3 ) = KNT
00365 NINFO( 3 ) = NINFO( 3 ) + 1
00366 GO TO 240
00367 END IF
00368 DO 140 I = 1, N
00369 IF( STMP( I ).NE.DUM( 1 ) )
00370 $ VMAX = ONE / EPS
00371 IF( SEPTMP( I ).NE.SEP( I ) )
00372 $ VMAX = ONE / EPS
00373 140 CONTINUE
00374
00375
00376
00377 DO 150 I = 1, N
00378 SELECT( I ) = .TRUE.
00379 150 CONTINUE
00380 CALL DCOPY( N, DUM, 0, STMP, 1 )
00381 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00382 CALL DTRSNA( 'Bothcond', 'Some', SELECT, N, T, LDT, LE, LDT,
00383 $ RE, LDT, STMP, SEPTMP, N, M, WORK, N, IWORK,
00384 $ INFO )
00385 IF( INFO.NE.0 ) THEN
00386 LMAX( 3 ) = KNT
00387 NINFO( 3 ) = NINFO( 3 ) + 1
00388 GO TO 240
00389 END IF
00390 DO 160 I = 1, N
00391 IF( SEPTMP( I ).NE.SEP( I ) )
00392 $ VMAX = ONE / EPS
00393 IF( STMP( I ).NE.S( I ) )
00394 $ VMAX = ONE / EPS
00395 160 CONTINUE
00396
00397
00398
00399 CALL DCOPY( N, DUM, 0, STMP, 1 )
00400 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00401 CALL DTRSNA( 'Eigcond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
00402 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
00403 IF( INFO.NE.0 ) THEN
00404 LMAX( 3 ) = KNT
00405 NINFO( 3 ) = NINFO( 3 ) + 1
00406 GO TO 240
00407 END IF
00408 DO 170 I = 1, N
00409 IF( STMP( I ).NE.S( I ) )
00410 $ VMAX = ONE / EPS
00411 IF( SEPTMP( I ).NE.DUM( 1 ) )
00412 $ VMAX = ONE / EPS
00413 170 CONTINUE
00414
00415
00416
00417 CALL DCOPY( N, DUM, 0, STMP, 1 )
00418 CALL DCOPY( N, DUM, 0, SEPTMP, 1 )
00419 CALL DTRSNA( 'Veccond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
00420 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
00421 IF( INFO.NE.0 ) THEN
00422 LMAX( 3 ) = KNT
00423 NINFO( 3 ) = NINFO( 3 ) + 1
00424 GO TO 240
00425 END IF
00426 DO 180 I = 1, N
00427 IF( STMP( I ).NE.DUM( 1 ) )
00428 $ VMAX = ONE / EPS
00429 IF( SEPTMP( I ).NE.SEP( I ) )
00430 $ VMAX = ONE / EPS
00431 180 CONTINUE
00432 IF( VMAX.GT.RMAX( 1 ) ) THEN
00433 RMAX( 1 ) = VMAX
00434 IF( NINFO( 1 ).EQ.0 )
00435 $ LMAX( 1 ) = KNT
00436 END IF
00437
00438
00439
00440 IF( WI( 1 ).EQ.ZERO ) THEN
00441 LCMP( 1 ) = 1
00442 IFND = 0
00443 DO 190 I = 2, N
00444 IF( IFND.EQ.1 .OR. WI( I ).EQ.ZERO ) THEN
00445 SELECT( I ) = .FALSE.
00446 ELSE
00447 IFND = 1
00448 LCMP( 2 ) = I
00449 LCMP( 3 ) = I + 1
00450 CALL DCOPY( N, RE( 1, I ), 1, RE( 1, 2 ), 1 )
00451 CALL DCOPY( N, RE( 1, I+1 ), 1, RE( 1, 3 ), 1 )
00452 CALL DCOPY( N, LE( 1, I ), 1, LE( 1, 2 ), 1 )
00453 CALL DCOPY( N, LE( 1, I+1 ), 1, LE( 1, 3 ), 1 )
00454 END IF
00455 190 CONTINUE
00456 IF( IFND.EQ.0 ) THEN
00457 ICMP = 1
00458 ELSE
00459 ICMP = 3
00460 END IF
00461 ELSE
00462 LCMP( 1 ) = 1
00463 LCMP( 2 ) = 2
00464 IFND = 0
00465 DO 200 I = 3, N
00466 IF( IFND.EQ.1 .OR. WI( I ).NE.ZERO ) THEN
00467 SELECT( I ) = .FALSE.
00468 ELSE
00469 LCMP( 3 ) = I
00470 IFND = 1
00471 CALL DCOPY( N, RE( 1, I ), 1, RE( 1, 3 ), 1 )
00472 CALL DCOPY( N, LE( 1, I ), 1, LE( 1, 3 ), 1 )
00473 END IF
00474 200 CONTINUE
00475 IF( IFND.EQ.0 ) THEN
00476 ICMP = 2
00477 ELSE
00478 ICMP = 3
00479 END IF
00480 END IF
00481
00482
00483
00484 CALL DCOPY( ICMP, DUM, 0, STMP, 1 )
00485 CALL DCOPY( ICMP, DUM, 0, SEPTMP, 1 )
00486 CALL DTRSNA( 'Bothcond', 'Some', SELECT, N, T, LDT, LE, LDT,
00487 $ RE, LDT, STMP, SEPTMP, N, M, WORK, N, IWORK,
00488 $ INFO )
00489 IF( INFO.NE.0 ) THEN
00490 LMAX( 3 ) = KNT
00491 NINFO( 3 ) = NINFO( 3 ) + 1
00492 GO TO 240
00493 END IF
00494 DO 210 I = 1, ICMP
00495 J = LCMP( I )
00496 IF( SEPTMP( I ).NE.SEP( J ) )
00497 $ VMAX = ONE / EPS
00498 IF( STMP( I ).NE.S( J ) )
00499 $ VMAX = ONE / EPS
00500 210 CONTINUE
00501
00502
00503
00504 CALL DCOPY( ICMP, DUM, 0, STMP, 1 )
00505 CALL DCOPY( ICMP, DUM, 0, SEPTMP, 1 )
00506 CALL DTRSNA( 'Eigcond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
00507 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
00508 IF( INFO.NE.0 ) THEN
00509 LMAX( 3 ) = KNT
00510 NINFO( 3 ) = NINFO( 3 ) + 1
00511 GO TO 240
00512 END IF
00513 DO 220 I = 1, ICMP
00514 J = LCMP( I )
00515 IF( STMP( I ).NE.S( J ) )
00516 $ VMAX = ONE / EPS
00517 IF( SEPTMP( I ).NE.DUM( 1 ) )
00518 $ VMAX = ONE / EPS
00519 220 CONTINUE
00520
00521
00522
00523 CALL DCOPY( ICMP, DUM, 0, STMP, 1 )
00524 CALL DCOPY( ICMP, DUM, 0, SEPTMP, 1 )
00525 CALL DTRSNA( 'Veccond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
00526 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
00527 IF( INFO.NE.0 ) THEN
00528 LMAX( 3 ) = KNT
00529 NINFO( 3 ) = NINFO( 3 ) + 1
00530 GO TO 240
00531 END IF
00532 DO 230 I = 1, ICMP
00533 J = LCMP( I )
00534 IF( STMP( I ).NE.DUM( 1 ) )
00535 $ VMAX = ONE / EPS
00536 IF( SEPTMP( I ).NE.SEP( J ) )
00537 $ VMAX = ONE / EPS
00538 230 CONTINUE
00539 IF( VMAX.GT.RMAX( 1 ) ) THEN
00540 RMAX( 1 ) = VMAX
00541 IF( NINFO( 1 ).EQ.0 )
00542 $ LMAX( 1 ) = KNT
00543 END IF
00544 240 CONTINUE
00545 GO TO 10
00546
00547
00548
00549 END