加权Jacobi的迭代优化

mifan · 2018 年10 月 17 日 13:00

我自己写了一个加权 Jacobi 的迭代求法，可以帮我看下还可以再优化吗？(目前只是 2-dim 中，dim 可以忽略)

function solvedByItSolvers!(::weightedJaboci{dim}, uh::Vector{Float64}, fh::Vector{Float64}, Ninter::Int, weight::Float64 = 0.5) where {dim}
    uh = reshape(uh, Ninter+1, Ninter+1)
    fh = reshape(fh, Ninter+1, Ninter+1)

    # update Jacobi iteration
    uhJ = uh
    for ii = 2:Ninter
        @simd for jj = 2:Ninter
            @inbounds uhJ[ii,jj] = ( fh[ii,jj] + uh[ii-1,jj] + uh[ii+1,jj] + uh[ii,jj-1] + uh[ii,jj+1] )/4
        end
    end

    uh = weight.*uh .+ (1-weight).*uhJ
    return uh[:]
end

woclass · 2018 年10 月 18 日 11:09

能具体说下是什么算法么，我查 Jacobi 迭代 查出来是解线性方程组的那个，加上加权并不能查到有用的东西。

另外最好给一下测试数据。

Scheme · 2018 年10 月 18 日 12:58

楼主应该就是在解拉普拉斯方程？测试数据的话，画一下就好了，大概。虽然还没怎么仔细看，不过我觉得代码没什么问题。把后面的uh[:]应该写成uh或者vec(uh)？broadcast可能没有直接写for循环快。

mifan · 2018 年10 月 18 日 13:11

特别感谢，我是求解 Au = f，A 为对称正定，大概可以参考(其实也是老师留的作业)，Programming of Finite Difference Methods 第5页，我有 matlab 程序，但是学了几天 Julia，现在想转 Julia 实现，所以又把算法用 Julia 写了一下。

其实我是想知道，上面的循环中还有什么优化的。

mifan · 2018 年10 月 18 日 13:15

我可以整理一下贴出比较完整的代码

Scheme · 2018 年10 月 18 日 13:17

自己profile一下就知道哪里慢了。优化的话试试我之前说的，或者直接return nothing。

mifan · 2018 年10 月 18 日 13:21

为啥？
broadcast 看着很简洁

mifan · 2018 年10 月 18 日 13:36

恩，优化了一些，非常感谢

Scheme · 2018 年10 月 18 日 13:39

emmmm，简洁不代表快。因为Julia不能告诉LLVM数组的alias信息，目前。

Scheme · 2018 年10 月 19 日 03:36

刚刚又看了一下，发现了问题。你把不该加@simd的地方加@simd了。去掉以后应该更快。

而且你代码写的也不对啊。你写的那个是Gauss-Seidel，不应该有加权的。我把你的代码稍微改了一下。

using LinearAlgebra
using PyPlot
len = 300
x = y = range(-2pi,stop=2pi,length=len)
uh = [sin(xx^2)+cos(yy)^2 for xx in x, yy in y]
fh = [sin(xx)+cos(yy) for xx in x, yy in y]
function solve_poisson!(uh, fh, iter=1000)
    @inbounds begin
        n = LinearAlgebra.checksquare(uh)
        II = LinearIndices((n, n))
        for _ in 1:iter, ii in 2:n-1, jj in 2:n-1
            uh[II[ii,jj]] = ( fh[II[ii,jj]] + uh[II[ii-1,jj]] + uh[II[ii+1,jj]] + uh[II[ii,jj-1]] + uh[II[ii,jj+1]] )/4
        end
        return uh
    end
end
plt[:subplot](2, 1, 1)
contour(x, y, uh)
plt[:title]("Initial Condition")
solve_poisson!(uh, fh)
plt[:subplot](2, 1, 2)
contour(x, y, uh)
plt[:title]("After 1000 Gauss-Seidel Iterations")
plt[:savefig]("Poisson.png", dpi=600)

mifan · 2018 年10 月 19 日 07:08

受教了，特别感谢！
不过我写的确实是 weighted Jacobi，最新的结果，我用 uhJ 表示

uhJ[ii,jj] = ( fh[ii,jj] + uh[ii-1,jj] + uh[ii+1,jj] + uh[ii,jj-1] + uh[ii,jj+1] )/4

然后 uh 与 uhJ 加权得到最终的结果：uh = weight.*uh .+ (1-weight).*uhJ，

这都不重要，那是

for ii in 2:n-1, jj in 2:n-1
    #---
end

这种写法快，还是

for ii in 2:n-1
    @simd for jj in 2:n-1
            #---
    end
end

更好点？

我需要测试一下？

mifan · 2018 年10 月 19 日 07:19

我有点模糊的地方是，在 matlab 中，全部向量化，如

ii = 2:n-1
jj = 2:n-1
uhJ[ii,jj] = ( fh[ii,jj] + uh[ii-1,jj] + uh[ii+1,jj] + uh[ii,jj-1] + uh[ii,jj+1] )/4

这差不多就是很好的性能了。
julia 不用向量化，但是 for 循环的写法上还挺多的，一般是怎么确定哪一种比较好呢？根据自己的经验？

Scheme · 2018 年10 月 19 日 16:23

你写的根本不是weighted Jacobi，uhJ = uh 以后，uhJ和uh全等。

julia> uh = zeros(5);

julia> uhJ = uh;

julia> pointer(uhJ)
Ptr{Float64} @0x00007fdd153fe7d0

julia> pointer(uh)
Ptr{Float64} @0x00007fdd153fe7d0

julia> uh[1] = 1000
1000

julia> uhJ
5-element Array{Float64,1}:
 1000.0
    0.0
    0.0
    0.0
    0.0

你改变uh以后，uhJ也改变了。因为两个数组alias，所以下一个迭代要等上一个迭代完成以后才能运行，进而不能用 @simd。

mifan:

这都不重要，那是
for ii in 2:n-1, jj in 2:n-1
    #---
end
这种写法快，还是
for ii in 2:n-1
    @simd for jj in 2:n-1
            #---
    end
end
更好点？

我需要测试一下？

for ii in 2:n-1, jj in 2:n-1
    #---
end

和

for ii in 2:n-1
  for jj in 2:n-1
    # ---
  end
end

等价。

就像之前说的，你写的两个数组完全一样，所以实际上运行的是Gauss-Seidel迭代。因为Gauss-Seidel迭代是有顺序的，既下一个迭代依赖上一个迭代，所以用@simd会更慢。详情可以看性能优化常见问题 - #2，来自 Scheme 。

一般情况下我不怎么喜欢加 @simd 因为大部分时间它不会让你的for循环快多少，但是如果你没用好那么一定变慢，或者给出错误结果。做为一个用户，考虑这种超级优化没什么意义。

是的，你写的Gauss-Seidel迭代不但可以更快，而且可以快很多。只是，要付出的代价太大了。你可以考虑做些更有趣的事情。解Poisson’s equation要看东西南北四个点，在column major的矩阵上南北两点相邻，东西两点不相邻，所以读取东西两点要用的时间长一些。你可以想想怎么迭代或者如何重新排列数组，让它只读取在内存相邻的元素。你还可以想想如何写个buffer来打破一些前后依赖，等等。。。这是一般优化这种循环的思想，而不是想怎么写for循环本身更快。你优化思路本身就有问题。我在线性代数教程1: 矩阵乘法 & 性能优化提到过一点优化循环。当然，重新排列Gauss-Seidel迭代可比矩阵乘法难多了。如果你真的非常感兴趣，可以按照 http://people.csail.mit.edu/jrk/halide-pldi13.pdf 的3.1那个部分试一试，他的例子是照片模糊处理。解Poisson’s equation和照片模糊处理要用的循环几乎完全一样。不过你也不用愁，我认识一个人他正在给Julia写这种优化器。到时候你直接用现成的优化器就好了。

如果你要优化的是加权Jacobi迭代，可以用 @simd，因为下一个迭代完全不依赖于上一个迭代。我在下面写了正确的加权Jacobi迭代写法。你可以自己试试。

我就不说怎么优化Jacobi迭代了，它收敛比Gauss-Seidel迭代慢。你要是能优化Gauss-Seidel迭代，就没必要牵挂Jacobi迭代了。而且简单情况下，加权Jacobi迭代还用Gauss-Seidel的两倍内存。你从下面的例子就能看出来Jacobi迭代远不如Gauss-Seidel迭代。当然，数学上你也可以证明它更慢。如果你好奇怎么证明的话，你可以看看 A. Ralston and P. Rabinowitz, A First Course in Numerical Analysis , 2nd edition, McGraw-Hill, New York, 1978. 的445页，Theorem 9.3。看不懂的话可以问我。

Scheme · 2018 年10 月 19 日 16:28

这才是weighted Jacobi，肉眼都能看出结果不一样。

using LinearAlgebra
using PyPlot
len = 300
x = y = range(-2pi,stop=2pi,length=len)
uh = [sin(xx^2)+cos(yy)^2 for xx in x, yy in y]
fh = [sin(xx)+cos(yy) for xx in x, yy in y]
function solve_poisson(uh, fh, iter=1000, weight=0.5)
    @inbounds begin
        n = LinearAlgebra.checksquare(uh)
        uhJ = zero(uh)
        II = LinearIndices((n, n))
        for _ in 1:iter, ii in 2:n-1, jj in 2:n-1
            uhJ[II[ii,jj]] = ( fh[II[ii,jj]] + uh[II[ii-1,jj]] + uh[II[ii+1,jj]] + uh[II[ii,jj-1]] + uh[II[ii,jj+1]] )/4
        end
        for ii in eachindex(uhJ)
            uhJ[ii] = weight*uh[ii] + (1-weight)*uhJ[ii]
        end
        return uhJ
    end
end
plt[:subplot](2, 1, 1)
contour(x, y, uh)
plt[:title]("Initial Condition")
uhJ = solve_poisson(uh, fh)
plt[:subplot](2, 1, 2)
contour(x, y, uhJ)
plt[:title]("After 1000 weighted Jacobi Iterations")
plt[:savefig]("Poisson.png", dpi=400)

Scheme · 2018 年10 月 19 日 16:35

我再过一段时间写一点解PDE的教程好了。

mifan · 2018 年10 月 20 日 03:19

晕了，还是 matlab 的习惯，完全没意识到，感谢！

赞，赞

还有特别感谢这么详细的解答。

Scheme · 2018 年10 月 20 日 04:11

你最近在学FDM？

mifan · 2018 年10 月 20 日 05:15

算是学习吧，我其实做 FEM, DG 的(也在考虑将自己之前写的 matlab 程序转到 julia 上)，只是最近老师让学 multigrid 所以看 FDM ，并试着用 julia 写，也算是一并学习了 julia.

Scheme · 2018 年10 月 20 日 16:47

要是写PDE教程的话，你想看什么方程和什么算法？

mifan · 2018 年10 月 21 日 02:24

DG 求解 Poisson？我看了 JuAFEM 的包，打算后面写个 DG 的