function [U,Umax,Xmax]=theBigFish(alpha)
% A 3D Fast Poisson Solver for the equation Laplace(u)+f=0 in a cube with the
% Dirichlet boundary condition.
%
% This uses the FFTW library to choose the fastest fourrier transfrom algorithm to
% perform matrix multiplications. Indeed, we know the eigen vectors and
% values of this particular problem so we can use fft to speed up the
% performances.
% 
% DON'T FORGET TO RUN IT TWICE IF YOU CHANGE THE PROBLEM SIZE. See fftw
% below)
%
% Author: François HEREMANS adapted from Weiguo Gao and Yuan Cao and with help from
% http://www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html and
% http://ocw.mit.edu/courses/mathematics/18-086-mathematical-methods-for-en
% gineers-ii-spring-2006/readings/am35.pdf and
% http://www.fftw.org/
% and open lectures from professor Gilbert Strand (MIT)
%
% Version: 2.0 (18/12/2011)
%
% theBigFish(alpha):
% - alpha is the multiplicating factor for the number of points.
% - U is the 3D result tensor
% - Umax is the maximum value found 
% - Xmax is the coordinates of that maximum. (if multiple, returns one of
% them)

%% SETTINGS

fftw('planner', 'exhaustive') % Dertermines the fastest algorithm for 
%                            the current problem size. Only lauched at first
%                            call of the function. Don't forget to lauch
%                            the function twice when you change the problem
%                            size.

n = 4*alpha+1;
h=2/(n-1);
h2 = h^2;
hx2=h2; % just in case we want to change the form of the cube.
hy2=h2;
hz2=h2;
f=zeros(n,n,n);
f(round(n/2):n,round(n/2):n,round(n/2):n)=10; % Initial condition
[nx ny nz]=size(f);
U=zeros(n,n,n);

%% BOUNDARY CONDITIONS

f(2,:,:)=f(2,:,:)+U(1,:,:)/hx2;
f(nx-1,:,:)=f(nx-1,:,:)+U(n,:,:)/hx2;

f(:,2,:)=f(:,2,:)+U(:,1,:)/hy2;
f(:,ny-1,:)=f(:,ny-1,:)+U(:,n,:)/hy2;

f(:,:,2)=f(:,:,2)+U(:,:,1)/hz2;
f(:,:,nz-1)=f(:,:,nz-1)+U(:,:,n)/hz2;

u=f(2:nx-1,2:ny-1,2:nz-1);
clear f

nx=nx-2;
ny=ny-2;
nz=nz-2;

%% COMPUTE FFT

%inverse discrete sine transform in two directions.
u=my3diDST(u);
u=permute(my3diDST(permute(u,[2 1 3])),[2 1 3]);

%solve the tridiagonal linear equations in one direction.
lambdax=4/hx2*(sin((1:nx)*pi/2/(nx+1))).^2;
lambday=4/hy2*(sin((1:ny)*pi/2/(ny+1))).^2;
trid=sparse([1:nz-1 2:nz 1:nz]',[2:nz 1:nz-1 1:nz]',-ones(3*nz-2,1)/hz2);
for i=1:nx
    for j=1:ny
        trid(1:nz+1:nz*nz)=lambdax(i)+lambday(j)+2/hz2;
        u(i,j,:)=trid\reshape(u(i,j,:),nz,1);
    end
end

%call discrete sine transform to get the solution.
u=my3dDST(u);
u=permute(my3dDST(permute(u,[2 1 3])),[2 1 3]);

U(2:n-1,2:n-1,2:n-1)=u;
M2=max(max(U));
[Umax,I3]=max(M2);
X=-1+h*(I3-1);
Xmax=[X,X,X];

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y=my3diDST(x)
%Discrete sin transform for 3D arrays. (c.f. dst)
y=2/(size(x,1)+1)*my3dDST(x);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y=my3dDST(x)
%Discrete sin transform for 3D arrays. (c.f. dst)
n=size(x,1);
sizex=size(x);
y=zeros(sizex);
xx=zeros([2*(n+1),sizex(2:3)]);
xx(2:n+1,:,:)=x;
xx(n+3:2*(n+1),:,:)=-x(n:-1:1,:,:);
xx=imag(fft(xx));
y=xx(2:n+1,:,:)/(-2);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% JUST IN CASE, THE VERY SLOW VERSION USING THE TEACHER'S METHOD

% function [U,Umax,Xmax] = theBigFish2(alpha)
% %theBigFish2 Joli programme qui resout
% %  une equation de poisson par des
% %  differences finies sur un cube
% %
% %   Programme dedie a Gregoire et Jean-Francois
% %   
% 
% %close all;
% nx = 4*alpha+1;
% ny = 4*alpha+1;
% nz = 4*alpha+1;
% n = nx*ny*nz;
% h = 2/(ny-1);
% 
% A = diag(ones(n,1));
% B = zeros(n,1);
% for i = 2:nx-1
%     for j = 2:ny-1
%         for k=2:nz-1
%             index = i + (j-1)*nx +(k-1)*nx^2;
%             A(index,index) = 6;
%             A(index,index+1) = -1;
%             A(index,index-1) = -1;
%             A(index,index+nx) = -1;
%             A(index,index-nx) = -1;
%             A(index,index-nx^2) = -1;
%             A(index,index+nx^2) = -1;
%             B(index) = 10*(i>=nx/2 && j>=ny/2 && k>=nz/2);
%         end
%     end
% end
% 
% A2 = sparse(A)./(h*h);
% 
% U = A2\B;
% 
% U = reshape(U,nx,ny,nz);
% 
% [M1,I1]=max(U);
% [M2,I2]=max(M1);
% [Umax,I3]=max(M2);
% 
% Xmax=[-1+h*(I1(M1==Umax)-1),-1+h*(I2(M2==Umax)-1),-1+h*(I3-1)];
% 
% end