00001 SUBROUTINE SGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
00002 $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
00003
00004
00005
00006
00007
00008
00009
00010 CHARACTER JOBVL, JOBVR
00011 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00012
00013
00014 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00015 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
00016 $ VR( LDVR, * ), WORK( * )
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
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 REAL ZERO, ONE
00225 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
00226
00227
00228 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
00229 CHARACTER CHTEMP
00230 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
00231 $ IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
00232 $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
00233 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
00234 $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
00235 $ SALFAI, SALFAR, SBETA, SCALE, TEMP
00236
00237
00238 LOGICAL LDUMMA( 1 )
00239
00240
00241 EXTERNAL SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLACPY,
00242 $ SLASCL, SLASET, SORGQR, SORMQR, STGEVC, XERBLA
00243
00244
00245 LOGICAL LSAME
00246 INTEGER ILAENV
00247 REAL SLAMCH, SLANGE
00248 EXTERNAL ILAENV, LSAME, SLAMCH, SLANGE
00249
00250
00251 INTRINSIC ABS, INT, MAX
00252
00253
00254
00255
00256
00257 IF( LSAME( JOBVL, 'N' ) ) THEN
00258 IJOBVL = 1
00259 ILVL = .FALSE.
00260 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
00261 IJOBVL = 2
00262 ILVL = .TRUE.
00263 ELSE
00264 IJOBVL = -1
00265 ILVL = .FALSE.
00266 END IF
00267
00268 IF( LSAME( JOBVR, 'N' ) ) THEN
00269 IJOBVR = 1
00270 ILVR = .FALSE.
00271 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
00272 IJOBVR = 2
00273 ILVR = .TRUE.
00274 ELSE
00275 IJOBVR = -1
00276 ILVR = .FALSE.
00277 END IF
00278 ILV = ILVL .OR. ILVR
00279
00280
00281
00282 LWKMIN = MAX( 8*N, 1 )
00283 LWKOPT = LWKMIN
00284 WORK( 1 ) = LWKOPT
00285 LQUERY = ( LWORK.EQ.-1 )
00286 INFO = 0
00287 IF( IJOBVL.LE.0 ) THEN
00288 INFO = -1
00289 ELSE IF( IJOBVR.LE.0 ) THEN
00290 INFO = -2
00291 ELSE IF( N.LT.0 ) THEN
00292 INFO = -3
00293 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00294 INFO = -5
00295 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00296 INFO = -7
00297 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
00298 INFO = -12
00299 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
00300 INFO = -14
00301 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00302 INFO = -16
00303 END IF
00304
00305 IF( INFO.EQ.0 ) THEN
00306 NB1 = ILAENV( 1, 'SGEQRF', ' ', N, N, -1, -1 )
00307 NB2 = ILAENV( 1, 'SORMQR', ' ', N, N, N, -1 )
00308 NB3 = ILAENV( 1, 'SORGQR', ' ', N, N, N, -1 )
00309 NB = MAX( NB1, NB2, NB3 )
00310 LOPT = 2*N + MAX( 6*N, N*(NB+1) )
00311 WORK( 1 ) = LOPT
00312 END IF
00313
00314 IF( INFO.NE.0 ) THEN
00315 CALL XERBLA( 'SGEGV ', -INFO )
00316 RETURN
00317 ELSE IF( LQUERY ) THEN
00318 RETURN
00319 END IF
00320
00321
00322
00323 IF( N.EQ.0 )
00324 $ RETURN
00325
00326
00327
00328 EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
00329 SAFMIN = SLAMCH( 'S' )
00330 SAFMIN = SAFMIN + SAFMIN
00331 SAFMAX = ONE / SAFMIN
00332 ONEPLS = ONE + ( 4*EPS )
00333
00334
00335
00336 ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
00337 ANRM1 = ANRM
00338 ANRM2 = ONE
00339 IF( ANRM.LT.ONE ) THEN
00340 IF( SAFMAX*ANRM.LT.ONE ) THEN
00341 ANRM1 = SAFMIN
00342 ANRM2 = SAFMAX*ANRM
00343 END IF
00344 END IF
00345
00346 IF( ANRM.GT.ZERO ) THEN
00347 CALL SLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
00348 IF( IINFO.NE.0 ) THEN
00349 INFO = N + 10
00350 RETURN
00351 END IF
00352 END IF
00353
00354
00355
00356 BNRM = SLANGE( 'M', N, N, B, LDB, WORK )
00357 BNRM1 = BNRM
00358 BNRM2 = ONE
00359 IF( BNRM.LT.ONE ) THEN
00360 IF( SAFMAX*BNRM.LT.ONE ) THEN
00361 BNRM1 = SAFMIN
00362 BNRM2 = SAFMAX*BNRM
00363 END IF
00364 END IF
00365
00366 IF( BNRM.GT.ZERO ) THEN
00367 CALL SLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
00368 IF( IINFO.NE.0 ) THEN
00369 INFO = N + 10
00370 RETURN
00371 END IF
00372 END IF
00373
00374
00375
00376
00377
00378 ILEFT = 1
00379 IRIGHT = N + 1
00380 IWORK = IRIGHT + N
00381 CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
00382 $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
00383 IF( IINFO.NE.0 ) THEN
00384 INFO = N + 1
00385 GO TO 120
00386 END IF
00387
00388
00389
00390
00391
00392 IROWS = IHI + 1 - ILO
00393 IF( ILV ) THEN
00394 ICOLS = N + 1 - ILO
00395 ELSE
00396 ICOLS = IROWS
00397 END IF
00398 ITAU = IWORK
00399 IWORK = ITAU + IROWS
00400 CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00401 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
00402 IF( IINFO.GE.0 )
00403 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00404 IF( IINFO.NE.0 ) THEN
00405 INFO = N + 2
00406 GO TO 120
00407 END IF
00408
00409 CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00410 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00411 $ LWORK+1-IWORK, IINFO )
00412 IF( IINFO.GE.0 )
00413 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00414 IF( IINFO.NE.0 ) THEN
00415 INFO = N + 3
00416 GO TO 120
00417 END IF
00418
00419 IF( ILVL ) THEN
00420 CALL SLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
00421 CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00422 $ VL( ILO+1, ILO ), LDVL )
00423 CALL SORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
00424 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00425 $ IINFO )
00426 IF( IINFO.GE.0 )
00427 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00428 IF( IINFO.NE.0 ) THEN
00429 INFO = N + 4
00430 GO TO 120
00431 END IF
00432 END IF
00433
00434 IF( ILVR )
00435 $ CALL SLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
00436
00437
00438
00439 IF( ILV ) THEN
00440
00441
00442
00443 CALL SGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
00444 $ LDVL, VR, LDVR, IINFO )
00445 ELSE
00446 CALL SGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
00447 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
00448 END IF
00449 IF( IINFO.NE.0 ) THEN
00450 INFO = N + 5
00451 GO TO 120
00452 END IF
00453
00454
00455
00456
00457
00458 IWORK = ITAU
00459 IF( ILV ) THEN
00460 CHTEMP = 'S'
00461 ELSE
00462 CHTEMP = 'E'
00463 END IF
00464 CALL SHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
00465 $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
00466 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
00467 IF( IINFO.GE.0 )
00468 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00469 IF( IINFO.NE.0 ) THEN
00470 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00471 INFO = IINFO
00472 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00473 INFO = IINFO - N
00474 ELSE
00475 INFO = N + 6
00476 END IF
00477 GO TO 120
00478 END IF
00479
00480 IF( ILV ) THEN
00481
00482
00483
00484 IF( ILVL ) THEN
00485 IF( ILVR ) THEN
00486 CHTEMP = 'B'
00487 ELSE
00488 CHTEMP = 'L'
00489 END IF
00490 ELSE
00491 CHTEMP = 'R'
00492 END IF
00493
00494 CALL STGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
00495 $ VR, LDVR, N, IN, WORK( IWORK ), IINFO )
00496 IF( IINFO.NE.0 ) THEN
00497 INFO = N + 7
00498 GO TO 120
00499 END IF
00500
00501
00502
00503 IF( ILVL ) THEN
00504 CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
00505 $ WORK( IRIGHT ), N, VL, LDVL, IINFO )
00506 IF( IINFO.NE.0 ) THEN
00507 INFO = N + 8
00508 GO TO 120
00509 END IF
00510 DO 50 JC = 1, N
00511 IF( ALPHAI( JC ).LT.ZERO )
00512 $ GO TO 50
00513 TEMP = ZERO
00514 IF( ALPHAI( JC ).EQ.ZERO ) THEN
00515 DO 10 JR = 1, N
00516 TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
00517 10 CONTINUE
00518 ELSE
00519 DO 20 JR = 1, N
00520 TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
00521 $ ABS( VL( JR, JC+1 ) ) )
00522 20 CONTINUE
00523 END IF
00524 IF( TEMP.LT.SAFMIN )
00525 $ GO TO 50
00526 TEMP = ONE / TEMP
00527 IF( ALPHAI( JC ).EQ.ZERO ) THEN
00528 DO 30 JR = 1, N
00529 VL( JR, JC ) = VL( JR, JC )*TEMP
00530 30 CONTINUE
00531 ELSE
00532 DO 40 JR = 1, N
00533 VL( JR, JC ) = VL( JR, JC )*TEMP
00534 VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
00535 40 CONTINUE
00536 END IF
00537 50 CONTINUE
00538 END IF
00539 IF( ILVR ) THEN
00540 CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
00541 $ WORK( IRIGHT ), N, VR, LDVR, IINFO )
00542 IF( IINFO.NE.0 ) THEN
00543 INFO = N + 9
00544 GO TO 120
00545 END IF
00546 DO 100 JC = 1, N
00547 IF( ALPHAI( JC ).LT.ZERO )
00548 $ GO TO 100
00549 TEMP = ZERO
00550 IF( ALPHAI( JC ).EQ.ZERO ) THEN
00551 DO 60 JR = 1, N
00552 TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
00553 60 CONTINUE
00554 ELSE
00555 DO 70 JR = 1, N
00556 TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
00557 $ ABS( VR( JR, JC+1 ) ) )
00558 70 CONTINUE
00559 END IF
00560 IF( TEMP.LT.SAFMIN )
00561 $ GO TO 100
00562 TEMP = ONE / TEMP
00563 IF( ALPHAI( JC ).EQ.ZERO ) THEN
00564 DO 80 JR = 1, N
00565 VR( JR, JC ) = VR( JR, JC )*TEMP
00566 80 CONTINUE
00567 ELSE
00568 DO 90 JR = 1, N
00569 VR( JR, JC ) = VR( JR, JC )*TEMP
00570 VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
00571 90 CONTINUE
00572 END IF
00573 100 CONTINUE
00574 END IF
00575
00576
00577
00578 END IF
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588 DO 110 JC = 1, N
00589 ABSAR = ABS( ALPHAR( JC ) )
00590 ABSAI = ABS( ALPHAI( JC ) )
00591 ABSB = ABS( BETA( JC ) )
00592 SALFAR = ANRM*ALPHAR( JC )
00593 SALFAI = ANRM*ALPHAI( JC )
00594 SBETA = BNRM*BETA( JC )
00595 ILIMIT = .FALSE.
00596 SCALE = ONE
00597
00598
00599
00600 IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
00601 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
00602 ILIMIT = .TRUE.
00603 SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
00604 $ MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
00605
00606 ELSE IF( SALFAI.EQ.ZERO ) THEN
00607
00608
00609
00610
00611 IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
00612 ALPHAI( JC-1 ) = ZERO
00613 ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
00614 ALPHAI( JC+1 ) = ZERO
00615 END IF
00616 END IF
00617
00618
00619
00620 IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
00621 $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
00622 ILIMIT = .TRUE.
00623 SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
00624 $ MAX( ONEPLS*SAFMIN, ANRM2*ABSAR ) )
00625 END IF
00626
00627
00628
00629 IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
00630 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
00631 ILIMIT = .TRUE.
00632 SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / BNRM1 ) /
00633 $ MAX( ONEPLS*SAFMIN, BNRM2*ABSB ) )
00634 END IF
00635
00636
00637
00638 IF( ILIMIT ) THEN
00639 TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
00640 $ ABS( SBETA ) )
00641 IF( TEMP.GT.ONE )
00642 $ SCALE = SCALE / TEMP
00643 IF( SCALE.LT.ONE )
00644 $ ILIMIT = .FALSE.
00645 END IF
00646
00647
00648
00649 IF( ILIMIT ) THEN
00650 SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
00651 SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
00652 SBETA = ( SCALE*BETA( JC ) )*BNRM
00653 END IF
00654 ALPHAR( JC ) = SALFAR
00655 ALPHAI( JC ) = SALFAI
00656 BETA( JC ) = SBETA
00657 110 CONTINUE
00658
00659 120 CONTINUE
00660 WORK( 1 ) = LWKOPT
00661
00662 RETURN
00663
00664
00665
00666 END