加权Jacobi的迭代优化


#1

我自己写了一个加权 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


库函数的正确用法?
#2

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

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


#3


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


#4

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

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


#5

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


#6

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


#7

:astonished: 为啥?
broadcast 看着很简洁 :sweat_smile:


#8

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


#9

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


#10

刚刚又看了一下,发现了问题。你把不该加@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)


#11

受教了,特别感谢!
不过我写的确实是 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

然后 uhuhJ 加权得到最终的结果: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

更好点?

我需要测试一下?


#12

我有点模糊的地方是,在 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 循环的写法上还挺多的,一般是怎么确定哪一种比较好呢?根据自己的经验?


#13

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

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

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会更慢。详情可以看 性能优化常见问题

一般情况下我不怎么喜欢加 @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。看不懂的话可以问我。


#14

这才是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)


#15

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


#16

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

赞,赞

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


#17

你最近在学FDM?


#18

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


#19

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


#20

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