function [U,Umax,Xmax] = theBigFish(alpha)
% THEBIGFISH
%   Resoud un tres gros poisson sur un cube d'arete 2.
%
%   Plus serieusement, cette fonction resoud l'equation de Poisson,
%
%       d^2(U)/dx^2 + d^2(U)/dy^2 + d^2(U)/dz^2 = f(x,y,z)
%
%   dans un cube defini entre (-1,-1,-1) et (1,1,1). La fonction f est
%   definie par 
%
%       f(x,y,z) = -10      si 0 <= x <= 1
%                              0 <= y <= 1
%                              0 <= z <= 1
%       f(x,y,z) = 0        sinon
%
%   et la condition aux limites vaut
%
%       U(-1,y,z) = U(1,y,z) = 0
%       U(x,-1,z) = U(x,1,z) = 0
%       U(x,y:-1) = U(x,y,1) = 0
%
%
%   in  : alpha, un parametre lie au maillage du domaine de U
%
%   out : U la solution de l'equation
%            ( size(U)=[4*alpha+1,4*alpha+1,4*alpha+1] )
%         Umax la valeur maximale de U sur son domaine
%         Xmax les coordonnees de ce maximum
%

% Laurent Van Eesbeeck
% 17/12/11
% Version 3 : nettoyage de pcg en la redefinissant en interne

global m M

m=4*alpha+1; % le nombre de noeuds sur une arete du maillage
M=1:m^3; % les indices absolus associes au maillage
h=2/(m-1);

%m=4

% indices des variables sur le bord et a l'interieur du cube + ceux ou le
% laplacien est non-nul et nul
[border,nonborder]=borderindexes;
[nonnul,nul]=findb;

% Ecriture de A
% la matrice du systeme lineaire a resoudre est de la forme
% size(A)=[m^3 m^3]
% et une ligne ne correspondant pas a une condition frontiere est
% 
% indice     i-m^2     i-m   i-1  i i+1   i+m     i+m^2
% ligne  [ - - - 1 - - - 1 - - 1 -6 1 - - 1 - - - 1 - - - ]
% 
% la matrice est donc creuse de forme 7-diagonale.
% Les conditions de bord n'interviendront pas dans cette matrice, il faudra
% donc en supprimer certaines lignes et colonnes. On optimise le temps de
% calcul en prealouant A de taille la plus petite possible. (Definir A est
% l'operation la plus lourde en temps de calcul apres la resolution du
% systeme)
nonborder2=nonborder-nonborder(1)+1; % nonborder2 commence ainsi a 1

fill=ones(nonborder2(end),1); % vecteur de remplissage de la matrice
    % contenant suffisamment d'elements pour considerer toutes les
    % equations des points interieurs
A=spdiags([fill fill fill -6*fill fill fill fill],[-m^2 -m -1 0 1 m m^2],length(fill),length(fill));

% Ecriture de b
% Obtention des indices pour lesquels le laplacien est non-nul et ecriture
% de b
b=[-10*ones(size(nonnul)) zeros(size(nul))];
% tri de b pour le faire correspondre aux indices
[~,idx]=sort([nonnul nul]);
b=b(idx);
b=b*h^2;


% On ne doit pas resoudre l'equation sur son bord puisqu'on sait qu'il vaut
% 0! On enleve donc les lignes et inconnues inutiles du systeme, ce qui
% fera gagner en temps de resolution
b=b(nonborder); b=b';
A=A(nonborder2,nonborder2);


% Resolution du systeme
% A est definie negative... pour pouvoir utiliser pcg, A doit etre definie
% positive, ce qui s'obtient en considerant le systeme sous la forme
% -Ax'=-b.
% La tolerance de la methode est fixee a 1e-4.
% Remarque: une pauvre copie de pcg (nettoyee de toutes ses options
% inutiles dans notre cas) est utilisee.
U1=mypcg(-A,-b,1e-4,400);
%[U1,~]=pcg(-A,-b,1e-4,400);
%U1=A\b;

% reecriture de U en ajoutant le bord du domaine
U=[U1' zeros(size(border))];
[~,idx]=sort([nonborder border]);
U=reshape(U(idx),m,m,m);

% extraction des coordonnees du maximum.
[Umax,I]=max(U(:));
[x,y,z]=ind2sub([m m m],I);
Xmax=[x y z];

% on normalise Xmax sur [-1,1]
Xmax=2*(Xmax-1)/(m-1)-1;


end





function [nonnul,nul]=findb
% vectorise la recherche des indices appartenant au sous-cube
% [0,1]x[0,1]x[0,1] pour lesquels le laplacien est non-nul

% Mon raisonnement est trop long a expliquer, donc "croyez-moi sur parole,
% ma fonction est bonne"

global m M

half=ceil(m^3/2);
d=floor(m/2);

temp1=(half:half+d)';
temp2=temp1*ones(1,d+1)+m*ones(d+1,1)*(0:d);
temp2=temp2(:);
temp3=temp2*ones(1,d+1)+m^2*ones((d+1)^2,1)*(0:d);

nonnul=temp3(:)';
nul=M(~ismember(M,nonnul));
end



function [border,nonborder]=borderindexes
% deux vectorisations differentes pour trouver les indices des variables
% correspondant au bord du domaine et aux point interieurs. La premiere est
% plus rapide pour des grands m, la seconde pour des plus petits m. Apres
% des tests, la valeur seuil de m est approximativement 23

global m M

if m<4
    % pour masquer les defauts de ma vectorisation generique
    fprintf('Reprends ta vie en main, a sert a rien de vouloir resoudre\nune EDP avec un maillage aussi grossier...\n');
    error('Sois un peu raisonnable...')
end
if m>=23
    % border-indexes vectorisation for m>=4
    % just believe me, this works.
    temp1=1:m^2; % temp1 : first slice of the cube

    temp21=m^2+1:m^2+m+1;
    temp221=m^2+2*m+(0:m:m*(m-4));
    temp222=[temp221;(temp221+1)];
    temp22=temp222(:);
    temp23=2*m^2-m:2*m^2;
    temp2=[temp21 temp22' temp23];
    temp2=temp2'*ones(1,m-2)+ones(4*m-4,1)*(0:m-3).*m^2;
    temp2=temp2(:)'; % temp2 : edge of internal slices

    temp3=m^3-m^2+1:m^3; % temp3 : last slice of the cube

    border=[temp1 temp2 temp3]; % indices i tels que u(i) se trouve sur la frontiere
    nonborder=M(~ismember(M,border));
else
    % non-border indexes vectorisation for m>=3
    % like before, believe me.

    temp1=(m^2+m+2:m^2+2*m-1)';

    temp2=temp1*ones(1,m-2) +m*ones(m-2,1)*(0:m-3);
    temp2=temp2(:); % temp2 : interior indexes of the first non-border slice

    nonborder=temp2*ones(1,m-2) +m^2*ones((m-2)^2,1)*(0:m-3);
    nonborder=nonborder(:)';
    border=M(~ismember(M,nonborder));
end
end










function [x] = mypcg(A,b,tol,maxit)
%PCG   Preconditioned Conjugate Gradients Method. Version nettoyee pour ce
%      probleme.
%
%   X = PCG(A,B) attempts to solve the system of linear equations A*X=B for
%   X. The N-by-N coefficient matrix A must be symmetric and positive
%   definite and the right hand side column vector B must have length N.
%
%   X = PCG(A,B,TOL,MAXIT) specifies the maximum number of iterations. If
%   MAXIT is [] then PCG uses the default, min(N,20).

%   Copyright 1984-2010 The MathWorks, Inc.
%   $Revision: 1.18.4.8 $ $Date: 2010/02/25 08:11:56 $

n = length(A);

x = zeros(n,1);

% Set up for the method
tolb = tol*norm(b);                  % Relative tolerance
r = b - A*x;
normr = norm(r);                   % Norm of residual
normr_act = normr;

rho = 1;
moresteps = 0;
maxmsteps = min([floor(n/50),5,n-maxit]);

% loop over maxit iterations (unless convergence or failure)

for ii = 1 : maxit
    z = r;
    
    rho1 = rho;
    rho = r' * z;

    if (ii == 1)
        p = z;
    else
        beta = rho / rho1;
        p = z + beta * p;
    end
    q = A*p;
    pq = p' * q;
    alpha = rho / pq;
        
    x = x + alpha * p;             % form new iterate
    r = r - alpha * q;
    normr = norm(r);
    normr_act = normr;
    
    % check for convergence
    if (normr <= tolb || moresteps)
        r = b - A*x;
        normr_act = norm(r);
        if (normr_act <= tolb)
            break
        else
            moresteps = moresteps + 1;
            if moresteps >= maxmsteps
                break;
            end
        end
    end
end                                % for ii = 1 : maxit

end