使用MATLAB的广义特征值问题



我正试图解决一个广义特征值问题。我有两个矩阵H和S,这样:

HX=λSX

我需要找到特征值λ。矩阵H和S是实的、非对称的,包含正小数和负小数,并且(可以(是奇异的。我在MATLAB上使用命令eig(H,S(尝试了这一点,但问题是我正在获得真实和复杂的特征值,而在我跟随作者的研究论文中,我只获得了真实的特征值。

在阅读了这个问题之后,我知道对于这样的矩阵,像MATLAB和许多其他软件都使用QZ算法来解决广义本征值问题。我正在寻求以下问题的答案:

  1. 是否有任何标准来确定给定矩阵是否具有复数特征值?并且从MATLAB获得的结果是正确的还是不正确的?(尽管价值观与研究论文不一致(
  2. 本文使用了Fortran库EISPACK中的RGG子程序。我查阅了它的文档,了解到它也使用了QZ算法。所以我的问题是,即使MATLAB和Fortran都使用相同的算法,它们能对同一个问题给出不同的答案吗

以下是矩阵:

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]

据我所知,你给出的矩阵有两个复数eval,一个复数共轭对。我使用下面的Fortran程序进行了检查,该程序同时使用EISPACK和LAPACK,后者在几十年前已经取代了前者。在计算eval的地方,两者给出的答案相同。我这样说是为了注意http://www.netlib.org/lapack/explore-3.1.1-html/dggev.f.html这是LAPACK程序的文档:

*  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.

在你的情况下,贝塔在四种情况下都是零,我假设这是因为你的矩阵中有四列包含零。在这些情况下,你无法计算lambda,所以我在下面的结果中引用了alpha和beta的实部和虚部,然后在可能的情况下计算lambda。

LAPACK是apt-get在我的机器上获得的。EISPACK的相关部分可从http://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/eispack/rgg.f并用编译成一个小库

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$ 

数据文件是根据上面的内容创建的:

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

该程序,其编译是

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

结果是(注意复共轭对(

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

考虑到所有这些,我认为以下情况之一正在发生:

  • 原始论文错误
  • 你引用的数据是错误的
  • 我误解了你说的话

我不知道是否有可能先验地在不解决问题本身的情况下判断给定问题是否具有复杂的eval。

相关内容

  • 没有找到相关文章

最新更新