function [Lbnd,Abndeq,pXL,pXA,upperL,upperA,cpueigbd,cpuupper] = myprojeigbd(L,m)

% This computes projected eigenvalue bound, both -L and A.
% The input is Laplacian L of a graph, and the partition vector m.
% input
% L    - Laplacian
% m    - set size vector
% Output
% Lbnd, Abndeq   - projected eigenvalue lower bounds using -L and A
% pXL, pXA       - Partition matrix attaining upper bounds for -L and A
% upperL, upperA - projected eigenvalue upper bounds using -L and A
% cpueigbd       - cpu time for lower bounds
% cpuupper       - cpu time for upper bounds

cputimestarteig = tic;   % start time for eigbnd
% Initialization
m = m(:);
n = sum(m);
k = length(m);

mtilde = sqrt(m);
Mtilde = diag(mtilde);

en = ones(n,1);
ek = ones(k,1);

[q,~] = qr(en);
V = q(:,2:end);
[q,~] = qr(mtilde);
W = q(:,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%
Xhat = (1/n)*en*m';
J = ones(k-1);
B = [ J-eye(k-1)   zeros(k-1,1)
    zeros(1,k-1)         0        ] ;
A = -L + diag(diag(L)); % compute adjacency matrix

%% solve for linear part using simplex method
v1 = sum(A,2);
v2 = zeros(n,1);
count = 0;
for i = 1:k-1
  v2(count+1:count+m(i)) = n - m(end) - m(i);
  count = count + m(i);
end
fvalsimplexeq = minscalarprod(v1,v2)/n;

% compute the rest: scalar product and constant
Bhat = W'*Mtilde*B*Mtilde*W;
Bhat = (Bhat + Bhat')/2;
mL = -L; %negative Laplacian

Ahat = V'*A*V; 
Lhat = V'*(mL)*V; 
Ahat = (Ahat + Ahat')/2;
Lhat = (Lhat + Lhat')/2;

[Ahatv,eigvav] = eig(Ahat); 
[Lhatv,eigvlv] = eig(Lhat);
[eigvav,indseiga] = sort(diag(eigvav),'descend');
[eigvlv,indseigl] = sort(diag(eigvlv),'descend');

Ahatv = Ahatv(:,indseiga);    % sorted eigenvectors
Lhatv = Lhatv(:,indseigl);    % sorted eigenvectors

[Bhatv,eigbtilde] = eig(Bhat);  
[eigbtilde,indsbtilde] = sort(diag(eigbtilde),'ascend');
Bhatv = Bhatv(:,indsbtilde);

Aminscalarprod = minscalarprod(eigbtilde,eigvav);
Aconst = (sum(sum(A)))*(m'*B*m)/(n*n);
Abndeq = (Aminscalarprod + 2*fvalsimplexeq - Aconst)/2;
Lbnd = (minscalarprod(eigbtilde,eigvlv))/2;

cpueigbd = toc(cputimestarteig);

%%%Now find closest feasible solutions for upper bounds
%%% Note first attainment for A,L bounds:
%%% choose the largest k-2 eigs in A,-L for neg eigs in Bhat
%%% choose the smallest    eig  in A,-L for pos eig in Bhat
%%% check orthogonality
if norm(Bhatv'*Bhatv - eye(k-1),'fro') > 1e-6
  disp('error in orthogonality for eigenvectors in Bhatv: reorthogonalize')
  [Bhatv,~] = qr(Bhatv);
end
if norm(Lhatv(:,[1:k-2 n-1])'*Lhatv(:,[1:k-2 n-1]) - eye(k-1),'fro') > 1e-6
  disp('error in orthogonality for eigenvectors in Lhatv: reorthogonalize')
  [Lhatv,~] = qr(Lhatv);
end
if norm(Ahatv(:,[1:k-2 n-1])'*Ahatv(:,[1:k-2 n-1]) - eye(k-1),'fro') > 1e-6
  disp('error in orthogonality for eigenvectors in Ahatv: reorthogonalize')
  [Ahatv,~] = qr(Ahatv);
end

ZoptpqA = Ahatv(:,[1:k-2 n-1])*Bhatv';
testAbnd = (1/2)*trace(Ahat*ZoptpqA*Bhat*ZoptpqA');
if abs(testAbnd - (1/2)*Aminscalarprod) > .0001,
	disp('error in testing Abndeq attainment')
        disp([num2script([(1/2)*Aminscalarprod (1/2)*trace(Ahat*ZoptpqA*Bhat*ZoptpqA')])]);
end
ZoptpqL = Lhatv(:,[1:k-2 n-1])*Bhatv';
testLbnd = (1/2)*trace(Lhat*ZoptpqL*Bhat*ZoptpqL');
if abs(testLbnd - Lbnd) > .0001,
	disp('error in testing Lbnd attainment')
        disp([num2str([Lbnd (1/2)*trace(Lhat*ZoptpqL*Bhat*ZoptpqL')])]);
end

%%%% We now have two approx solutions
XA = Xhat + (V*ZoptpqA)*(W'*Mtilde);
XL = Xhat + (V*ZoptpqL)*(W'*Mtilde);

% project onto partition matrices
cpuupperstart = tic;
pXL = proj_partition(XL,m);
pXA = proj_partition(XA,m);
cpuupper = toc(cpuupperstart);

upperL = trace(A*pXL*B*pXL')/2;
upperA = trace(A*pXA*B*pXA')/2;