function [optvalt3A,pX,upperSDP,cpusdp,cpuupper,cpusdp0] = mysdpbd(L,m)

% This computes sdp bound
% The input is Laplacian L of a graph, and the partition vector m.
% input
% L    - Laplacian
% m    - set size vector
% Output
% optval3A     - SDP lower bound (two values from output of SDPT3)
% pX           - Partition matrix attaining upper bound
% upperSDP     - SDP upper bound
% cpusdp       - cpu time for lower bound (including time for constructing input to SDPT3 and solving the problem by SDPT3)
% cpuupper     - cpu time for upper bound
% cpusdp0      - cpu time for constructing into to SDPT3

cpusdpstart = tic;
% Initialization
m = m(:);
n = sum(m);
k = length(m);

en = ones(n,1);
en1 = ones(n-1,1);
ek1 = ones(k-1,1);

Vk1 = [eye(k-1); -ek1'];
Vn1 = [eye(n-1); -en1'];
Vhat = [1 zeros(1,(n-1)*(k-1)); kron(m/n,en)  kron(Vk1,Vn1)];
%%%????calculate size of Jbar first!! will save time!!!???
% record positions for the gangster operator
Jbar = [];         %   top left ignored for  now
for q = 1:n,
	for r = 1:k,
		for p = 1:r-1,
			i = (p-1)*n + q;
			j = (r-1)*n + q;
			Jbar = [Jbar; i j ];
		end
	end
end
Jbar = Jbar + 1;    % correct for 0 row

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

LA = [0 zeros(1,k*n); zeros(k*n,1)  kron(B,A)]; 
VLAV = Vhat'*LA*Vhat;

% solve for sdp relaxation
nk1 = (n-1)*(k-1) + 1;
tnk1 = nk1*(nk1 + 1)/2;

% construct input to SDPT3
blk{1,1} = 's'; blk{1,2} = nk1;

%profile on

J0 = Jbar';
At3 = sparse(size(J0,2),tnk1);

for kk = 1:size(J0,2)
  ii = J0(1,kk);
  jj = J0(2,kk);
  vii = (Vhat(ii,:))';
  vjj = (Vhat(jj,:))';
  At3(kk,:) = svec(blk,vii*vjj'+vjj*vii');
  %keyboard
end

Atfull{1} = [sparse(1,1,1,tnk1,1,1) At3'];   % one block, so tnk1 by #constr. 
bt3 = sparse(1,1,1,size(J0,2)+1,1,1);
Ct3A = VLAV;
cpusdp0 = toc(cpusdpstart); % time for writing input to SDPT3

OPTIONS.gaptol = 1e-6;
[optvalt3A,Xt3,~,~,~] = sqlp(blk,Atfull,Ct3A,bt3,OPTIONS); 
cpusdp = toc(cpusdpstart); % total time

Yoptt3A = Vhat*Xt3{1}*Vhat';
Yoptt3A = (Yoptt3A + Yoptt3A')/2;

optvalt3A = optvalt3A/2; 

cpuupperstart = tic;
Xoptsdpt3A = reshape(Yoptt3A(2:end,1),n,k); % use first column of Y
pX1 = proj_partition(Xoptsdpt3A,m); % projection to find feasible solution for upper bound

[Xsdpcolv,~] = eig(Yoptt3A); % use eigenvector
Xsdpcol = reshape(Xsdpcolv(2:end,end),n,k);
temp = sum(Xsdpcol,2);
Xsdpcol = Xsdpcol/temp(1);   % normalize column sums
if norm(Xsdpcol*ones(k,1)-ones(n,1)) +norm(Xsdpcol'*ones(n,1)-m) > .00001,
  disp('error: Xsdpcol not EE feasible?')
  keyboard
end
pX2 = proj_partition(Xsdpcol,m); % projection to find feasible solution for upper bound
cpuupper = toc(cpuupperstart);

upperSDP1 = trace(A*pX1*B*pX1')/2;
upperSDP2 = trace(A*pX2*B*pX2')/2;
if upperSDP1 > upperSDP2
  upperSDP = upperSDP2; pX = pX2;
else upperSDP = upperSDP1; pX = pX1;
end
