function E=RIMCv0(A,B,s)
% The code computes generalized eigenvalues of
%         Ax = lambda Bx
% in a rectangle s on the complex plane.
% Input: 
%        A    -- N x N matrix
%        R    -- N x N matrix
%        s    -- 4 x 1 vector [xmin ymin xmax ymax]
% Output: l
%        lambda -- generalized eigenvalues in s
% Jiguang Sun, 05/09/2017, jiguangs@mtu.edu
% Please report bugs to jiguangs@mtu.edu
% Copyright (c) 2017, Jiguang Sun, all rights reserved. 
% THIS SOFTWARE IS PROVIDED "AS IS".
% Redistribution and use in source and binary forms, with or without 
% modification, for academic purpose only are permitted.
% References:
% 1. J.Sun and A.Zhou, Finite Element Methods for Eigenvalue Problems, CRC Press, 2016.
% 2. R.Huang, A.Struthers, J.Sun and R.Zhang,  
%    Recursive integral method for transmission eigenvalues, 
%    Journal of Computational Physics, Vol. 327, 830-840, 2016.
% 3. R Huang, J Sun, C Yang, 
%    Recursive Integral Method with Cayley Transformation
%     - arXiv preprint arXiv:1705.01646, 2017 - arxiv.org

l1=s(1); r1=s(2); l2=s(3); r2=s(4);
tol = 1.0e-12;
size_Krylov_space=50;        
delta_len = min(min((l2-l1)/50,(r2-r1)/50),0.04);


global f_test;
N = length(A);
f_test=rand(N,1)+2*exp(1i*rand(N,1));
Square_point=zeros(10000,5);
Selected_point=zeros(10000,5);
final_point=zeros(10000,1);
tol_robust=0.1*tol;

global p;
global p_test;
p=[0,1/2,1/4,3/4,1/8,5/8,3/8,7/8,1/16,9/16,5/16,13/16,3/16,11/16,7/16,15/16];
p_test=[0,1/2,1/4,3/4,1/8,5/8,3/8,7/8];

f_test=f_test/norm(f_test);
len_l=l2-l1; len_r=r2-r1;
size_length=ceil(len_l/delta_len);
size_width=ceil(len_r/delta_len);

global k;
global number_reduce_sys;
global V_h;
global D_h;
global V_h_inv_v;
global H_para;
global shift_phi;

k=size_Krylov_space;
V_h=zeros(k,k,100);
D_h=zeros(k,100);
V_h_inv_v=zeros(k,100);
H_para=zeros(1,100);
shift_phi=zeros(1,100);
number_reduce_sys=0;

for i=1:size_length
    for j=1:size_width    
        Square_point((j-1)*size_length+i,1)=l1+(i-1)*delta_len;
        Square_point((j-1)*size_length+i,2)=r1+(j-1)*delta_len;
        Square_point((j-1)*size_length+i,3)=l1+i*delta_len;
        Square_point((j-1)*size_length+i,4)=r1+j*delta_len;
    end
end
Square_number=size_length*size_width;

for i=1:Square_number
    ind=[i-1,i-size_length,i-size_length+1,i-size_length-1];
    ind=ind(ind>0);
    if isempty(ind)==1 | number_reduce_sys == 0
       shift=(Square_point(i,1)+Square_point(i,3))/2+1i*(Square_point(i,2)+Square_point(i,4))/2;
       arnoldi(A,B,shift);
       Square_point(i,5)=number_reduce_sys;
    else       
       flag=check_ind(Square_point(ind,5),Square_point(i,1:4),tol_robust); 
       if (flag~=0) 
           Square_point(i,5)=flag;
       else
          shift=(Square_point(i,1)+Square_point(i,3))/2+1i*(Square_point(i,2)+Square_point(i,4))/2; 
          arnoldi(A,B,shift);           
          Square_point(i,5)=number_reduce_sys;
       end
    end
end

tol_iter=ceil(log2(delta_len/tol_robust));

for i=1:tol_iter
    Selected_number=0;    
    for j=1:Square_number
        if (Check_square(Square_point(j,:))==1)
           Selected_number= Selected_number+1;
           Selected_point(Selected_number,:)=Square_point(j,:);
        end                        
    end
    Square_number=Selected_number*4;
    for j=1:Selected_number
        l1=Selected_point(j,1);
        r1=Selected_point(j,2);
        l2=Selected_point(j,3);
        r2=Selected_point(j,4);
        index=Selected_point(j,5);         
        
        Square_point((j-1)*4+1,1:5)=[l1,r1,(l1+l2)/2,(r1+r2)/2,index];
        Square_point((j-1)*4+2,1:5)=[(l1+l2)/2,r1,l2,(r1+r2)/2,index];
        Square_point((j-1)*4+3,1:5)=[(l1+l2)/2,(r1+r2)/2,l2,r2,index];
        Square_point((j-1)*4+4,1:5)=[l1,(r1+r2)/2,(l1+l2)/2,r2,index];                       
    end
end

Num_Eigenvalue=0;
for i=1:Selected_number
    if (Selected_point(i,5)~=0)
       tmp=Selected_point(i,1:4);
       num=1;
       for j=1:Selected_number
           if (norm(Selected_point(i,1:4)-Selected_point(j,1:4))<16*tol)
              num=num+1;
              tmp=tmp+Selected_point(j,1:4);
              Selected_point(j,5)=0;
           end
       end
       Num_Eigenvalue=Num_Eigenvalue+1;
       tmp=tmp./num;
       final_point(Num_Eigenvalue)=(tmp(1)+tmp(3))/2+1i*(tmp(2)+tmp(4))/2;
    end
end

E=final_point(1:Num_Eigenvalue); 

function flag=Check_square(point)
flag=0;
center=(point(1)+point(3))/2+1i*(point(2)+point(4))/2;
radius=sqrt((real(center)-point(1))^2+(imag(center)-point(2))^2);
ind=point(5);

global p;
global k;
global D_h;
global V_h_inv_v;
global shift_phi;

Para_p=center+radius*exp(1i*2*pi*p);
iter=[0,2,4,8,16];
w=[1/2,1/4,1/8,1/16];
res=zeros(1,4);
u=rand(k,1);
level = 3;
for i=1:level 
    for j=iter(i)+1:iter(i+1)
        tmp_v=u'*((V_h_inv_v(:,ind)./(ones(k,1)+(shift_phi(ind)-Para_p(j))*D_h(:,ind))));
        res(i)=res(i)+w(i)*(Para_p(j)-center)*tmp_v;
    end
    res(i+1)=res(i)*1/2;
end

if (abs(norm(res(level-1))/norm(res(level)))<15)
   flag=1; 
end

function ind=check_ind(index,point,tol)

global p_test;
global k;
global V_h;
global D_h;
global V_h_inv_v;
global H_para;
global shift_phi;

center=(point(1)+point(3))/2+1i*(point(2)+point(4))/2;
radius=sqrt((real(center)-point(1))^2+(imag(center)-point(2))^2);
Para_p=center+radius*exp(1i*2*pi*p_test);

ind=0;
for i=1:length(index)
    flag=1;
    for j=1:length(Para_p)
        if abs(H_para(index(i))*V_h(k,:,index(i))*(V_h_inv_v(:,index(i))./((ones(k,1)+(shift_phi(index(i))-Para_p(j))*D_h(:,index(i))))))>tol
            flag=0;
        end    
    end
    if flag==1
       ind=index(i);
       break;
    end
end

function arnoldi(A,B,phi)

global f_test;
global k;
global number_reduce_sys;
global V_h;
global D_h;
global V_h_inv_v;
global H_para;
global shift_phi;

I=sparse(k,1);
I(1)=1;
number_reduce_sys=number_reduce_sys+1;
[L,U,P,Q]=lu(A-phi*B);
b=Q*(U\(L\(P*f_test)));
V1=zeros(length(A),k+1);
H1=zeros(k+2,k+1);
V1(:,1)=b/norm(b);
for i=1:k
    v1=Q*(U\(L\(P*(B*V1(:,i)))));
    for j=1:i
        H1(j,i)=V1(:,j)'*v1;       
        v1=v1-H1(j,i)*V1(:,j);        
    end
    H1(i+1,i)=norm(v1);    
    V1(:,i+1)=v1/H1(i+1,i);            
end    
V1=V1(:,1:k);
[V,D]=eig(H1(1:k,1:k));
V_h(:,:,number_reduce_sys)=V;
D_h(:,number_reduce_sys)=diag(D);
V_h_inv_v(:,number_reduce_sys)=V\I;
shift_phi(number_reduce_sys)=phi;
H_para(number_reduce_sys)=H1(k+1,k);
