diff --git a/data/gleich/data/Erdos02-cc.smat.gz b/data/gleich/data/Erdos02-cc.smat.gz
new file mode 100644
index 0000000..2fd84d4
Binary files /dev/null and b/data/gleich/data/Erdos02-cc.smat.gz differ
diff --git a/src/index_estrada.m b/src/index_estrada.m
new file mode 100644
index 0000000..fb0ee65
--- /dev/null
+++ b/src/index_estrada.m
@@ -0,0 +1,42 @@
+% EE = index_estrada(c, ab )
+% EE = index_estrada(Afuns, n, nZ, N, ab)
+% EE = index_estrada(As, nZ, N, ab)
+%
+% Compute the Estrada index tr(exp(A)) by the Chebyshev expansion of 
+% exponential function and moments of A; the spectrum of As should 
+% already lie in [-1,1], and the rescaling parameters should be given
+% as the last argument.
+%
+% Inputs:
+%    c: A column vector of N Chebyshev moment estimates of A
+%    As: Rescaled matrix or function to apply matrix (to multiple RHS)
+%    n: Dimension of the space (if A is a function)
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    EE: Estrada index of the form tr(exp(A))
+%
+function EE = index_estrada(varargin)
+    % If rescaling parameters are given
+    if length(varargin{end})==2
+        ab = varargin{end};
+        varargin = varargin(1:end-1);
+    else
+        ab = [1,0];
+    end
+    
+    % Use the moments if they are given; otherwise compute them
+    if min(size(varargin{1}))==1 && length(varargin{1})>1
+        c = varargin{1};
+    else
+        [c,~] = moments_cheb_dos(varargin{:});
+    end
+    
+    % Get the coefficients of Chebyshev expansion of exponential
+    N = length(c);
+    w = moments_exponential(N,ab(1));
+    EE=exp(ab(2))*w'*c;
+    
+end
\ No newline at end of file
diff --git a/src/index_inform.m b/src/index_inform.m
new file mode 100755
index 0000000..9158f7e
--- /dev/null
+++ b/src/index_inform.m
@@ -0,0 +1,52 @@
+% IC = index_inform(A, nZ, N)
+%
+% Compute the information centrality by the Chebyshev expansion of 
+% 1/(ab(1)*x_ab(2)) and moments of (L+ee^T/n).
+%
+%    A: Adjacency matrix
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%
+% Output:
+%    IC: Information centrality
+%
+function IC = index_inform(A, nZ, N)
+    if nargin<3
+        N = 150;
+    end
+    
+    if nargin<2
+        nZ=100;
+    end
+    
+    n = size(A,1);
+   
+    % Compute the eigenrange of Laplacian
+    L = matrix_laplacian(A);
+    opts = [];
+    opts.isreal = 1;
+    opts.issym = 1;
+    range = sort(eigs(L, 2, 'sm', opts));
+    range(3) = eigs(L, 1, 'lm', opts);
+    range(1) = 1;
+    range = sort(range);
+    range = range([1,end]);
+    
+    % Linear rescaling of (L+ee^T/n)
+    Lfun = @(X) L*X;
+    Bfun = @(X) (bsxfun(@plus,mean(X,1),Lfun(X)));
+    [Bfuns,ab] = rescale_mfunc(Bfun,n,range);
+    
+    % Compute the ldos moments
+    [c,~] = moments_cheb_ldos(Bfuns,n,nZ,N);
+    
+    % Compute the Chebyshev expansion of 1/(ab(1)x+ab(2))
+    z = -ab(2)/ab(1);
+    w = 1./(sqrt(z^2-1)*(z-sqrt(z^2-1)).^(0:N-1))/ab(1);
+    w(2:N)=2*w(2:N);
+    
+    % Compute information centrality
+    d = c'*w'-(n-1)/n^2;
+    IC = 1./(d+mean(d)-2/n^2);
+    
+end
\ No newline at end of file
diff --git a/src/index_sub_exp.m b/src/index_sub_exp.m
new file mode 100644
index 0000000..2522264
--- /dev/null
+++ b/src/index_sub_exp.m
@@ -0,0 +1,54 @@
+% ESC = index_sub_exp(c, ab)
+% ESC = index_sub_exp(beta, c, ab)
+% ESC = index_sub_exp(beta, Afuns, n, nZ, N, ab)
+% ESC = index_sub_exp(beta, As, nZ, N, ab)
+%
+% Compute the exponential subgraph centrality diag(exp(beta*A)) by the
+% Chebyshev expansion of exponential function and moments of A; the 
+% spectrum of As should already lie in [-1,1], and the rescaling 
+% parameters should be given as the last argument.
+%
+% Inputs:
+%    beta: Exponential coefficient
+%    c: An N-by-n matrix of ldos moment estimates of A
+%    As: Rescaled matrix or function to apply matrix (to multiple RHS)
+%    n: Dimension of the space (if A is a function)
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    ESC: Exponential subgraph centraility of the form diag(exp(beta*A))
+%
+function ESC = index_sub_exp(varargin)
+    % If rescaling parameters are given
+    if length(varargin{end})==2
+        ab = varargin{end};
+        varargin = varargin(1:end-1);
+    else
+        ab = [1,0];
+    end
+    
+    % Check if an exponential coefficient is given (default:1)
+    if isnumeric(varargin{1}) && isscalar(varargin{1})
+        beta = varargin{1};
+        varargin = varargin(2:end);
+    else
+        beta = 1;
+    end
+    
+    ab = beta*ab;
+    
+    % Check if ldos moments are given; otherwise calculate an estimate
+    if size(varargin{1},1)~=size(varargin{1},2)
+        c = varargin{1};
+    else
+        [c,~] = moments_cheb_ldos(varargin{:});
+    end
+    
+    % Get the Chebyshev expansi0n of exp(beta*x)
+    N = size(c,1);
+    w = moments_exponential(N,ab(1));
+    ESC = exp(ab(2))*(c'*w);
+    
+end
\ No newline at end of file
diff --git a/src/index_sub_res.m b/src/index_sub_res.m
new file mode 100644
index 0000000..9d1b9a5
--- /dev/null
+++ b/src/index_sub_res.m
@@ -0,0 +1,59 @@
+% ESC = index_sub_res(c, ab)
+% ESC = index_sub_res(alpha, c, ab)
+% ESC = index_sub_res(alpha, Afuns, n, nZ, N, ab)
+% ESC = index_sub_res(alpha, As, nZ, N, ab)
+%
+% Compute the resolvent subgraph centrality diag((I-alpha*A)^-1) by 
+% the Chebyshev expansion of exponential function and moments of A; 
+% the spectrum of As should already lie in [-1,1], and the rescaling 
+% parameters should be given as the last argument. 
+%
+% Inputs:
+%    alpha: Resolvent coefficient
+%    c: An N-by-n matrix of ldos moment estimates of A
+%    As: Rescaled matrix or function to apply matrix (to multiple RHS)
+%    n: Dimension of the space (if A is a function)
+%    nZ: Number of probe vectors with which we want to compute moments
+%    N: Number of moments to compute
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    RSC: Resolvent subgraph centraility of the form diag((I-alpha*A)^-1)
+%
+function RSC = index_sub_res(varargin)
+    % If rescaling parameters are given
+    if length(varargin{end})==2
+        ab = varargin{end};
+        varargin = varargin(1:end-1);
+    else
+        ab = [1,0];
+    end
+    
+    % Check if an exponential coefficient is given (default:0.9)
+    if isnumeric(varargin{1}) && isscalar(varargin{1})
+        alpha = varargin{1};
+        varargin = varargin(2:end);
+    else
+        alpha = 0.9;
+    end
+    
+    % Throw a warning if choice of alpha does not guarantee convergence
+    if alpha>=1/sum(ab)
+        warning(['Chebyshev expansion will not converge because'...
+                ' function has pole in [-1,1]. alpha must be less'...
+                ' than 1/sum(ab)']);
+    end
+    
+    % Check if ldos moments are given; otherwise calculate an estimate
+    if size(varargin{1},1)~=size(varargin{1},2)
+        c = varargin{1};
+    else
+        [c,~] = moments_cheb_ldos(varargin{:});
+    end
+    
+    % Get the Chebyshev expansion of 1/(1-alpha*x)
+    N = size(c,1);
+    w = moments_resolvent(N,alpha,ab);
+    RSC = c'*w;
+    
+end
\ No newline at end of file
diff --git a/src/moments_exponential.m b/src/moments_exponential.m
new file mode 100644
index 0000000..fd6e20c
--- /dev/null
+++ b/src/moments_exponential.m
@@ -0,0 +1,22 @@
+% [c] = moments_exponential(N)
+% [c] = moments_exponential(N, beta)
+%
+% Compute Chebyshev moments 0 through N-1 of exponential function
+% f(x)=exp(beta*x) on [-1,1]
+%
+% Input:
+%    N: Number of moments
+%    beta: Exponential coefficient
+%
+% Output:
+%    c: A column vector of N moments
+%
+function c = moments_exponential(N,beta)
+    if nargin<2
+        beta = 1;
+    end
+    
+    c = besseli(0:N-1,beta);
+    c(2:N) = c(2:N)*2;
+    c=c';
+end
\ No newline at end of file
diff --git a/src/moments_resolvent.m b/src/moments_resolvent.m
new file mode 100644
index 0000000..e855929
--- /dev/null
+++ b/src/moments_resolvent.m
@@ -0,0 +1,30 @@
+% [c] = moments_resolvent(N)
+% [c] = moments_resolvent(N, alpha)
+%
+% Compute Chebyshev moments 0 through N-1 of the function
+% f(x)=1/(1-alpha*x) on [-1,1]
+%
+% Input:
+%    N: Number of moments
+%    alpha: Resolvent coefficient in (0,1)
+%    ab: Rescaling parameters. Original matrix should be ab(1)*A+ab(2)
+%
+% Output:
+%    c: A column vector of N moments
+%
+function c = moments_resolvent(N,alpha,ab)
+    if nargin<3
+        ab = [1,0];
+    end
+    if nargin<2
+        alpha = 0.9;
+    end
+    
+    ab = alpha*ab;
+    
+    z = (1-ab(2))/ab(1);
+    c = 1./(sqrt(z^2-1)*(z+sqrt(z^2-1)).^(0:N-1))/ab(1);
+    c(2:N) = c(2:N)*2;
+    c=c';
+    
+end
\ No newline at end of file
diff --git a/test/test_index.m b/test/test_index.m
new file mode 100644
index 0000000..68985e1
--- /dev/null
+++ b/test/test_index.m
@@ -0,0 +1,255 @@
+clear
+clc
+
+addpath([fileparts(pwd),'/src']);
+addpath([fileparts(pwd),'/data/gleich']);
+
+pass = 'Test passed.';
+fail = 'Test failed.';
+lb = '------------------------------------------------------------\n';
+
+% test moments_exponential, which computes Chebyshev coefficient of
+% exp(beta*x) bewteen [-1,1]
+
+beta = 0.9;
+N = 20;
+fprintf(['Test momemts_expoential with beta = ',num2str(beta),...
+        ' and N = ',num2str(N),'.\n']);
+w = moments_exponential(N,beta);
+xx = linspace(-1+1e-8,1-1e-8,1001);
+yy1 = plot_chebp(w,xx);
+yy2 = exp(beta*xx);
+relerr = norm(yy1-yy2)/norm(yy2);
+if relerr < 1e-12
+    fprintf(['The relative error is ',num2str(relerr),'. ',pass,'\n']);
+else
+    fprintf(['The relative error is ',num2str(relerr),'. ',fail,'\n']);
+end
+fprintf(lb);
+
+% test moments_resolvent, which computes Chebyshev coefficient of 
+% 1/(1-alpha*x) on [-1,1].
+
+alpha = 0.9;
+N = 70;
+fprintf(['Test momemts_resolvent with alpha = ',num2str(alpha),...
+        ' and N = ',num2str(N),'.\n']);
+w = moments_resolvent(N,alpha);
+xx = linspace(-1+1e-8,1-1e-8,1001);
+yy1 = plot_chebp(w,xx);
+yy2 = 1./(1-alpha*xx);
+relerr = norm(yy1-yy2)/norm(yy2);
+if relerr < 1e-12
+    fprintf(['The relative error is ',num2str(relerr),'. ',pass,'\n']);
+else
+    fprintf(['The relative error is ',num2str(relerr),'. ',fail,'\n']);
+end
+fprintf(lb);
+
+% Computing eigenvalue of Erdos02-cc for test purposes
+fprintf(['Loading Erdos02-cc and computing its eigen-decomposition \n'...
+        'for testing purposes. Takes 50s.\n']);
+
+A = load_graph('Erdos02-cc');
+L = matrix_normalize(A);
+[VA,eA] = eig(full(A),'vector');
+[VL,eL] = eig(full(L),'vector');
+fprintf(lb)
+
+% Test on index_estrada, which computes the Estrada index of matrix
+fprintf('Test index_estrada on normalized adjacency. No rescaling.\n');
+EE2 = sum(exp(eL));
+N = 10;
+for nZ =10:10:50
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    EE1 = index_estrada(L,nZ,N);
+    relerr = abs(EE1-EE2)/EE2;
+    if relerr < 0.01
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+fprintf(['Test index_estrada on normalized adjacency with '... 
+    'precomputed dos. No rescaling.\n']);
+N = 10;
+for nZ =10:10:50
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    c = moments_cheb_dos(L,nZ,N);
+    EE1 = index_estrada(c);
+    relerr = abs(EE1-EE2)/EE2;
+    if relerr < 0.01
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+fprintf('Test index_estrada on adjacency. With rescaling.\n');
+EE2 = sum(exp(eA));
+N = 30;
+[As,ab] = rescale_matrix(A);
+for nZ =200:200:1000
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    EE1 = index_estrada(As,nZ,N,ab);
+    relerr = abs(EE1-EE2)/EE2;
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+% Test on index_sub_exp, which computes the exponential subgraph
+% centrality
+
+fprintf(['Test index_sub_exp on normalized adjacency. No rescaling.'...
+        'beta = 0.9\n']);
+beta = 0.9;
+ESC2 = VL.^2*exp(beta*eL);
+N = 10;
+for nZ =200:200:1000
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    ESC1 = index_sub_exp(beta,L,nZ,N);
+    relerr = norm(ESC1-ESC2)/norm(ESC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+fprintf(['Test index_sub_exp on normalized adjacency with '... 
+    'precomputed ldos. No rescaling. beta = 0.9\n']);
+beta = 0.9;
+N = 10;
+for nZ =200:200:1000
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    c = moments_cheb_ldos(L,nZ,N);
+    ESC1 = index_sub_exp(beta,c);
+    relerr = norm(ESC1-ESC2)/norm(ESC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+fprintf('Test index_sub_exp on adjacency. With rescaling. beta = 0.9\n');
+beta = 0.9;
+ESC2 = VA.^2*exp(beta*eA);
+N = 20;
+[As,ab] = rescale_matrix(A);
+for nZ =400:400:2000
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    ESC1 = index_sub_exp(beta,As,nZ,N,ab);
+    relerr = norm(ESC1-ESC2)/norm(ESC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+% Test on index_sub_res, which computes the resolvent subgraph
+% centrality
+
+fprintf(['Test index_sub_res on normalized adjacency. No rescaling.'...
+        ' alpha = 0.9\n']);
+alpha = 0.9;
+RSC2 = diag((eye(5534)-alpha*L)^(-1));
+N = 10;
+for nZ =300:300:1500
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    RSC1 = index_sub_res(alpha,L,nZ,N);
+    relerr = norm(RSC1-RSC2)/norm(RSC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+fprintf(['Test index_sub_res on normalized adjacency with '... 
+    'precomputed ldos. No rescaling. alpha = 0.9\n']);
+alpha = 0.9;
+N = 10;
+for nZ =300:300:1500
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    c = moments_cheb_ldos(L,nZ,N);
+    RSC1 = index_sub_res(alpha,c);
+    relerr = norm(RSC1-RSC2)/norm(RSC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+[As,ab] = rescale_matrix(A);
+alpha = 0.9/sum(ab);
+fprintf(['Test index_sub_res on adjacency. With rescaling. alpha = ',...
+    num2str(alpha),'\n']);
+RSC2 = diag((eye(5534)-alpha*A)^(-1));
+N = 20;
+for nZ =300:300:1500
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    RSC1 = index_sub_res(alpha,As,nZ,N,ab);
+    relerr = norm(RSC1-RSC2)/norm(RSC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+
+% Test on index_inform, which computes the information centrality
+
+fprintf('Test index_inform on adjacency.\n');
+B = matrix_laplacian(A)+1;
+K = B^(-1);
+I = repmat(diag(K),[1,5534])+repmat(diag(K)',[5534,1])-2*K; 
+IC2 = 1./(mean(I,2));
+N = 150;
+for nZ =100:100:500
+    fprintf(['With nZ = ',num2str(nZ),' and N = ',num2str(N),':\n']);
+    IC1 = index_inform(A,nZ,N);
+    relerr = norm(IC1-IC2)/norm(IC2);
+    if relerr < 0.05
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                pass,'\n']);
+    else
+        fprintf(['The relative error is ',num2str(relerr),'. '...
+                fail,'\n']);
+    end
+end
+fprintf(lb)
+fprintf('Test Over\n')
\ No newline at end of file