clear all
delete result.dat
diary result.dat
diary on

fprintf('Choose example \n');
fprintf('1. Two dimensional Laplacian equation \n');
fprintf('2. Three dimensional pyramid quantum dot \n');
fprintf('3. Simple cubic photonic crystal \n');
example = input('==> ');

switch example
    case 1
        n       = input('Mesh number = ');

        mtx_A   = mtx_generate_2DLaplacian(n);
        rhs_f   = Construct_rhs(n);
        
        tic;
        sol     = mtx_A \ rhs_f;
        toc_Lap = toc;
        
        fprintf('CPUTime of solving Laplacian = %14.4e; matrix size = %8.0f \n', toc_Lap, n^2);
    case 2
        [mtx_A,rows,cols,entries,rep,field,symm] = mmread('mtx_A0_SEP_pyramid_77315.mat');
        n              = size(mtx_A,1);
        
        s              = RandStream.create('mt19937ar','seed',5489);
        RandStream.setDefaultStream(s);

        rhs_vec        = randn(n,1);
        fprintf('Solving linear system \n');
        
        tic;
        sol_q          = mtx_A \ rhs_vec; 
        toc_qdot       = toc;
        
        fprintf('CPUTime of solving pyramid qdot = %14.4e; matrix size = %8.0f \n', toc_qdot, n);
        
    case 3
        n_1            = input('Mesh number = ');
        vec_k          = [ 0.1; 0.00; 0.00 ];
        [mtx_A, mtx_B] = mtx_GEP_Maxwell(n_1, vec_k);
        sigma          = 0.1;
        mtx_A_new      = mtx_A - sigma * mtx_B;
        n              = size(mtx_A,1);
        
        s              = RandStream.create('mt19937ar','seed',5489);
        RandStream.setDefaultStream(s);

        rhs_vec        = randn(n,1);
        tic;
        sol_m          = mtx_A_new \ rhs_vec; 
        toc_Max         = toc;
        
        fprintf('CPUTime of solving Maxwell = %14.4e; matrix size = %8.0f \n', toc_Max, n);

end

diary off