使用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，结果似乎具有相似的准确性。