使用Picard迭代中的矩陣列表優化計算

Question

當前我正在使用一些Mathematica代碼來執行Picard迭代。代碼本身工作正常，但我試圖使它更有效率。我取得了一些成功，但正在尋找建議。它可能無法加速了，但我已經用完了想法，希望擁有比我更多編程/ Mathematica經驗的人可能會提出一些建議。我只發佈迭代本身，但可以根據需要提供其他信息。下面

的代碼被編輯成完全可執行的要求

此外，我改變了它從當一個Do循環，使測試，並不需要收斂容易。

Clear["Global`*"] 

ngrid = 2048; 
delr = 4/100; 
delk = \[Pi]/delr/ngrid; 
rvalues = Table[(i - 1/2) delr, {i, 1, ngrid}]; 
kvalues = Table[(i - 1/2) delk, {i, 1, ngrid}]; 
wa[x_] := (19 + .5 x) Exp[-.7 x] + 1 
wb[x_] := (19 + .1 x) Exp[-.2 x] + 1 
wd = SetPrecision[ 
    Table[{{wa[(i - 1/2) delk], 0}, {0, wb[(i - 1/2) delk]}}, {i, 1, 
    ngrid}], 26]; 
sigmaAA = 1; 
hcloseAA = {}; 
i = 1; 
While[(i - 1/2)*delr < sigmaAA, hcloseAA = Append[hcloseAA, -1]; i++] 
hcloselenAA = Length[hcloseAA]; 
hcloseAB = hcloseAA; 
hcloselenAB = hcloselenAA; 
hcloseBB = hcloseAA; 
hcloselenBB = hcloselenAA; 
ccloseAA = {}; 
i = ngrid; 
While[(i - 1/2)*delr >= sigmaAA, ccloseAA = Append[ccloseAA, 0]; i--] 
ccloselenAA = Length[ccloseAA]; 
ccloselenAA = Length[ccloseAA]; 
ccloseAB = ccloseAA; 
ccloselenAB = ccloselenAA; 
ccloseBB = ccloseAA; 
ccloselenBB = ccloselenAA; 
na = 20; 
nb = 20; 
pa = 27/(1000 \[Pi]); 
pb = 27/(1000 \[Pi]); 
p = {{na pa, 0}, {0, nb pb}}; 
id = {{1, 0}, {0, 1}}; 
AFD = 1; 
AFDList = {}; 
timelist = {}; 
gammainitial = Table[{{0, 0}, {0, 0}}, {ngrid}]; 
gammafirst = gammainitial; 
step = 1; 
tol = 10^-7; 
old = 95/100; 
new = 1 - old; 

Do[ 
t = AbsoluteTime[]; 
extractgAA = Table[Extract[gammafirst, {i, 1, 1}], {i, hcloselenAA}]; 
extractgBB = Table[Extract[gammafirst, {i, 2, 2}], {i, hcloselenBB}]; 
extractgAB = Table[Extract[gammafirst, {i, 1, 2}], {i, hcloselenAB}]; 
csolutionAA = (Join[hcloseAA - extractgAA, ccloseAA]) rvalues; 
csolutionBB = (Join[hcloseBB - extractgBB, ccloseBB]) rvalues; 
csolutionAB = (Join[hcloseAB - extractgAB, ccloseAB]) rvalues; 
chatAA = FourierDST[SetPrecision[csolutionAA, 32], 4]; 
chatBB = FourierDST[SetPrecision[csolutionBB, 32], 4]; 
chatAB = FourierDST[SetPrecision[csolutionAB, 32], 4]; 
chatmatrix = 
    2 \[Pi] delr Sqrt[2*ngrid]* 
    Transpose[{Transpose[{chatAA, chatAB}], 
     Transpose[{chatAB, chatBB}]}]/kvalues; 
gammahat = 
    Table[(wd[[i]].chatmatrix[[i]].(Inverse[ 
     id - p.wd[[i]].chatmatrix[[i]]]).wd[[i]] - 
     chatmatrix[[i]]) kvalues[[i]], {i, ngrid}]; 
gammaAA = 
    FourierDST[SetPrecision[Table[gammahat[[i, 1, 1]], {i, ngrid}], 32], 
    4]; 
gammaBB = 
    FourierDST[SetPrecision[Table[gammahat[[i, 2, 2]], {i, ngrid}], 32], 
    4]; 
gammaAB = 
    FourierDST[SetPrecision[Table[gammahat[[i, 1, 2]], {i, ngrid}], 32], 
    4]; 
gammasecond = 
    Transpose[{Transpose[{gammaAA, gammaAB}], 
    Transpose[{gammaAB, gammaBB}]}]/(rvalues 2 \[Pi] delr Sqrt[ 
     2*ngrid]); 
AFD = Sqrt[ 
    1/ngrid Sum[((gammafirst[[i, 1, 1]] - 
      gammasecond[[i, 1, 1]])/(gammafirst[[i, 1, 1]] + 
      gammasecond[[i, 1, 1]]))^2 + ((gammafirst[[i, 2, 2]] - 
      gammasecond[[i, 2, 2]])/(gammafirst[[i, 2, 2]] + 
      gammasecond[[i, 2, 2]]))^2 + ((gammafirst[[i, 1, 2]] - 
      gammasecond[[i, 1, 2]])/(gammafirst[[i, 1, 2]] + 
      gammasecond[[i, 1, 2]]))^2 + ((gammafirst[[i, 2, 1]] - 
      gammasecond[[i, 2, 1]])/(gammafirst[[i, 2, 1]] + 
      gammasecond[[i, 2, 1]]))^2, {i, 1, ngrid}]]; 
gammafirst = old gammafirst + new gammasecond; 
time2 = AbsoluteTime[] - t; 
timelist = Append[timelist, time2], {1}] 
Print["Mean time per calculation = ", Mean[timelist]] 
Print["STD time per calculation = ", StandardDeviation[timelist]] 

只是對事物的一些注意事項
ngrid，delr，德爾克，右值，kvalues都只是使問題離散使用的值。典型地，它們是

ngrid = 2048; 
delr = 4/100; 
delk = \[Pi]/delr/ngrid; 
rvalues = Table[(i - 1/2) delr, {i, 1, ngrid}]; 
kvalues = Table[(i - 1/2) delk, {i, 1, ngrid}]; 

所有正在使用的矩陣是2×2具有相同斷開對角線

單位矩陣和所述P矩陣（它實際上是用於密度）是

p = {{na pa, 0}, {0, nb pb}}; 
id = {{1, 0}, {0, 1}}; 

我已確定的計算中的主要慢點是FourierDST計算（前向和後向變換佔計算時間的近40％）。伽瑪計算佔40％的時間，其餘時間由AFD ca控制） 在我的i7處理器上，每個週期的平均計算時間爲1.52秒。我的希望是不到一秒鐘，但這可能是不可能的。 我的希望是引入一些並行計算，這是與ParallelTable命令以及使用ParallelSubmitWaitAll兩個嘗試。然而，我發現任何來自並行計算的加速都被從主內核到其他內核的通信時間所抵消（至少這是我的假設，因爲對新數據的計算需要重新計算現有數據的兩倍。我認爲這意味着減速是在傳播新名單）我玩DistributDefinitions以及SetSharedVariable，然而，無法做到這一點。

我想知道的一件事是，如果使用Table做離散計算是做到這一點的最好方法？

我也曾想過我可以重寫這樣的方式，以便能夠編譯它，但我的理解是，只有當你正在處理機器精度，我需要以更高的精度工作以獲得收斂。

非常感謝您的任何建議。

Answer 1

我會等待代碼ACL建議，但關上，我懷疑這個結構：

Table[Extract[gammafirst, {i, 1, 1}], {i, hcloselenAA}] 

可以寫入，並且將執行速度更快，因爲：

gammafirst[[hcloselenAA, 1, 1]] 

但我被迫猜測你的數據的形狀。

Answer 2

在幾行使用：

FourierDST[SetPrecision[Table[gammahat[[i, 1, 1]], {i, ngrid}], 32], 4]; 

你可以刪除Table：

FourierDST[SetPrecision[gammahat[[All, 1, 1]], 32], 4]; 

而且，如果你真的，真的需要這個SetPrecision，不能你做的一次伽瑪的計算？

AFAI可以看到，伽瑪計算中使用的所有數字都是確切的。這可能是故意的，但速度很慢。您可以考慮使用近似數字。

編輯
隨着您的最新編輯的完整代碼只是增加一個//N你的第二和第三線切割線時間至少一半沒有太大的降低數值精度。如果我比較res = {gammafirst，gammasecond，AFD}中的所有數字，原始和添加的// N之間的差值爲res1 - res2 // Flatten // Total ==> 1.88267 * 10^-13

刪除所有SetPrecision填充物會加速編碼爲7，結果似乎具有相似的準確性。