function [LAM,Err,LAM_Dir,Err_Dir,N] = Eigen_Multigrid(Max_Level,mesh)

% Compue the second order elliptic eigenvalue problem by the Multilevel 
% correction method with the conforming finite element method.
% The program computes the first 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(7);
order = 1;
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,40);
xv = xv(:,1);
LAM = lm(1);
LAM_Dir = lm(1);
%-----------to do the multilevel correction method-----
mesh = Init_Multi_Mesh(mesh);
meshs = [mesh,mesh];
clear mesh;
Cur_Level = 2;
%---------Multigrid eigenvalue solution------------------------------------
%---set the iteration times in each level of space-------------------------
Num_Correction = ones(Max_Level,1);
Num_Correction(end) = 1;
Num_Eig = 1;
N(1) = size(meshs(1).elem,1);
A1 = A;
M1 = M;
while Cur_Level <= Max_Level
    Cur_Level
    %------------propolongation---------------------
    [meshs(1),xv] = refinement_multi(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);
    %---Compute the eigenpair approximation by the direct eigensolver------
    [xv_dir,lm_dir,iresult] = sptarn(A,M,19,40);
    lm_dir = lm_dir(1)
    xv_dir = xv_dir(:,1);
    [v_max,index] = max(abs(xv_dir(:,1)));
    v_max = xv_dir(index,1);
    xv_dir = xv_dir/v_max;
    %--------------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);
        %-- solve the eigenvalue problem ----------------------------------
        [xv2,lm,iresult] = sptarn(AA,MM,18,100);        
        lm = lm(1)
        xv2 = xv2(:,1);
        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(:,1)));
        v_max = xv(index,1);
        xv = xv/v_max;
    end 
    %---------------------------------------------------------------------- 
    LAM(Cur_Level) = lm;
    LAM_Dir(Cur_Level) = lm_dir;
    Err(Cur_Level) = error_estimate(meshs(1),Base,coord3,xv,@Exact_Grad,'gradient');
    Err_Dir(Cur_Level) = error_estimate(meshs(1),Base,coord3,xv_dir,@Exact_Grad,'gradient');              
    Cur_Level = Cur_Level + 1;    
end
%----plot the convergence information--------------------------------------
loglog(N,Err,'-o',N,Err_Dir,'-*',N, 0.75*Err(2)*N(2)^((1+order)/2)*N.^(-(order+1)/2),'-.',...
       N, 0.75*Dis_0(2)*N(2)^((order+1)/2)*N.^(-(order+1)/2),'-.','LineWidth',3,'markersize',14)
hold off
if mod(order,2) == 1
    word1 = ['slope=-', num2str((1+order)/2)];
    word2 = ['slope=-', num2str((order+1)/2)];
else
    word1 = ['slope=-', num2str(order+1), '/2'];
    word2 = ['slope=-', num2str(order+1), '/2'];
end
legend('|u_h-u_h^{dir}|_1','||u_h-u_h^{dir}||_0', word1,word2);
set(gca, 'FontSize', 14)
title('Errors of the finite element method','fontsize',18)
xlabel('Number of elements','fontsize',18)
ylabel('Errors','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);