function P4()
% Checks ThreeDPoisson
clc
clear
randn('seed',0)


% A Small example..
n = [10 20 30]';
f = randn(n(1)*n(2)*n(3),1);
uCVL = ThreeDPoissonCVL(f,n);
uYou = ThreeDPoisson(f,n);
relErrSmall = norm(uCVL-uYou)/norm(uCVL)

% A Big example..
n = [100 200 50]';
f = randn(n(1)*n(2)*n(3),1);
uCVL = ThreeDPoissonCVL(f,n);
uYou = ThreeDPoisson(f,n);
relErrBig = norm(uCVL-uYou)/norm(uCVL)

% -2 if yur code takes longer than a minute to solve the big example

 function u = ThreeDPoissonCVL(f,n)
% Solves (In1 x In2 x Tn3  + In1 x Tn2 x In3 + Tn1 x In2 x In3)u = f
% where u and f are column N vectors with N = n(1)*n(2)*n(3)

% The problem transforms to
% 
% (In1 x In2 x Dn3  + In1 x Dn2 x In3 +  Dn1 x In2 x In3)* utilde = ftilde
%
% where
%         ftilde = (inv(Sn1) x inv(Sn2) x inv(Sn3))*f
%                = (8/((n1+1)(n2+1)(n3+1)))*(Sn1 x Sn2 x Sn3)*f
% and 
%              u = (Sn1 x Sn2 x Sn3)*utilde
%
% See the function ThreeSProd for carrying out (Sn1 x Sn2 x Sn3)*vector

% Compute ftilde...

  ftilde = (8/((n(1)+1)*(n(2)+1)*(n(3)+1)))*ThreeSprod(n,f);
  
% Now solve the diagonal system. Note that if
%    D = In1 x In2 x Dn3  + In1 x Dn2 x In3 +  Dn1 x In2 x In3
% then d = diag(d) where
%    d = ones(n1*n2,1) x d3  + ones(n1,1) x d2 x ones(n3,1) + d1 x ones(n2*n3,1)
% with
   d1 = 4*sin(pi*(1:n(1))'/(2*(n(1)+1))).^2;
   d2 = 4*sin(pi*(1:n(2))'/(2*(n(2)+1))).^2;
   d3 = 4*sin(pi*(1:n(3))'/(2*(n(3)+1))).^2;
   d = kron(ones(n(1)*n(2),1),d3) + kron(kron(ones(n(1),1),d2),ones(n(3),1))+...
       kron(d1,ones(n(2)*n(3),1));
% Solving a diagonal system is pointwise vector division...
   utilde = ftilde ./d;
   
% Finally, compute u from utilde...
   u = ThreeSprod(n,utilde);
   
   function z = ThreeSprod(n,y)
 % n is an integer 3-vector [n1 n2 n3].
 % y is a column vector with length n1*n2*n3 and 
 % z = (Sn1 x Sn2 x Sn3)y. 
 
 % The implementation involves repeated use of the following Kronecker product fact:
 %       If A is m1-by-m1 and B is m2-by-m2 and y = (A x B)x, then
 %       reshape(y,m2,m1) = B*reshape(x,m2,m1)*A'

 % Applying this to our problem gives:
 %        reshape(z,n2*n3,n1) = (Sn2 x Sn3)*reshape(y,n2*n3,n1)*Sn1
 
 % Overwrite
 Y = reshape(y,n(2)*n(3),n(1));
 % with Y*Sn1...
 for i=1:n(2)*n(3)
     Y(i,:) = MySineCVL(Y(i,:)')';
 end
 % Now overwrite Y  with (Sn2 x Sn3)*Y. We do this one column at a time..
 for j=1:n(1)
     %  Reshape Y(:,j)...
     Yj = reshape(Y(:,j),n(3),n(2));
     % Now overwrite Yj with  Sn3*Yj*Sn2
     % First overwrite Yj with Sn3*Yj
     for i=1:n(2)
         Yj(:,i) = MySineCVL(Yj(:,i));
     end
     % Now overwrite Yj with Yj*Sn2
     for i=1:n(3)
         Yj(i,:) = MySineCVL(Yj(i,:)')';
     end
     % Reshape Yj into the jth column of Y...
     Y(:,j) = reshape(Yj,n(2)*n(3),1);
 end
 % Y houses the solution. Reshape it into z..
 z = reshape(Y,n(1)*n(2)*n(3),1);
 
 
 function y = MySineCVL(x)
% x is a real column n-vector
% y = Sn*x where Sn is then-by-n  matrix defined by
%               Sn(k,j) = sin(k*j*pi/(n+1))
  n = length(x);
  xTilde = [0; x;0; -x(n:-1:1)];
  yTilde = fft(xTilde);
  y = real(sqrt(-1)*yTilde(2:n+1)/2);
     

   
   