function fast_poisson

m = 2047;
h = 1/(m+1);
% test matrix, size m^2 x m^2
A = -gallery('poisson',m)/h^2;
%
% A has a very special structure, which can be exploited using FFT
%
% Generate a random right-hand-side
b = randn(m*m,1);
% Solve the system in Matlab (more-or-less Gaussian elimination)
tic,
  x_ref = A\b;
toc;

tic,
%
% Now use FFT-based solver, so calles Fast Poisson Solver
%
% S=Q*D_S*Q, T=Q*D_T*Q
D_T=repmat(1/h^2,m,1);
D_S=-4/h^2+2/h^2*cos(pi*h*(1:m)');

c=zeros(m*m,1);
% c := Q*b
for i=1:m,
  row_ind    = (i-1)*m + (1:m);
  c(row_ind) = Q_mult(m,b(row_ind));
end
y=zeros(m*m,1);
% y := [D_T D_S D_T]^{-1} C
for j=1:m,
  Gamma_j    = spdiags(repmat([D_T(j), D_S(j), D_T(j)],m,1), -1:1, m, m);
  col_ind    = j:m:(m*m);
  y(col_ind) = Gamma_j\c(col_ind);
end
x=zeros(m*m,1);
% x := Q*y
for i=1:m,
  row_ind    = (i-1)*m + (1:m);
  x(row_ind) = Q_mult(m,y(row_ind));
end
toc;

% calculate the difference:
rel_error = norm(x-x_ref)/norm(x_ref)

function Qc=Q_mult(m,f)
% Use FFT instead of multiplying with the symmetric orthogonal
% matrix Q_ij = sqrt(2*h) * sin(pi*h*i*j)
h=1/(m+1);
% Take 2*m+2-long fft of [0;f;0;0;0;...;0]
w=fft([0;-sqrt(2*h)*f],2*m+2);
Qc=imag(w(2:m+1));