00001 SUBROUTINE DDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
00002 $ ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE,
00003 $ RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK,
00004 $ IWORK, LIWORK, RESULT, BWORK, INFO )
00005
00006
00007
00008
00009
00010
00011 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
00012 $ NSIZE
00013 DOUBLE PRECISION THRESH
00014
00015
00016 LOGICAL BWORK( * )
00017 INTEGER IWORK( * )
00018 DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
00019 $ ALPHAR( * ), B( LDA, * ), BETA( * ),
00020 $ BI( LDA, * ), DIF( * ), DIFTRU( * ), DTRU( * ),
00021 $ LSCALE( * ), RESULT( 4 ), RSCALE( * ), S( * ),
00022 $ VL( LDA, * ), VR( LDA, * ), WORK( * )
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 DOUBLE PRECISION ZERO, ONE, TEN, TNTH
00189 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
00190 $ TNTH = 1.0D-1 )
00191
00192
00193 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
00194 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
00195 DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
00196 $ ULP, ULPINV
00197
00198
00199 DOUBLE PRECISION WEIGHT( 5 )
00200
00201
00202 INTEGER ILAENV
00203 DOUBLE PRECISION DLAMCH, DLANGE
00204 EXTERNAL ILAENV, DLAMCH, DLANGE
00205
00206
00207 EXTERNAL ALASVM, DGET52, DGGEVX, DLACPY, DLATM6, XERBLA
00208
00209
00210 INTRINSIC ABS, MAX, SQRT
00211
00212
00213
00214
00215
00216 INFO = 0
00217
00218 NMAX = 5
00219
00220 IF( NSIZE.LT.0 ) THEN
00221 INFO = -1
00222 ELSE IF( THRESH.LT.ZERO ) THEN
00223 INFO = -2
00224 ELSE IF( NIN.LE.0 ) THEN
00225 INFO = -3
00226 ELSE IF( NOUT.LE.0 ) THEN
00227 INFO = -4
00228 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
00229 INFO = -6
00230 ELSE IF( LIWORK.LT.NMAX+6 ) THEN
00231 INFO = -26
00232 END IF
00233
00234
00235
00236
00237
00238
00239
00240
00241 MINWRK = 1
00242 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
00243 MINWRK = 2*NMAX*NMAX + 12*NMAX + 16
00244 MAXWRK = 6*NMAX + NMAX*ILAENV( 1, 'DGEQRF', ' ', NMAX, 1, NMAX,
00245 $ 0 )
00246 MAXWRK = MAX( MAXWRK, 2*NMAX*NMAX+12*NMAX+16 )
00247 WORK( 1 ) = MAXWRK
00248 END IF
00249
00250 IF( LWORK.LT.MINWRK )
00251 $ INFO = -24
00252
00253 IF( INFO.NE.0 ) THEN
00254 CALL XERBLA( 'DDRGVX', -INFO )
00255 RETURN
00256 END IF
00257
00258 N = 5
00259 ULP = DLAMCH( 'P' )
00260 ULPINV = ONE / ULP
00261 THRSH2 = TEN*THRESH
00262 NERRS = 0
00263 NPTKNT = 0
00264 NTESTT = 0
00265
00266 IF( NSIZE.EQ.0 )
00267 $ GO TO 90
00268
00269
00270
00271 WEIGHT( 1 ) = SQRT( SQRT( ULP ) )
00272 WEIGHT( 2 ) = TNTH
00273 WEIGHT( 3 ) = ONE
00274 WEIGHT( 4 ) = ONE / WEIGHT( 2 )
00275 WEIGHT( 5 ) = ONE / WEIGHT( 1 )
00276
00277 DO 80 IPTYPE = 1, 2
00278 DO 70 IWA = 1, 5
00279 DO 60 IWB = 1, 5
00280 DO 50 IWX = 1, 5
00281 DO 40 IWY = 1, 5
00282
00283
00284
00285 CALL DLATM6( IPTYPE, 5, A, LDA, B, VR, LDA, VL,
00286 $ LDA, WEIGHT( IWA ), WEIGHT( IWB ),
00287 $ WEIGHT( IWX ), WEIGHT( IWY ), DTRU,
00288 $ DIFTRU )
00289
00290
00291
00292
00293
00294 CALL DLACPY( 'F', N, N, A, LDA, AI, LDA )
00295 CALL DLACPY( 'F', N, N, B, LDA, BI, LDA )
00296
00297 CALL DGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI,
00298 $ LDA, ALPHAR, ALPHAI, BETA, VL, LDA,
00299 $ VR, LDA, ILO, IHI, LSCALE, RSCALE,
00300 $ ANORM, BNORM, S, DIF, WORK, LWORK,
00301 $ IWORK, BWORK, LINFO )
00302 IF( LINFO.NE.0 ) THEN
00303 RESULT( 1 ) = ULPINV
00304 WRITE( NOUT, FMT = 9999 )'DGGEVX', LINFO, N,
00305 $ IPTYPE
00306 GO TO 30
00307 END IF
00308
00309
00310
00311 CALL DLACPY( 'Full', N, N, AI, LDA, WORK, N )
00312 CALL DLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ),
00313 $ N )
00314 ABNORM = DLANGE( 'Fro', N, 2*N, WORK, N, WORK )
00315
00316
00317
00318 RESULT( 1 ) = ZERO
00319 CALL DGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA,
00320 $ ALPHAR, ALPHAI, BETA, WORK,
00321 $ RESULT( 1 ) )
00322 IF( RESULT( 2 ).GT.THRESH ) THEN
00323 WRITE( NOUT, FMT = 9998 )'Left', 'DGGEVX',
00324 $ RESULT( 2 ), N, IPTYPE, IWA, IWB, IWX, IWY
00325 END IF
00326
00327 RESULT( 2 ) = ZERO
00328 CALL DGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA,
00329 $ ALPHAR, ALPHAI, BETA, WORK,
00330 $ RESULT( 2 ) )
00331 IF( RESULT( 3 ).GT.THRESH ) THEN
00332 WRITE( NOUT, FMT = 9998 )'Right', 'DGGEVX',
00333 $ RESULT( 3 ), N, IPTYPE, IWA, IWB, IWX, IWY
00334 END IF
00335
00336
00337
00338 RESULT( 3 ) = ZERO
00339 DO 10 I = 1, N
00340 IF( S( I ).EQ.ZERO ) THEN
00341 IF( DTRU( I ).GT.ABNORM*ULP )
00342 $ RESULT( 3 ) = ULPINV
00343 ELSE IF( DTRU( I ).EQ.ZERO ) THEN
00344 IF( S( I ).GT.ABNORM*ULP )
00345 $ RESULT( 3 ) = ULPINV
00346 ELSE
00347 WORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ),
00348 $ ABS( S( I ) / DTRU( I ) ) )
00349 RESULT( 3 ) = MAX( RESULT( 3 ), WORK( I ) )
00350 END IF
00351 10 CONTINUE
00352
00353
00354
00355 RESULT( 4 ) = ZERO
00356 IF( DIF( 1 ).EQ.ZERO ) THEN
00357 IF( DIFTRU( 1 ).GT.ABNORM*ULP )
00358 $ RESULT( 4 ) = ULPINV
00359 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
00360 IF( DIF( 1 ).GT.ABNORM*ULP )
00361 $ RESULT( 4 ) = ULPINV
00362 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
00363 IF( DIFTRU( 5 ).GT.ABNORM*ULP )
00364 $ RESULT( 4 ) = ULPINV
00365 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
00366 IF( DIF( 5 ).GT.ABNORM*ULP )
00367 $ RESULT( 4 ) = ULPINV
00368 ELSE
00369 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
00370 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
00371 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
00372 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
00373 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
00374 END IF
00375
00376 NTESTT = NTESTT + 4
00377
00378
00379
00380 DO 20 J = 1, 4
00381 IF( ( RESULT( J ).GE.THRSH2 .AND. J.GE.4 ) .OR.
00382 $ ( RESULT( J ).GE.THRESH .AND. J.LE.3 ) )
00383 $ THEN
00384
00385
00386
00387
00388 IF( NERRS.EQ.0 ) THEN
00389 WRITE( NOUT, FMT = 9997 )'DXV'
00390
00391
00392
00393
00394
00395 WRITE( NOUT, FMT = 9995 )
00396 WRITE( NOUT, FMT = 9994 )
00397 WRITE( NOUT, FMT = 9993 )
00398
00399
00400
00401 WRITE( NOUT, FMT = 9992 )'''',
00402 $ 'transpose', ''''
00403
00404 END IF
00405 NERRS = NERRS + 1
00406 IF( RESULT( J ).LT.10000.0D0 ) THEN
00407 WRITE( NOUT, FMT = 9991 )IPTYPE, IWA,
00408 $ IWB, IWX, IWY, J, RESULT( J )
00409 ELSE
00410 WRITE( NOUT, FMT = 9990 )IPTYPE, IWA,
00411 $ IWB, IWX, IWY, J, RESULT( J )
00412 END IF
00413 END IF
00414 20 CONTINUE
00415
00416 30 CONTINUE
00417
00418 40 CONTINUE
00419 50 CONTINUE
00420 60 CONTINUE
00421 70 CONTINUE
00422 80 CONTINUE
00423
00424 GO TO 150
00425
00426 90 CONTINUE
00427
00428
00429
00430
00431 READ( NIN, FMT = *, END = 150 )N
00432 IF( N.EQ.0 )
00433 $ GO TO 150
00434 DO 100 I = 1, N
00435 READ( NIN, FMT = * )( A( I, J ), J = 1, N )
00436 100 CONTINUE
00437 DO 110 I = 1, N
00438 READ( NIN, FMT = * )( B( I, J ), J = 1, N )
00439 110 CONTINUE
00440 READ( NIN, FMT = * )( DTRU( I ), I = 1, N )
00441 READ( NIN, FMT = * )( DIFTRU( I ), I = 1, N )
00442
00443 NPTKNT = NPTKNT + 1
00444
00445
00446
00447
00448
00449 CALL DLACPY( 'F', N, N, A, LDA, AI, LDA )
00450 CALL DLACPY( 'F', N, N, B, LDA, BI, LDA )
00451
00452 CALL DGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, LDA, ALPHAR,
00453 $ ALPHAI, BETA, VL, LDA, VR, LDA, ILO, IHI, LSCALE,
00454 $ RSCALE, ANORM, BNORM, S, DIF, WORK, LWORK, IWORK,
00455 $ BWORK, LINFO )
00456
00457 IF( LINFO.NE.0 ) THEN
00458 RESULT( 1 ) = ULPINV
00459 WRITE( NOUT, FMT = 9987 )'DGGEVX', LINFO, N, NPTKNT
00460 GO TO 140
00461 END IF
00462
00463
00464
00465 CALL DLACPY( 'Full', N, N, AI, LDA, WORK, N )
00466 CALL DLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), N )
00467 ABNORM = DLANGE( 'Fro', N, 2*N, WORK, N, WORK )
00468
00469
00470
00471 RESULT( 1 ) = ZERO
00472 CALL DGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHAR, ALPHAI,
00473 $ BETA, WORK, RESULT( 1 ) )
00474 IF( RESULT( 2 ).GT.THRESH ) THEN
00475 WRITE( NOUT, FMT = 9986 )'Left', 'DGGEVX', RESULT( 2 ), N,
00476 $ NPTKNT
00477 END IF
00478
00479 RESULT( 2 ) = ZERO
00480 CALL DGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHAR, ALPHAI,
00481 $ BETA, WORK, RESULT( 2 ) )
00482 IF( RESULT( 3 ).GT.THRESH ) THEN
00483 WRITE( NOUT, FMT = 9986 )'Right', 'DGGEVX', RESULT( 3 ), N,
00484 $ NPTKNT
00485 END IF
00486
00487
00488
00489 RESULT( 3 ) = ZERO
00490 DO 120 I = 1, N
00491 IF( S( I ).EQ.ZERO ) THEN
00492 IF( DTRU( I ).GT.ABNORM*ULP )
00493 $ RESULT( 3 ) = ULPINV
00494 ELSE IF( DTRU( I ).EQ.ZERO ) THEN
00495 IF( S( I ).GT.ABNORM*ULP )
00496 $ RESULT( 3 ) = ULPINV
00497 ELSE
00498 WORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ),
00499 $ ABS( S( I ) / DTRU( I ) ) )
00500 RESULT( 3 ) = MAX( RESULT( 3 ), WORK( I ) )
00501 END IF
00502 120 CONTINUE
00503
00504
00505
00506 RESULT( 4 ) = ZERO
00507 IF( DIF( 1 ).EQ.ZERO ) THEN
00508 IF( DIFTRU( 1 ).GT.ABNORM*ULP )
00509 $ RESULT( 4 ) = ULPINV
00510 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
00511 IF( DIF( 1 ).GT.ABNORM*ULP )
00512 $ RESULT( 4 ) = ULPINV
00513 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
00514 IF( DIFTRU( 5 ).GT.ABNORM*ULP )
00515 $ RESULT( 4 ) = ULPINV
00516 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
00517 IF( DIF( 5 ).GT.ABNORM*ULP )
00518 $ RESULT( 4 ) = ULPINV
00519 ELSE
00520 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
00521 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
00522 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
00523 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
00524 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
00525 END IF
00526
00527 NTESTT = NTESTT + 4
00528
00529
00530
00531 DO 130 J = 1, 4
00532 IF( RESULT( J ).GE.THRSH2 ) THEN
00533
00534
00535
00536
00537 IF( NERRS.EQ.0 ) THEN
00538 WRITE( NOUT, FMT = 9997 )'DXV'
00539
00540
00541
00542
00543
00544 WRITE( NOUT, FMT = 9996 )
00545
00546
00547
00548 WRITE( NOUT, FMT = 9992 )'''', 'transpose', ''''
00549
00550 END IF
00551 NERRS = NERRS + 1
00552 IF( RESULT( J ).LT.10000.0D0 ) THEN
00553 WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J )
00554 ELSE
00555 WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J )
00556 END IF
00557 END IF
00558 130 CONTINUE
00559
00560 140 CONTINUE
00561
00562 GO TO 90
00563 150 CONTINUE
00564
00565
00566
00567 CALL ALASVM( 'DXV', NOUT, NERRS, NTESTT, 0 )
00568
00569 WORK( 1 ) = MAXWRK
00570
00571 RETURN
00572
00573 9999 FORMAT( ' DDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00574 $ I6, ', JTYPE=', I6, ')' )
00575
00576 9998 FORMAT( ' DDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
00577 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
00578 $ 'N=', I6, ', JTYPE=', I6, ', IWA=', I5, ', IWB=', I5,
00579 $ ', IWX=', I5, ', IWY=', I5 )
00580
00581 9997 FORMAT( / 1X, A3, ' -- Real Expert Eigenvalue/vector',
00582 $ ' problem driver' )
00583
00584 9996 FORMAT( ' Input Example' )
00585
00586 9995 FORMAT( ' Matrix types: ', / )
00587
00588 9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
00589 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
00590 $ / ' YH and X are left and right eigenvectors. ', / )
00591
00592 9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
00593 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
00594 $ / ' YH and X are left and right eigenvectors. ', / )
00595
00596 9992 FORMAT( / ' Tests performed: ', / 4X,
00597 $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4X,
00598 $ ' r is a right eigenvector and ', A, ' means ', A, '.',
00599 $ / ' 1 = max | ( b A - a B )', A, ' l | / const.',
00600 $ / ' 2 = max | ( b A - a B ) r | / const.',
00601 $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
00602 $ ' over all eigenvalues', /
00603 $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
00604 $ ' over the 1st and 5th eigenvectors', / )
00605
00606 9991 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
00607 $ I2, ', IWY=', I2, ', result ', I2, ' is', 0P, F8.2 )
00608 9990 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
00609 $ I2, ', IWY=', I2, ', result ', I2, ' is', 1P, D10.3 )
00610 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00611 $ ' result ', I2, ' is', 0P, F8.2 )
00612 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00613 $ ' result ', I2, ' is', 1P, D10.3 )
00614 9987 FORMAT( ' DDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00615 $ I6, ', Input example #', I2, ')' )
00616
00617 9986 FORMAT( ' DDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
00618 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
00619 $ 'N=', I6, ', Input Example #', I2, ')' )
00620
00621
00622
00623
00624 END