Solve Poisson on a square using Chebyshev collocation and sum factorization

logg.jl

## Julia implementation of solving Poisson on a square grid using Chebyshev collocation
# julia> load("logg.jl")
# julia> loggProfile(10000)
# error: 3.310154100518143e-7
# Average time (us): 132.39421844482422

function logg(m)
    (D,xb) = cheb(m+1) # Chebyshev differentiation matrix and nodes
    D = 2*D; xb = 0.5 + 0.5*xb   # Remap to interval [0,1]
    D2 = -(D*D)[2:end-1,2:end-1] # Negative 1D Laplacian, with Dirichlet BCs
    xd = xb[2:end-1]             # Remove Dirichlet BCs
    sinx = sin(pi*xd'')
    xexact = kron(sinx', sinx)   # sin(pi*x)*sin(pi*y)
    b = 2*pi^2*xexact
    # Compute action of inverse of L using sum factorization
    # See Lynch, Rice, and Thomas (1964)
    (Lambda, R) = eig(D2)
    e = ones(m)
    Diag = kron(Lambda,e) + kron(e,Lambda)
    iDiag = 1./Diag
    Rinv = inv(R)
    # L = kron(D2,eye(m)) + kron(eye(m),D2)
    # Linv = kron(R,R) * diagm(reshape(iDiag,m^2,1)) * kron(Rinv,Rinv)
    x = R * (iDiag .* (Rinv*b*Rinv')) * R' # x = Linv * b  (with reshaping)
    norm(reshape(x - xexact,m^2,1))        # L2 norm of the vector
end

function cheb(N)
    x = cos(pi*(0:N)/N)
    c = ones(N+1); c[1] = 2; c[end] = 2; c[2:2:] = -c[2:2:]
    X = repmat(x,1,N+1)
    dX = X-X'
    D = (c*(1/c')) ./ (dX+eye(N+1))  # off-diagonal entries
    D = D - diagm(sum(D,2))          # diagonal entries
    (D, x)
end

function loggProfile(count)
    gc()                        # improves timing reproducibility
    tic()
    l2error = 0
    for i=1:count
        l2error = logg(7)
    end
    avgtime = toq()/count
    println("L2 error: ", l2error)
    println("Average time (us): ", 1e6*avgtime)
end

logg.m

function error = logg(m)
    [D,xb] = Cheb(m+1); % Chebyshev differentiation matrix and nodes
    D = 2*D; xb = 0.5 + 0.5*xb;   % Remap to interval [0,1]
    D2 = -(D*D)(2:end-1,2:end-1); % Negative 1D Laplacian
    xd = xb(2:end-1);             % Remove Dirichlet BCs
    sinx = sin(pi*xd);
    xexact = kron(sinx', sinx);   % sin(pi*x)*sin(pi*y)
    b = 2*pi^2*xexact;
    % Compute action of inverse of L using sum factorization
    % See Lynch, Rice, and Thomas (1964)
    [R, Lambda] = eig(D2);
    Lambda = diag(Lambda);
    e = ones(m,1);
    Diag = kron(Lambda,e') + kron(e,Lambda');
    iDiag = 1./Diag;
    Rinv = inv(R);
    % L = kron(D2,eye(m)) + kron(eye(m),D2);
    % Linv = kron(R, R) * diag(reshape(iDiag,m^2,1)) * kron(Rinv, Rinv);
    x = R * (iDiag .* (Rinv*b*Rinv')) * R'; % x = Linv * b  (with reshaping)
    error = norm(x - xexact);
end

function [D,x] = Cheb(N)
    if N==0, D=0; x=1; return, end
    x = cos(pi*(0:N)/N)';
    c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
    X = repmat(x,1,N+1);
    dX = X-X';
    D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
    D  = D - diag(sum(D'));                 % diagonal entries
end