I'm trying to solve a generalized eigenvalue problem. I have two matrices H and S such that:
HX=λSX
I need to find the eigenvalues λ. The matrices H and S are real, asymmetric, contains positive and negative decimal numbers and (can be) singular. I tried this on MATLAB by using the command eig(H,S), but the problem is that I am obtaining real and complex eigenvalues, while in the research paper that I am following the authors have got only real eigenvalues.
After reading about this problem I have got to know that for such matrices, softwares like MATLAB and many other use the QZ Algorithm to solve generalized eigenavalue problem. I am seeking answers to the following questions:
- Is there any criteria to determine whether the given matrices will have complex eigenvalues or not? and that the result obtained from MATLAB is correct or not? (though the values don't agree with research paper).
- In the research paper the authors have used RGG subroutine which is a part of Fortran library EISPACK. I went through its documentation and got to know that it also uses the QZ algorithm. So my question is can MATLAB and Fortran give different answers to the same problem even when both of them are using the same algorithm?
Here are the matrices:
H=[0,0,0,0,192,1917.04064,10332.51505,40092.51227,125681.1486,338350.2206,811892.8294,1779728.921,3625982.355,6953387.916,12670976.81,22100000,37132930.27,60353006.25,95276316.19,146559937.4,220274060.5,324208411.9,468219308,664618721.8;
0,0,0,0,192,1893.475124,10051.90014,38308.22391,117609.6433,309187.9535,722364.2569,1537115.973,3030677.025,5606681.841,9824567.083,16426180.74,26355891.41,40770017.02,61031174.08,88683247.04,125403139.6,172926323.5,232944447.3,306974887;
0,0,0,0,192,1846.924351,9507.779872,34923.83828,102685.9158,256815.0258,566753.8301,1130503.089,2072244.588,3531888.862,5644710.273,8510616.36,12154554.65,16481872.51,21234784.98,25958063.51,29983201.26,32440241.94,32304892.86,28485362.66;
0,0,0,0,192,1778.534558,8732.953666,30281.20319,83090.50097,191433.6927,383369.9941,681577.7806,1089000.57,1571702.814,2043718.6,2360415.659,2327150.205,1728131.214,376286.8361,-1820993.208,-4791939.112,-8243585.409,-11638086.16,-14220529.21;
0,0,0,0,192,1689.989727,7773.203992,24831.36406,61515.88289,124689.9821,212051.9399,303583.89,356337.2627,308129.1939,94242.56924,-323226.8121,-920239.2345,-1583058.32,-2104105.201,-2214813.175,-1659107.74,-295993.1841,1796496.926,4254923.427;
0,0,0,0,192,1583.470126,6683.555129,19072.50505,40639.69537,66711.57632,81796.22718,60590.29663,-21861.2573,-171429.7577,-355062.5907,-492927.1661,-475398.2758,-208930.4549,320950.1747,997476.1882,1566265.895,1699768.472,1129529.762,-194208.6509;
0,0,0,0,192,1461.598622,5523.845413,13484.72812,22648.14037,23924.65349,4088.400074,-44523.67789,-110266.4937,-154796.0153,-124116.6074,18352.09343,250026.0416,467393.1673,513093.6957,258167.4211,-292568.7362,-936296.041,-1313188.428,-1070939.742;
0,0,0,0,192,1327.376095,4354.022806,8472.235861,8915.430969,-1775.330646,-26379.77787,-53287.95266,-55123.02039,-5131.095738,92963.75123,184430.9087,180646.3185,21137.78348,-252747.0452,-475735.34,-436097.8969,-33470.52706,583046.8515,1024123.527;
0,0,0,0,192,1184.107549,3229.595901,4321.651839,-104.2768198,-12455.61648,-25344.84529,-20852.91355,14642.78081,67162.133,88643.95935,29331.92183,-103721.9759,-215885.8219,-175668.2952,63339.97258,366040.1121,463007.7228,161145.8379,-426631.9431;
0,0,0,0,192,1035.320732,2197.653181,1181.837405,-4776.384006,-12614.95411,-10733.55999,10443.54925,39434.88701,40860.0242,-11919.57845,-90391.52336,-109544.5177,-4764.383231,169447.4417,238928.2291,63214.86651,-273067.0285,-452649.3406,-197417.4769;
0,0,0,0,192,884.6792681,1293.804398,-933.634351,-6033.445815,-7367.269353,3768.214436,22179.4179,22745.61092,-13409.71654,-59759.11695,-51183.66478,39553.74618,132447.0485,93359.81585,-98263.14966,-255824.9417,-144286.083,212696.4921,445204.5408;
0,0,0,0,192,735.8924507,540.3053094,-2124.343367,-5057.868211,-969.9812354,11405.43117,15940.52452,-5622.481002,-37623.47011,-29881.32575,35855.80045,87872.36995,28692.89466,-112944.5884,-154561.6584,24454.28757,254280.9736,202369.5839,-178659.2813;
0,0,0,0,192,592.6239049,-54.49033533,-2567.50857,-2988.296704,3941.876267,11526.38431,2233.953421,-22136.82236,-24139.16168,20459.96219,60268.45865,14532.0698,-88098.13253,-93277.95948,64910.42847,192595.8473,49803.13745,-244897.5404,-254777.9707;
0,0,0,0,192,458.4013779,-495.3380524,-2477.176593,-714.5839156,6392.684212,6873.986266,-9352.41825,-21134.06087,3844.754647,42020.47316,20505.10932,-58980.138,-70793.54333,51327.19794,141815.5672,6229.868714,-208390.0472,-130756.2533,224730.6889;
0,0,0,0,192,336.5298736,-798.086388,-2066.223865,1211.650343,6609.005295,832.2069828,-14542.95731,-9394.727768,23888.17206,29260.48864,-28610.53289,-63533.47519,18266.13986,110340.0964,20413.44929,-159630.1792,-99650.79156,191428.6434,224840.2275;
0,0,0,0,192,230.010273,-986.6727322,-1517.645248,2561.671612,5397.538212,-4248.815769,-13851.25499,3944.496017,28647.88724,2246.348802,-50351.46819,-19936.09886,77156.40402,55808.14847,-103790.0999,-116346.81,120768.0917,206007.3722,-114289.9057;
0,0,0,0,192,141.4654422,-1089.461027,-968.7547191,3346.025684,3627.717551,-7474.097775,-9965.085871,13492.9252,22500.21035,-20544.91549,-44251.15549,26525.3978,78409.54708,-27774.42335,-127872.1365,18879.12551,194650.9903,7368.517829,-279197.6607;
0,0,0,0,192,73.07564879,-1135.312342,-508.5655728,3706.949976,1946.428626,-9020.436512,-5497.089494,18315.25514,12841.39554,-32819.98602,-26305.72752,53592.96081,48907.83539,-81344.97478,-84369.59991,116249.2281,137093.641,-157744.2065,-212101.4947;
0,0,0,0,192,26.52487642,-1149.801347,-185.5323873,3822.418606,714.8963246,-9523.125445,-2036.182875,19914.17359,4805.440196,-36987.23947,-9962.723932,63039.66842,18780.08487,-100646.6018,-32907.74918,152629.087,54418.22672,-222018.2854,-85848.08848;
0,0,0,0,192,2.959359616,-1151.972632,-20.71532046,3839.781056,79.90093794,-9599.042127,-227.8620289,20156.93483,538.572741,-37623.95402,-1118.549352,64493.6094,2112.7619,-103642.0699,-3710.543262,158327.7533,6151.495274,-232190.8611,-9731.39151;
1,0,-1,0,1,0,-1,0,1,0,-1,0,1,0,-1,0,1,0,-1,0,1,0,-1,0;
0,0,4,0,-16,0,36,0,-64,0,100,0,-144,0,196,0,-256,0,324,0,-400,0,484,0;
0,0,0,24,192,840,2688,7056,16128,33264,63360,113256,192192,312312,489216,742560,1096704,1581408,2232576,3093048,4213440,5653032,7480704,9775920;
0,1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,441,484,529]
S=[1,0.998458667,0.993839419,0.986156496,0.975433581,0.96170373,0.945009268,0.925401657,0.902941342,0.87769756,0.84974813,0.819179209,0.786085032,0.750567618,0.712736454,0.672708162,0.630606134,0.586560159,0.540706014,0.493185053,0.444143767,0.393733335,0.342109153,0.289430364;
1,0.98618496,0.945121551,0.877944359,0.786509494,0.673343309,0.541572595,0.394838187,0.237194368,0.07299685,-0.093217576,-0.256856394,-0.413398249,-0.558517877,-0.688205612,-0.798878172,-0.887477663,-0.951556077,-0.98934292,-0.999794139,-0.982620967,-0.938297899,-0.868049586,-0.773816993;
1,0.961939766,0.850656228,0.67462034,0.447232036,0.18580022,-0.089774795,-0.35851611,-0.599967012,-0.795748145,-0.930956556,-0.99530012,-0.983880972,-0.897568346,-0.742932397,-0.531744087,-0.280079168,-0.007094493,0.266430219,0.519674138,0.733360219,0.891222577,0.981244656,0.996573933;
1,0.926320082,0.716137789,0.400425549,0.025706666,-0.352800347,-0.679318759,-0.90573287,-0.998678335,-0.944458724,-0.751063831,-0.446992295,-0.077052048,0.304242576,0.640704064,0.882751506,0.994716832,0.960100849,0.784004562,0.492377492,0.128193756,-0.254880591,-0.600395776,-0.857436738;
1,0.880202983,0.549514582,0.087165765,-0.396067449,-0.784405265,-0.984804259,-0.949250027,-0.686261152,-0.258848199,0.230583238,0.664768308,0.939678856,0.989447956,0.802151229,0.422663852,-0.058091262,-0.524928056,-0.86599522,-0.999575095,-0.89366274,-0.573634124,-0.116166194,0.369134463;
1,0.824724024,0.360339432,-0.230362851,-0.740310987,-0.990741662,-0.893865914,-0.483643725,0.096120716,0.642189852,0.963138082,0.946456378,0.597992543,0.039901255,-0.532177495,-0.917700386,-0.981521616,-0.701268527,-0.175184388,0.412310981,0.85526993,0.998412337,0.79155935,0.307223688;
1,0.761249282,0.15900094,-0.51917058,-0.949437402,-0.926346503,-0.460923818,0.224590651,0.802862762,0.997766752,0.716235686,0.092701052,-0.575098468,-0.968287643,-0.899118079,-0.400618342,0.289177228,0.840890257,0.991076981,0.668023024,0.025987114,-0.62845768,-0.98281303,-0.867873748;
1,0.691341716,-0.044093263,-0.75230874,-0.996111568,-0.624998222,0.131936882,0.807425162,0.984476513,0.553794202,-0.218754445,-0.856262349,-0.965185319,-0.4782834,0.303870785,0.8984405,0.93838801,0.399053054,-0.386623964,-0.933631603,-0.904292986,-0.316719326,0.46637042,0.96156198;
1,0.616722682,-0.239306267,-0.911893888,-0.885465021,-0.180278837,0.663100925,0.998177599,0.568096607,-0.297461473,-0.934999082,-0.855808809,-0.120594327,0.707062296,0.992717038,0.517399932,-0.354532491,-0.954696389,-0.823033344,-0.060470273,0.748446566,0.98363822,0.464817436,-0.410311308;
1,0.539229548,-0.418462989,-0.990524765,-0.649777453,0.289766361,0.96227862,0.74801177,-0.155578523,-0.915796843,-0.832070912,0.0184424,0.851960286,0.90036192,0.119043216,-0.771978681,-0.951590646,-0.254272907,0.677367717,0.984786282,0.384684007,-0.569920316,-0.999319756,-0.507805164;
1,0.460770452,-0.575381181,-0.991007746,-0.337872993,0.679643962,0.964192705,0.208899055,-0.771683681,-0.920037132,-0.07616817,0.849845048,0.859335144,-0.057932562,-0.91272237,-0.783178435,0.190991406,0.959184828,0.692936648,-0.320615363,-0.98839682,-0.590232736,0.444473211,0.99983298;
1,0.383277318,-0.706196995,-0.924615899,-0.002571609,0.92264462,0.70982912,-0.378521817,-0.999986774,-0.38802268,0.702546189,0.926562719,0.007714758,-0.920648935,-0.713442468,0.373756304,0.999947095,0.392757778,-0.698876799,-0.928485029,-0.012857704,0.918628896,0.717036943,-0.368980903;
1,0.308658284,-0.809460128,-0.808351431,0.310451397,0.999998222,0.306864073,-0.810565945,-0.807239861,0.312243405,0.999992888,0.305068772,-0.811668881,-0.80612542,0.314034304,0.999983998,0.303272386,-0.81276893,-0.805008112,0.315824085,0.999971552,0.301474921,-0.813866089,-0.803887941;
1,0.238750718,-0.88599619,-0.66181517,0.569978496,0.93398072,-0.124001362,-0.993191548,-0.350249028,0.825947135,0.74463997,-0.47038048,-0.969247324,0.007563492,0.972858903,0.456978031,-0.754651237,-0.81732508,0.364377338,0.991315782,0.10897737,-0.939278931,-0.557484408,0.673079326;
1,0.175275976,-0.938556665,-0.504288846,0.761777225,0.771331339,-0.491385519,-0.943587492,0.160609082,0.999889319,0.18990407,-0.933318077,-0.517080544,0.752054483,0.78071471,-0.478373418,-0.948409446,0.145906635,0.999557301,0.204490127,-0.927872888,-0.529757779,0.742165265,0.789925261;
1,0.119797017,-0.971297349,-0.352514068,0.886837082,0.564994942,-0.751467664,-0.745042111,0.572960019,0.882319914,-0.361561431,-0.968947876,0.1294073,0.999953093,0.110175495,-0.973555702,-0.343433634,0.891271052,0.556976861,-0.757822719,-0.738546663,0.580871344,0.877719972,-0.370574875;
1,0.073679918,-0.989142539,-0.2194398,0.956805927,0.360434564,-0.903692348,-0.49360252,0.830955162,0.616051936,-0.74017385,-0.725123833,0.633319721,0.818449723,-0.512713105,-0.894003042,0.380972963,0.950143155,-0.240960024,-0.985650985,0.095714657,0.999755481,0.051609146,-0.992150366;
1,0.038060234,-0.997102837,-0.113960168,0.988428136,0.18919978,-0.97402616,-0.263343106,0.95398036,0.335960537,-0.928406887,-0.406631304,0.897453922,0.474945916,-0.861300817,-0.540508536,0.820157054,0.602939275,-0.774261035,-0.661876387,0.723878695,0.716978371,-0.669301966,-0.76792595;
1,0.01381504,-0.999618289,-0.041434573,0.998473449,0.069022474,-0.996566352,-0.096557681,0.993898456,0.124019175,-0.990471796,-0.151385989,0.986288989,0.178637233,-0.981353228,-0.2057521,0.975668281,0.232709893,-0.969238489,-0.259490029,0.962068758,0.286072066,-0.954164565,-0.312435708;
1,0.001541333,-0.999995249,-0.004623985,0.999980994,0.007706592,-0.999957238,-0.010789127,0.999923978,0.013871559,-0.999881217,-0.016953859,0.999828954,0.020035998,-0.999767189,-0.023117946,0.999695925,0.026199675,-0.99961516,-0.029281155,0.999524896,0.032362357,-0.999425133,-0.035443251;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
CodePudding user response：
As far as I can tell the matrices you give have 2 complex evals, a complex conjugate pair. I checked using the Fortran program below, which uses both EISPACK and LAPACK, the latter having superseded the former decades ago. Both give the same answers where you can calculate an eval. I say this as note the comment at http://www.netlib.org/lapack/explore-3.1.1-html/dggev.f.html which is the documentation for the LAPACK routine:
* A generalized eigenvalue for a pair of matrices (A,B) is a scalar
* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
* singular. It is usually represented as the pair (alpha,beta), as
* there is a reasonable interpretation for beta=0, and even for both
* being zero.
In your case beta is zero in four cases, I am assuming due to the 4 columns of your matrices which contain zeros. In these cases you can't calculate lambda, so I quote the real and imaginary parts of alpha, and beta in the results below, before calculating lambda where possible.
LAPACK was the one acquired by apt get on my machine. The relevant parts of EISPACK were downloaded from http://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/eispack/rgg.f and the compiled into a little library with
ijb@ijb-Latitude-5410:~/Downloads/eispack$ unzip netlibfiles.zip
Archive: netlibfiles.zip
inflating: eispack/epslon.f
inflating: eispack/qzhes.f
inflating: eispack/qzit.f
inflating: eispack/qzval.f
inflating: eispack/qzvec.f
inflating: eispack/rgg.f
ijb@ijb-Latitude-5410:~/Downloads/eispack$ cd eispack/
ijb@ijb-Latitude-5410:~/Downloads/eispack/eispack$ gfortran -c -O *f
qzvec.f:92:72:
92 | 610 r = r (betm * a(i,j) - alfm * b(i,j)) * b(j,en)
| 1
Warning: Fortran 2018 deleted feature: DO termination statement which is not END DO or CONTINUE with label 610 at (1)
qzvec.f:207:72:
207 | do 880 i = 1, n
| 1
Warning: Fortran 2018 deleted feature: Shared DO termination label 880 at (1)
qzvec.f:211:72:
211 | 860 zz = zz z(i,k) * b(k,j)
| 1
Warning: Fortran 2018 deleted feature: DO termination statement which is not END DO or CONTINUE with label 860 at (1)
qzvec.f:228:72:
228 | 900 z(i,j) = z(i,j) / d
| 1
Warning: Fortran 2018 deleted feature: DO termination statement which is not END DO or CONTINUE with label 900 at (1)
ijb@ijb-Latitude-5410:~/Downloads/eispack/eispack$ ar r eispack.a *.o
ar: creating eispack.a
ijb@ijb-Latitude-5410:~/Downloads/eispack/eispack$
The data file was created from what you have above:
ijb@ijb-Latitude-5410:~/work/stack$ cat eig.dat
0 0 0 0 192 1917.04064 10332.51505 40092.51227 125681.1486 338350.2206 811892.8294 1779728.921 3625982.355 6953387.916 12670976.81 22100000 37132930.27 60353006.25 95276316.19 146559937.4 220274060.5 324208411.9 468219308 664618721.8
0 0 0 0 192 1893.475124 10051.90014 38308.22391 117609.6433 309187.9535 722364.2569 1537115.973 3030677.025 5606681.841 9824567.083 16426180.74 26355891.41 40770017.02 61031174.08 88683247.04 125403139.6 172926323.5 232944447.3 306974887
0 0 0 0 192 1846.924351 9507.779872 34923.83828 102685.9158 256815.0258 566753.8301 1130503.089 2072244.588 3531888.862 5644710.273 8510616.36 12154554.65 16481872.51 21234784.98 25958063.51 29983201.26 32440241.94 32304892.86 28485362.66
0 0 0 0 192 1778.534558 8732.953666 30281.20319 83090.50097 191433.6927 383369.9941 681577.7806 1089000.57 1571702.814 2043718.6 2360415.659 2327150.205 1728131.214 376286.8361 -1820993.208 -4791939.112 -8243585.409 -11638086.16 -14220529.21
0 0 0 0 192 1689.989727 7773.203992 24831.36406 61515.88289 124689.9821 212051.9399 303583.89 356337.2627 308129.1939 94242.56924 -323226.8121 -920239.2345 -1583058.32 -2104105.201 -2214813.175 -1659107.74 -295993.1841 1796496.926 4254923.427
0 0 0 0 192 1583.470126 6683.555129 19072.50505 40639.69537 66711.57632 81796.22718 60590.29663 -21861.2573 -171429.7577 -355062.5907 -492927.1661 -475398.2758 -208930.4549 320950.1747 997476.1882 1566265.895 1699768.472 1129529.762 -194208.6509
0 0 0 0 192 1461.598622 5523.845413 13484.72812 22648.14037 23924.65349 4088.400074 -44523.67789 -110266.4937 -154796.0153 -124116.6074 18352.09343 250026.0416 467393.1673 513093.6957 258167.4211 -292568.7362 -936296.041 -1313188.428 -1070939.742
0 0 0 0 192 1327.376095 4354.022806 8472.235861 8915.430969 -1775.330646 -26379.77787 -53287.95266 -55123.02039 -5131.095738 92963.75123 184430.9087 180646.3185 21137.78348 -252747.0452 -475735.34 -436097.8969 -33470.52706 583046.8515 1024123.527
0 0 0 0 192 1184.107549 3229.595901 4321.651839 -104.2768198 -12455.61648 -25344.84529 -20852.91355 14642.78081 67162.133 88643.95935 29331.92183 -103721.9759 -215885.8219 -175668.2952 63339.97258 366040.1121 463007.7228 161145.8379 -426631.9431
0 0 0 0 192 1035.320732 2197.653181 1181.837405 -4776.384006 -12614.95411 -10733.55999 10443.54925 39434.88701 40860.0242 -11919.57845 -90391.52336 -109544.5177 -4764.383231 169447.4417 238928.2291 63214.86651 -273067.0285 -452649.3406 -197417.4769
0 0 0 0 192 884.6792681 1293.804398 -933.634351 -6033.445815 -7367.269353 3768.214436 22179.4179 22745.61092 -13409.71654 -59759.11695 -51183.66478 39553.74618 132447.0485 93359.81585 -98263.14966 -255824.9417 -144286.083 212696.4921 445204.5408
0 0 0 0 192 735.8924507 540.3053094 -2124.343367 -5057.868211 -969.9812354 11405.43117 15940.52452 -5622.481002 -37623.47011 -29881.32575 35855.80045 87872.36995 28692.89466 -112944.5884 -154561.6584 24454.28757 254280.9736 202369.5839 -178659.2813
0 0 0 0 192 592.6239049 -54.49033533 -2567.50857 -2988.296704 3941.876267 11526.38431 2233.953421 -22136.82236 -24139.16168 20459.96219 60268.45865 14532.0698 -88098.13253 -93277.95948 64910.42847 192595.8473 49803.13745 -244897.5404 -254777.9707
0 0 0 0 192 458.4013779 -495.3380524 -2477.176593 -714.5839156 6392.684212 6873.986266 -9352.41825 -21134.06087 3844.754647 42020.47316 20505.10932 -58980.138 -70793.54333 51327.19794 141815.5672 6229.868714 -208390.0472 -130756.2533 224730.6889
0 0 0 0 192 336.5298736 -798.086388 -2066.223865 1211.650343 6609.005295 832.2069828 -14542.95731 -9394.727768 23888.17206 29260.48864 -28610.53289 -63533.47519 18266.13986 110340.0964 20413.44929 -159630.1792 -99650.79156 191428.6434 224840.2275
0 0 0 0 192 230.010273 -986.6727322 -1517.645248 2561.671612 5397.538212 -4248.815769 -13851.25499 3944.496017 28647.88724 2246.348802 -50351.46819 -19936.09886 77156.40402 55808.14847 -103790.0999 -116346.81 120768.0917 206007.3722 -114289.9057
0 0 0 0 192 141.4654422 -1089.461027 -968.7547191 3346.025684 3627.717551 -7474.097775 -9965.085871 13492.9252 22500.21035 -20544.91549 -44251.15549 26525.3978 78409.54708 -27774.42335 -127872.1365 18879.12551 194650.9903 7368.517829 -279197.6607
0 0 0 0 192 73.07564879 -1135.312342 -508.5655728 3706.949976 1946.428626 -9020.436512 -5497.089494 18315.25514 12841.39554 -32819.98602 -26305.72752 53592.96081 48907.83539 -81344.97478 -84369.59991 116249.2281 137093.641 -157744.2065 -212101.4947
0 0 0 0 192 26.52487642 -1149.801347 -185.5323873 3822.418606 714.8963246 -9523.125445 -2036.182875 19914.17359 4805.440196 -36987.23947 -9962.723932 63039.66842 18780.08487 -100646.6018 -32907.74918 152629.087 54418.22672 -222018.2854 -85848.08848
0 0 0 0 192 2.959359616 -1151.972632 -20.71532046 3839.781056 79.90093794 -9599.042127 -227.8620289 20156.93483 538.572741 -37623.95402 -1118.549352 64493.6094 2112.7619 -103642.0699 -3710.543262 158327.7533 6151.495274 -232190.8611 -9731.39151
1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
0 0 4 0 -16 0 36 0 -64 0 100 0 -144 0 196 0 -256 0 324 0 -400 0 484 0
0 0 0 24 192 840 2688 7056 16128 33264 63360 113256 192192 312312 489216 742560 1096704 1581408 2232576 3093048 4213440 5653032 7480704 9775920
0 1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256 289 324 361 400 441 484 529
1 0.998458667 0.993839419 0.986156496 0.975433581 0.96170373 0.945009268 0.925401657 0.902941342 0.87769756 0.84974813 0.819179209 0.786085032 0.750567618 0.712736454 0.672708162 0.630606134 0.586560159 0.540706014 0.493185053 0.444143767 0.393733335 0.342109153 0.289430364
1 0.98618496 0.945121551 0.877944359 0.786509494 0.673343309 0.541572595 0.394838187 0.237194368 0.07299685 -0.093217576 -0.256856394 -0.413398249 -0.558517877 -0.688205612 -0.798878172 -0.887477663 -0.951556077 -0.98934292 -0.999794139 -0.982620967 -0.938297899 -0.868049586 -0.773816993
1 0.961939766 0.850656228 0.67462034 0.447232036 0.18580022 -0.089774795 -0.35851611 -0.599967012 -0.795748145 -0.930956556 -0.99530012 -0.983880972 -0.897568346 -0.742932397 -0.531744087 -0.280079168 -0.007094493 0.266430219 0.519674138 0.733360219 0.891222577 0.981244656 0.996573933
1 0.926320082 0.716137789 0.400425549 0.025706666 -0.352800347 -0.679318759 -0.90573287 -0.998678335 -0.944458724 -0.751063831 -0.446992295 -0.077052048 0.304242576 0.640704064 0.882751506 0.994716832 0.960100849 0.784004562 0.492377492 0.128193756 -0.254880591 -0.600395776 -0.857436738
1 0.880202983 0.549514582 0.087165765 -0.396067449 -0.784405265 -0.984804259 -0.949250027 -0.686261152 -0.258848199 0.230583238 0.664768308 0.939678856 0.989447956 0.802151229 0.422663852 -0.058091262 -0.524928056 -0.86599522 -0.999575095 -0.89366274 -0.573634124 -0.116166194 0.369134463
1 0.824724024 0.360339432 -0.230362851 -0.740310987 -0.990741662 -0.893865914 -0.483643725 0.096120716 0.642189852 0.963138082 0.946456378 0.597992543 0.039901255 -0.532177495 -0.917700386 -0.981521616 -0.701268527 -0.175184388 0.412310981 0.85526993 0.998412337 0.79155935 0.307223688
1 0.761249282 0.15900094 -0.51917058 -0.949437402 -0.926346503 -0.460923818 0.224590651 0.802862762 0.997766752 0.716235686 0.092701052 -0.575098468 -0.968287643 -0.899118079 -0.400618342 0.289177228 0.840890257 0.991076981 0.668023024 0.025987114 -0.62845768 -0.98281303 -0.867873748
1 0.691341716 -0.044093263 -0.75230874 -0.996111568 -0.624998222 0.131936882 0.807425162 0.984476513 0.553794202 -0.218754445 -0.856262349 -0.965185319 -0.4782834 0.303870785 0.8984405 0.93838801 0.399053054 -0.386623964 -0.933631603 -0.904292986 -0.316719326 0.46637042 0.96156198
1 0.616722682 -0.239306267 -0.911893888 -0.885465021 -0.180278837 0.663100925 0.998177599 0.568096607 -0.297461473 -0.934999082 -0.855808809 -0.120594327 0.707062296 0.992717038 0.517399932 -0.354532491 -0.954696389 -0.823033344 -0.060470273 0.748446566 0.98363822 0.464817436 -0.410311308
1 0.539229548 -0.418462989 -0.990524765 -0.649777453 0.289766361 0.96227862 0.74801177 -0.155578523 -0.915796843 -0.832070912 0.0184424 0.851960286 0.90036192 0.119043216 -0.771978681 -0.951590646 -0.254272907 0.677367717 0.984786282 0.384684007 -0.569920316 -0.999319756 -0.507805164
1 0.460770452 -0.575381181 -0.991007746 -0.337872993 0.679643962 0.964192705 0.208899055 -0.771683681 -0.920037132 -0.07616817 0.849845048 0.859335144 -0.057932562 -0.91272237 -0.783178435 0.190991406 0.959184828 0.692936648 -0.320615363 -0.98839682 -0.590232736 0.444473211 0.99983298
1 0.383277318 -0.706196995 -0.924615899 -0.002571609 0.92264462 0.70982912 -0.378521817 -0.999986774 -0.38802268 0.702546189 0.926562719 0.007714758 -0.920648935 -0.713442468 0.373756304 0.999947095 0.392757778 -0.698876799 -0.928485029 -0.012857704 0.918628896 0.717036943 -0.368980903
1 0.308658284 -0.809460128 -0.808351431 0.310451397 0.999998222 0.306864073 -0.810565945 -0.807239861 0.312243405 0.999992888 0.305068772 -0.811668881 -0.80612542 0.314034304 0.999983998 0.303272386 -0.81276893 -0.805008112 0.315824085 0.999971552 0.301474921 -0.813866089 -0.803887941
1 0.238750718 -0.88599619 -0.66181517 0.569978496 0.93398072 -0.124001362 -0.993191548 -0.350249028 0.825947135 0.74463997 -0.47038048 -0.969247324 0.007563492 0.972858903 0.456978031 -0.754651237 -0.81732508 0.364377338 0.991315782 0.10897737 -0.939278931 -0.557484408 0.673079326
1 0.175275976 -0.938556665 -0.504288846 0.761777225 0.771331339 -0.491385519 -0.943587492 0.160609082 0.999889319 0.18990407 -0.933318077 -0.517080544 0.752054483 0.78071471 -0.478373418 -0.948409446 0.145906635 0.999557301 0.204490127 -0.927872888 -0.529757779 0.742165265 0.789925261
1 0.119797017 -0.971297349 -0.352514068 0.886837082 0.564994942 -0.751467664 -0.745042111 0.572960019 0.882319914 -0.361561431 -0.968947876 0.1294073 0.999953093 0.110175495 -0.973555702 -0.343433634 0.891271052 0.556976861 -0.757822719 -0.738546663 0.580871344 0.877719972 -0.370574875
1 0.073679918 -0.989142539 -0.2194398 0.956805927 0.360434564 -0.903692348 -0.49360252 0.830955162 0.616051936 -0.74017385 -0.725123833 0.633319721 0.818449723 -0.512713105 -0.894003042 0.380972963 0.950143155 -0.240960024 -0.985650985 0.095714657 0.999755481 0.051609146 -0.992150366
1 0.038060234 -0.997102837 -0.113960168 0.988428136 0.18919978 -0.97402616 -0.263343106 0.95398036 0.335960537 -0.928406887 -0.406631304 0.897453922 0.474945916 -0.861300817 -0.540508536 0.820157054 0.602939275 -0.774261035 -0.661876387 0.723878695 0.716978371 -0.669301966 -0.76792595
1 0.01381504 -0.999618289 -0.041434573 0.998473449 0.069022474 -0.996566352 -0.096557681 0.993898456 0.124019175 -0.990471796 -0.151385989 0.986288989 0.178637233 -0.981353228 -0.2057521 0.975668281 0.232709893 -0.969238489 -0.259490029 0.962068758 0.286072066 -0.954164565 -0.312435708
1 0.001541333 -0.999995249 -0.004623985 0.999980994 0.007706592 -0.999957238 -0.010789127 0.999923978 0.013871559 -0.999881217 -0.016953859 0.999828954 0.020035998 -0.999767189 -0.023117946 0.999695925 0.026199675 -0.99961516 -0.029281155 0.999524896 0.032362357 -0.999425133 -0.035443251
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
The program, its compilation are
ijb@ijb-Latitude-5410:~/work/stack$ cat eig.f90
Program eig
Use, Intrinsic :: iso_fortran_env, Only : wp => real64
Implicit None
Integer, Parameter :: n = 24
Real( wp ), Dimension( 1:n, 1:n ) :: H
Real( wp ), Dimension( 1:n, 1:n ) :: S
Integer :: unit
Open( newunit = unit, file = 'eig.dat' )
Read( unit, * ) H
Read( unit, * ) S
Close( unit )
Write( *, * ) 'Lapack: '
Call lapack( H, S )
Write( *, * ) 'Eispack: '
Call eispack( H, S )
Contains
Subroutine lapack( H, S )
Use, Intrinsic :: iso_fortran_env, Only : wp => real64
Implicit None
Real( wp ), Dimension( :, : ), Intent( In ) :: H
Real( wp ), Dimension( :, : ), Intent( In ) :: S
Complex( wp ) :: lambda
Real( wp ), Dimension( :, : ), Allocatable :: A
Real( wp ), Dimension( :, : ), Allocatable :: B
Real( wp ), Dimension( :, : ), Allocatable :: QR
Real( wp ), Dimension( 1:1, 1:1 ) :: Ql_dummy
Real( wp ), Dimension( : ), Allocatable :: alphar
Real( wp ), Dimension( : ), Allocatable :: alphai
Real( wp ), Dimension( : ), Allocatable :: beta
Real( wp ), Dimension( : ), Allocatable :: work
Real( wp ), Dimension( 1:1 ) :: twork
Integer :: n
Integer :: lwork
Integer :: info
Integer :: i
n = Size( H, Dim = 1 )
! H and S would be overwritten
A = H
B = S
! Space for evals and evecs
Allocate( alphar( 1:n ) )
Allocate( alphai( 1:n ) )
Allocate( beta( 1:n ) )
Allocate( QR( 1:n, 1:n ) )
! Calculate and allocate the worksize
Call dggev( 'N', 'V', n, A, Size( A, Dim = 1 ), B, Size( B, Dim = 1 ), &
alphar, alphai, beta, &
QL_dummy, Size( QL_dummy, Dim = 1 ), &
QR, Size( QR, Dim = 1 ), &
twork, -1, info )
! Allocate workspace
lwork = Nint( twork( 1 ) )
Allocate( work( 1:lwork ) )
! Solve eval problem
Call dggev( 'N', 'V', n, A, Size( A, Dim = 1 ), B, Size( B, Dim = 1 ), &
alphar, alphai, beta, &
QL_dummy, Size( QL_dummy, Dim = 1 ), &
QR, Size( QR, Dim = 1 ), &
work, Size( work ), info )
Write( *, * ) 'info = ', info
! Report results
Write( *, * ) 'Evals - last column reports lambda where appropriate'
Do i = 1, n
If( Abs( beta( i ) ) > 1.0e-16_wp ) Then
lambda = Cmplx( alphar( i ), alphai( i ), wp ) / beta( i )
Write( *, '( i2, 1x, 3( e18.6, 1x ), 5x, e18.6, " ", e18.6, "i" )' ) &
i, alphar( i ), alphai( i ), beta( i ), lambda
Else
Write( *, '( i2, 1x, 3( e18.6, 1x ) )' ) i, alphar( i ), alphai( i ), beta( i )
End If
End Do
Write( *, * )
End Subroutine lapack
Subroutine eispack( H, S )
Use, Intrinsic :: iso_fortran_env, Only : wp => real64
Implicit None
Real( wp ), Dimension( :, : ), Intent( In ) :: H
Real( wp ), Dimension( :, : ), Intent( In ) :: S
Complex( wp ) :: lambda
Real( wp ), Dimension( :, : ), Allocatable :: A
Real( wp ), Dimension( :, : ), Allocatable :: B
Real( wp ), Dimension( :, : ), Allocatable :: QR
Real( wp ), Dimension( : ), Allocatable :: alphar
Real( wp ), Dimension( : ), Allocatable :: alphai
Real( wp ), Dimension( : ), Allocatable :: beta
Integer :: n
Integer :: info
Integer :: i
n = Size( H, Dim = 1 )
! H and S would be overwritten
A = H
B = S
! Space for evals and evecs
Allocate( alphar( 1:n ) )
Allocate( alphai( 1:n ) )
Allocate( beta( 1:n ) )
Allocate( QR( 1:n, 1:n ) )
! Calc evals and evecs by eispack
Call rgg( Size( A, Dim = 1 ), n, A, B, alphar, alphai, beta, 1, QR, info )
Write( *, * ) 'info = ', info
! Report results
Write( *, * ) 'Evals - last column reports lambda where appropriate'
Do i = 1, n
If( Abs( beta( i ) ) > 1.0e-16_wp ) Then
lambda = Cmplx( alphar( i ), alphai( i ), wp ) / beta( i )
Write( *, '( i2, 1x, 3( e18.6, 1x ), 5x, e18.6, " ", e18.6, "i" )' ) &
i, alphar( i ), alphai( i ), beta( i ), lambda
Else
Write( *, '( i2, 1x, 3( e18.6, 1x ) )' ) i, alphar( i ), alphai( i ), beta( i )
End If
End Do
Write( *, * )
End Subroutine eispack
End Program eig
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -O -g eig.f90 eispack.a -llapack
And the results are (note the complex conjugate pair)
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
Lapack:
info = 0
Evals - last column reports lambda where appropriate
1 -0.608116E 01 0.000000E 00 0.187802E-11 -0.323806E 13 0.000000E 00i
2 0.370614E 03 0.000000E 00 0.262593E-05 0.141137E 09 0.000000E 00i
3 -0.166285E 05 0.000000E 00 0.150764E-02 -0.110295E 08 0.000000E 00i
4 0.356929E 06 0.000000E 00 0.387651E-01 0.920748E 07 0.000000E 00i
5 0.376034E 06 0.000000E 00 0.654174E-01 0.574823E 07 0.000000E 00i
6 0.881731E 05 0.129854E 06 0.121455E 00 0.725976E 06 0.106916E 07i
7 0.154526E 06 -0.227573E 06 0.212853E 00 0.725976E 06 -0.106916E 07i
8 0.400747E 05 0.000000E 00 0.327119E-01 0.122508E 07 0.000000E 00i
9 0.353978E 05 0.000000E 00 0.623272E-01 0.567936E 06 0.000000E 00i
10 0.177191E 05 0.000000E 00 0.342541E-01 0.517282E 06 0.000000E 00i
11 0.703041E 04 0.000000E 00 0.224366E-01 0.313346E 06 0.000000E 00i
12 0.250241E 04 0.000000E 00 0.142464E-01 0.175653E 06 0.000000E 00i
13 0.816411E 02 0.000000E 00 0.815119E-03 0.100159E 06 0.000000E 00i
14 0.127044E 04 0.000000E 00 0.141856E-01 0.895583E 05 0.000000E 00i
15 0.396171E 04 0.000000E 00 0.991827E-01 0.399436E 05 0.000000E 00i
16 0.461457E 05 0.000000E 00 0.315689E 01 0.146175E 05 0.000000E 00i
17 0.109954E 05 0.000000E 00 0.288969E 01 0.380504E 04 0.000000E 00i
18 0.188721E 02 0.000000E 00 0.309986E 01 0.608807E 01 0.000000E 00i
19 0.123631E-03 0.000000E 00 0.139487E-06 0.886324E 03 0.000000E 00i
20 0.158220E 04 0.000000E 00 0.320845E 01 0.493134E 03 0.000000E 00i
21 0.899768E 03 0.000000E 00 0.000000E 00
22 0.376730E-01 0.000000E 00 0.000000E 00
23 0.366102E-01 0.000000E 00 0.000000E 00
24 0.960100E-01 0.000000E 00 0.000000E 00
Eispack:
info = 0
Evals - last column reports lambda where appropriate
1 -0.346410E 01 0.000000E 00 0.000000E 00
2 -0.546079E 03 0.000000E 00 0.000000E 00
3 0.144715E 08 0.000000E 00 0.000000E 00
4 0.528475E 03 0.000000E 00 0.000000E 00
5 -0.233202E 08 0.000000E 00 0.719984E-05 -0.323898E 13 0.000000E 00i
6 0.278909E 06 0.000000E 00 0.197617E-02 0.141136E 09 0.000000E 00i
7 -0.250582E 05 0.000000E 00 0.227192E-02 -0.110295E 08 0.000000E 00i
8 0.491054E 04 0.000000E 00 0.533321E-03 0.920746E 07 0.000000E 00i
9 0.777637E 05 0.000000E 00 0.135281E-01 0.574829E 07 0.000000E 00i
10 0.199724E 04 0.294137E 04 0.275110E-02 0.725977E 06 0.106916E 07i
11 0.335291E 04 -0.493790E 04 0.461849E-02 0.725977E 06 -0.106916E 07i
12 0.358837E 04 0.000000E 00 0.292908E-02 0.122508E 07 0.000000E 00i
13 0.981532E 03 0.000000E 00 0.172823E-02 0.567941E 06 0.000000E 00i
14 0.138597E 03 0.000000E 00 0.267934E-03 0.517279E 06 0.000000E 00i
15 0.286093E 03 0.000000E 00 0.913029E-03 0.313345E 06 0.000000E 00i
16 0.518162E 02 0.000000E 00 0.294991E-03 0.175653E 06 0.000000E 00i
17 0.957124E 00 0.000000E 00 0.955594E-05 0.100160E 06 0.000000E 00i
18 0.104269E 01 0.000000E 00 0.116426E-04 0.895582E 05 0.000000E 00i
19 0.874474E 02 0.000000E 00 0.218927E-02 0.399436E 05 0.000000E 00i
20 0.185986E 04 0.000000E 00 0.127235E 00 0.146175E 05 0.000000E 00i
21 0.968367E 04 0.000000E 00 0.254496E 01 0.380504E 04 0.000000E 00i
22 0.171890E 02 0.000000E 00 0.282339E 01 0.608807E 01 0.000000E 00i
23 0.100843E-04 0.000000E 00 0.113655E-07 0.887279E 03 0.000000E 00i
24 0.234040E 03 0.000000E 00 0.474598E 00 0.493134E 03 0.000000E 00i
Given all this I suggest one of the following is occurring:
- The original paper is wrong
- The data you quote is wrong
- I am misunderstanding something you are saying
I don't know if it is possible a priori to say whether a given problem has a complex eval without solving the problem itself.