function [LAM,LAM_Dir,Err_LAM,Err_Dir_LAM,N] = Eigen_Multigrid_Multi(Max_Level,mesh)

% Compue the second order elliptic eigenvalue problem by the Multigrid 
% method with the conforming finite element method. 
% The program computes the first several eigenpair and gives the error estimates
% between the approximations by the multigrid method and direct
% eigensolvers. 
% INPUT: Max_Level: the maxmum level of meshes
%        Mesh: the input mesh
% OUTPUT: LAM: the eigenvalue approximations by the multilevel correction
% method
%         LAM_Dir: the eigenvalue approximations by the direct eigensolver
%         Dis_1: \|u_h-u_h^{dir}\|_1
%         Dis_0: \|u_h-u_h^{dir}\|_0
%         N: the numbers of elements for the sequence of meshes
%
% Hehu Xie: 2010-12-20
%--------------------------------------------------------------------------

if nargin < 2
    uniformrefine = 3;
    %--------unit square: triangle-----------------------------------------
    node = [0 0; 1 0; 1 1; 0 1];
    elem = [ 1 2 4;
             2 3 4 ];  
    Dirichlet = [1 2;2 3; 3 4; 4 1];
    Neumann=[];
    %----------------------------------------------------------------------
    mesh=Initial_FEM_matlab(node,elem,Dirichlet,Neumann);
    for i=1:uniformrefine
        mesh = refinement(mesh);
    end
end
%--------------------------------------------------------------------------
coord3 = quad_triangle(13);
order = 2;
Base = Lagrange_basis_function(order);
Stiff = Init_Matrix(1,1,Base,Base,[0 1 1], [0 1 1],@Discrete_Form_Stiff,coord3);
%--------------------------------------------------------------------------
Mass = Init_Matrix(1,1,Base,Base,[1 0 0], [1 0 0],@Discrete_Form_Mass,coord3);
%------------- Assemble matrices --------------------------------------------
A = Matrix_Assemble(mesh,Stiff);
M = Matrix_Assemble(mesh,Mass);
%----------------handle the boundary condition-----------------------------
[A,M] = Poission_BoundaryTreatment(mesh,Base,@u_D,A,M);
%-----solve a linear eigenvalue to get the initial solution----------------
[xv,lm,iresult] = sptarn(A,M,19,100);
Num_Eig = length(lm);
xv = xv(:,1:Num_Eig);
LAM(1:Num_Eig,1) = lm(1:Num_Eig);
LAM_Dir(1:Num_Eig,1) = lm(1:Num_Eig);
%-----------to do the multilevel correction method-----
mesh = Init_Multi_Mesh(mesh);
meshs = [mesh,mesh];
clear mesh;
Cur_Level = 2;
%---------Multigrid eigenvalue solution------------------------------------
Num_Correction = 5*ones(Max_Level,1);
%Num_Correction(end) = 5;
N(1) = size(meshs(1).elem,1);
A1 = A;
M1 = M;
while Cur_Level <= Max_Level
    Cur_Level
    %------------propolongation---------------------
    [meshs(1),xv] = refinement(meshs(1),xv,Base);
    N(Cur_Level) = size(meshs(1).elem,1);
    %-----build the stiff matrix ------------------------------------------
    Stiff.R_Mesh = 1;
    Mass.R_Mesh = 1;
    A = Matrix_Assemble(meshs(1),Stiff);
    M = Matrix_Assemble(meshs(1),Mass);
    [A,M] = Poission_BoundaryTreatment(meshs(1),Base,@u_D,A,M);  
    %----Solve the eigenvalue problem by the direct solver-----------------
    [xv_dir,lm_dir,iresult] = sptarn(A,M,19,100);
    lm_dir = lm_dir(1:Num_Eig)
    xv_dir = xv_dir(:,1:Num_Eig);
    [v_max,index] = max(abs(xv_dir));
    xv_dir = xv_dir*diag(1./diag(xv_dir(index,:)));
    %----------Do the correction---------------
    for iter = 1 : Num_Correction(Cur_Level)
        RHS = M*xv*diag(lm);
        [A,RHS] = Poission_BoundaryTreatment(meshs(1),Base,@u_D,A,RHS);
        xv = A\RHS;
        Stiff.R_Mesh = 2;
        Mass.R_Mesh = 2;
        %---------to do the correction second time-------------------------
        ann = xv'*A*xv;
        AHh = Matrix_Assemble(meshs,Stiff);
        RHS = AHh*xv;
        AA = [A1,RHS; RHS',ann];
        RHS = Matrix_Assemble(meshs,Mass)*xv;
        MM = [M1,RHS;RHS',xv'*M*xv];
        [AA,MM] = Poission_BoundaryTreatment(meshs(2),Base,@u_D,AA,MM);
        %---eigenvalue-----------------------------------------------------
        [xv2,lm,iresult] = sptarn(AA,MM,18,100);  
        lm = lm(1:Num_Eig)
        xv2 = xv2(:,1:Num_Eig);
        iter
        xv = Get_Solution_Multilevel(Base,meshs(2),xv2(1:end-Num_Eig,:),Base,...
                                     meshs(1),xv2(end-Num_Eig+1:end,:),xv);  
        [v_max,index] = max(abs(xv));
        v_max = xv(index,1);
        xv = xv*diag(1./diag(xv(index,:)));
    end 
    %---------------------------------------------------------------------- 
    LAM(1:Num_Eig,Cur_Level) = lm
    LAM_Dir(1:Num_Eig,Cur_Level) = lm_dir
    Dis_1(Cur_Level) = 0;
    Dis_0(Cur_Level) = 0;
    
    Cur_Level = Cur_Level + 1;    
end
Dis_LAM = LAM-LAM_Dir;
for i = 1 : Num_Eig
    for j = 1 : Num_Eig
        Exact_LAM((i-1)*Num_Eig+j) = (i^2+j^2)*pi^2;
    end
end
Exact_LAM = sort(Exact_LAM);
for kk = 1 : Num_Eig
    Err_LAM(kk,:) = LAM(kk,:) - Exact_LAM(kk);
    Err_Dir_LAM(kk,:) = LAM_Dir(kk,:) - Exact_LAM(kk);
end

Notations = ['-o';'-*';'-s';'-+';'-x';'-d';'-p';'-v';'-^';'-<';'->'];
Ind = 0;
for kk = 1 : Num_Eig
    loglog(N,Err_LAM(kk,:),Notations(mod(Ind,11)+1,:),...
        N,Err_Dir_LAM(kk,:),Notations(mod(Ind+1,11)+1,:),...
        N,0.8*Err_LAM(kk,2)*N(2)^(order)*N.^(-(order)),'-.',... 
       'LineWidth',3,'markersize',14)
   Ind = Ind + 2;
   hold on
end
hold off
%--------------------------------------------------------------------------
if Num_Eig == 1
    word = {'1-st eigenvalue','1-st direct eigenvalue',['slope= -', num2str(order)]};
end
if Num_Eig >= 2
    word = {'1-st eigenvalue','1-st direct eigenvalue', ['slope= -', num2str(order)],'2-nd eigenvalue',... 
        '2-nd direct eigenvalue',['slope= -', num2str(order)]};
end
for i = 3 : Num_Eig
    word = [word,{[num2str(i) '-th eigenvalue']},{[num2str(i) '-th direct eigenvalue']},{['slope= -', num2str(order)]}];    
end
legend(word);
set(gca, 'FontSize', 14)
title('Convergence behavior of multigrid method for eigenvalue problem','fontsize',18)
xlabel('Number of iterations','fontsize',18)
ylabel('Eigenvalue error','fontsize',18)
toc
%--------------------------------------------------------------------------
%  Given conditions
%--------------------------------------------------------------------------
function v = Exact(x)
v = sin(pi*x(:,1)).*sin(pi*x(:,2));
%-------------------------------------------------------------------------
function v = Exact_Grad(x)
v = pi*[cos(pi*x(:,1)).*sin(pi*x(:,2)),sin(pi*x(:,1)).*cos(pi*x(:,2))];
%-----------------Sub functions-----------------------
function v = Discrete_Form_Stiff(VL,VR,X,Aux_FEM,Matrix)
v = VL(:,1).*VR(:,1) + VL(:,2).*VR(:,2);
%-----------------Mass Matrix---------------------
function v = Discrete_Form_Mass(VL,VR,X,Aux_FEM,Matrix)
v = VL(:,1).*VR(:,1);
%-------------------------------------------------------------------------
function v = Discrete_Form_Rh(VL,VR,X,Aux_FEM,Matrix)
v = VR(:,1).*2*pi^2.*sin(pi*X(:,1)).*sin(pi*X(:,2));
%----------boundary condition----------------------------------------------
function v=u_D(x)
v = sparse(size(x,1),1);
%-------------------------------------------------------------------------
function v = u_N(x)
v = sparse(size(x,1),1);