Some iterative solvers do not accept subarray argument?

General Usage
1

Hi,
I am surprised that IterativeSolvers.gauss_seidel can take views for vector arguments while jacobi, sor and ssor cannot. Any reason for this ?


function test()
    nx=ny=10
    dx=dy=2
    A=poissonSparse(nx,ny,dx,dy)
    x2D=ones(nx,ny)
    y2D=ones(nx,ny)
    xv=view(x2D,1:nx*ny)
    yv=view(y2D,1:nx*ny)
    h=1/nx
    w=2/(1+sin(π*h))

    IterativeSolvers.gauss_seidel!(xv,A,yv,maxiter=100) #OK

    # IterativeSolvers.jacobi!(xv,A,yv,maxiter=100) #KO
    # IterativeSolvers.ssor!(xv,A,yv,w,maxiter=100) #KO
    # IterativeSolvers.sor!(xv,A,yv,w,maxiter=100)  #KO

    #OK if I use 1D vectors (not views)
    x1D=Vector{Float64}(undef,nx*ny)
    y1D=Vector{Float64}(undef,nx*ny)
    copyto!(x1D,xv)
    copyto!(y1D,yv)

    IterativeSolvers.jacobi!(x1D,A,y1D,maxiter=100) #OK
    IterativeSolvers.ssor!(x1D,A,y1D,w,maxiter=100) #OK
    IterativeSolvers.sor!(x1D,A,y1D,w,maxiter=100)  #OK
end

Here is the

MWE 
using LinearAlgebra
using IterativeSolvers
using SparseArrays

function l1dpαI(nx,dx,α)
    dxm2=1.0/(dx^2)
    D=ones(nx)
    D .*= 2(1+α)dxm2
    U=ones(nx-1)
    U .*= -1dxm2
    SymTridiagonal(D,U)
end

function poissonSparse(nx,ny,dx,dy)
    Lx=l1dpαI(nx,dx,0.0)
    Ly=l1dpαI(ny,dy,0.0)

    sp=spzeros(Float64,nx*ny,nx*ny)
    for i=1:nx
        for j=1:ny
            I=i+nx*(j-1)

            sp[I,I]=Lx[i,i]
            i>1  && (sp[I,I-1]=Lx[i,i-1])
            i<nx && (sp[I,I+1]=Lx[i,i+1])
        end
    end
    for i=1:nx
        for j=1:ny
            I=i+nx*(j-1)

            sp[I,I]+=Ly[j,j]
            j<ny && (sp[I,I+ny]+=Ly[j,j+1])
            j>1 &&  (sp[I,I-ny]+=Ly[j,j-1])
        end
    end
    sp
end


function test()
    nx=ny=10
    dx=dy=2
    A=poissonSparse(nx,ny,dx,dy)
    x2D=ones(nx,ny)
    y2D=ones(nx,ny)
    xv=view(x2D,1:nx*ny)
    yv=view(y2D,1:nx*ny)
    h=1/nx
    w=2/(1+sin(π*h))

    IterativeSolvers.gauss_seidel!(xv,A,yv,maxiter=100) #OK

    # IterativeSolvers.jacobi!(xv,A,yv,maxiter=100) #OK
    # IterativeSolvers.ssor!(xv,A,yv,w,maxiter=100) #OK
    # IterativeSolvers.sor!(xv,A,yv,w,maxiter=100)  #OK

    x1D=Vector{Float64}(undef,nx*ny)
    y1D=Vector{Float64}(undef,nx*ny)
    copyto!(x1D,xv)
    copyto!(y1D,yv)

    IterativeSolvers.jacobi!(x1D,A,y1D,maxiter=100) #KO
    IterativeSolvers.ssor!(x1D,A,y1D,w,maxiter=100) #KO
    IterativeSolvers.sor!(x1D,A,y1D,w,maxiter=100)  #KO
end

test()
1 Like