% -------------------------------------------------------------------------
% Loads a finite element matrix 
% for the convection-diffusion-reaction problem
% d/dx(epsx du/dx) + d/dy(epsy du/dy) + ...
% v1 du/dx + v2 du/dy + c u = f;
% Parameters: 
% epsx and epsy may be constants and used to describe anisotropy
% epsx and epsy may depend on (x,y) and describe discontinuous coefficients
% (v1,v2) describes the velocity field
% c is the reaction term and may simulate the inverse of the time step
%   for time-dependent problems
% -----------------------------------------------------------------------
% For spd matrices the program compares three different multigird solvers
% -- Geometric MG from J. Demmel
% -- Algebraic MG fro the HSL library
% -- Algebraix MG (AGMG), developed by Yvan Notay
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%figure('renderer','painters'); % if OpelGL crashes

% addpath('/home/maya/student/NLA/Applications/Multigrid_Demmel')
% addpath('/home/maya/student/NLA/AGMG_Yvan_Notay/AGMG_3.2.0-aca/Matlab')
% addpath /home/maya/matlab/mi20_matlab
% javaaddpath /home/maya/matlab/mi20_matlab
addpath('/home/maya/student/NLA/Applications/Multigrid_Demmel')
addpath('/home/maya/student/NLA/AGMG_Yvan_Notay/AGMG_3.2.0-aca/Matlab')
addpath('/home/maya/matlab/AMG_HSL/mi20_matlab')
javaaddpath /home/maya/matlab/AMG_HSL/mi20_matlab
flag_Laplace=1;
% load a Laplace matrix K,sol,rhs,Node
load('Laplace_FEM_1089.mat')
n=size(K,1);
switch n
    case 1089
        levels=5;
    case 4225
        levels=6;
    case 16641
        levels=7;
    case 66049
        levels=8;
end
tic, sol_ex=K\rhs;tim_direct=toc;

% Run J. Demmel' test
tic,makemgdemo_time, tim_preMG=toc;
tic,testfmgv_time,   tim_exeMG=toc;
tim_totMG = tim_preMG+tim_exeMG; 

% Run AGMG 
tic,agmg(K,[],[],[],[],[],[],1);tim_preAGMG=toc;
tic,[x1,fl,relres,it_agmg,resv1] = agmg(K,rhs,[],[],[],[],[],2);tim_exeAGMG=toc;
tim_totAGMG = tim_preAGMG+tim_exeAGMG; 

% Run mi20
mi20_startup;
control = mi20_control();
tol = 1e-6; maxit=n;
tic,inform = mi20_setup(K,control);tim_preMI20=toc;
tic,[x2,flag,relres,it_mi20,resv2] = pcg(K,rhs,tol,maxit,@mi20_precondition);
mi20_finalize
tim_exeMI20=toc;
tim_totMI20 = tim_preMI20+tim_exeMI20; 

disp(['PCG-AMG  iterations: ' int2str(it_agmg)])
disp(['PCG-mi20 iterations: ' int2str(it_mi20)])

disp(['Direct solver time: ' num2str(tim_direct)])
disp(' ')
disp(['MG   preproc. time: ' num2str(tim_preMG)])
disp(['     solution time: ' num2str(tim_exeMG)])
disp(['     total time:    ' num2str(tim_totMG)])
disp(' ')
disp(['AGMG preproc. time: ' num2str(tim_preAGMG)])
disp(['     solution time: ' num2str(tim_exeAGMG)])
disp(['     total time:    ' num2str(tim_totAGMG)])
disp(' ')
disp(['MI20 preproc. time: ' num2str(tim_preMI20)])
disp(['     solution time: ' num2str(tim_exeMI20)])
disp(['     total time:    ' num2str(tim_totMI20)])
disp(' ')

figure(2)
semilogy(resv1./resv1(1),'rd-')
semilogy(resv2./resv2(1),'mp-')
if flag_Laplace==1,
   legend('MG','AMG','MI20')
end
title('residual(k)'), 
xlabel('iteration number k'), grid