我正在嘗試解決廣義特征值問題。我有兩個矩陣 H 和 S 使得:
HX=λSX
我需要找到特征值λ。矩陣 H 和 S 是實數、不對稱的,包含正負十進制數和(可以)單數。我使用命令eig(H,S)在 MATLAB 上進行了嘗試,但問題是我獲得了實數和復數特征值,而在我關注的作者的研究論文中,只有實數特征值。
在閱讀了這個問題之后,我了解到對于這樣的矩陣,像 MATLAB 和許多其他軟體使用 QZ 演算法來解決廣義特征值問題。我正在尋找以下問題的答案:
- 是否有任何標準來確定給定矩陣是否具有復雜的特征值?并且從 MATLAB 獲得的結果是否正確?(盡管這些值與研究論文不一致)。
- 在研究論文中,作者使用了 RGG 子程式,它是 Fortran 庫 EISPACK 的一部分。我瀏覽了它的檔案并知道它也使用了 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]
uj5u.com熱心網友回復:
據我所知,您給出的矩陣有 2 個復數 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.
在您的情況下,beta 在四種情況下為零,我假設由于您的矩陣的 4 列包含零。在這些情況下,您無法計算 lambda,因此在可能的情況下計算 lambda 之前,我在下面的結果中參考了 alpha 和 beta 的實部和虛部。
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
鑒于這一切,我建議正在發生以下情況之一:
- 原紙有誤
- 你參考的資料是錯誤的
- 我誤解了你在說什么
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.
轉載請註明出處,本文鏈接:https://www.uj5u.com/shujuku/437394.html
