function [optvaldnng,pX,upperDNN,cpuDNN,cpuupper,cpuDNN0] = myDNNbd(L,m)

% This computes DNN bound
% The input is Laplacian L of a graph, and the partition vector m.
% input
% L    - Laplacian
% m    - set size vector
% Output
% optvaldnng   - DNN lower bound (two values from output of SDPT3)
% pX           - Partition matrix attaining upper bound
% upperDNN     - DNN upper bound
% cpuDNN       - cpu time for lower bound (including time for constructing input to SDPT3 and solving the problem by SDPT3)
% cpuupper     - cpu time for upper bound
% cpuDNN0      - 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)];

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
nk  = n*k + 1;
nk1 = (n-1)*(k-1) + 1;
tnk = 0.5*nk*(nk + 1);
tnk1 = 0.5*nk1*(nk1 + 1);

blk{1,1} = 's'; blk{1,2} = nk1;

gangster = [zeros(1,nk); zeros(n*k,1) kron(ones(k)-eye(k), eye(n))];
blkgangster{1,1} = 's'; blkgangster{1,2} = n*k+1;
indgangster = svec(blkgangster, gangster)/2^0.5;
indnonneg = zeros(tnk,1);

coeff = sparse(tnk, tnk1);
counter = 0;
for jj = 1:nk,
  for ii = 1:jj,
    counter = counter + 1;
    if ii>=2 && jj>ii,
      indnonneg(counter) = 1;
      vii = (Vhat(ii,:))';
      vjj = (Vhat(jj,:))';
      coeff(counter,:) = svec(blk,vii*vjj' + vjj*vii');
    end
  end
end
% gangster operator
Agangster = coeff(logical(indgangster == 1),:);
% nonnegativity constraints
indnonneg = sparse(indnonneg);
indnonneg = indnonneg - indgangster;
Anonneg = coeff(logical(indnonneg == 1),:);
blk{2,1} = 'l'; blk{2,2} = sum(indnonneg);

Atfull{1} = [sparse(1,1,1,tnk1,1,1) Agangster' Anonneg'];
Atfull{2} = [sparse(sum(indnonneg), 1+size(Agangster,1)), ...
    -speye(sum(indnonneg))];
bt3 = sparse(1,1,1, (sum(indgangster)+1) + sum(indnonneg), 1, 1);
Ct3A = cell(1,2);
Ct3A{1} = VLAV;    
Ct3A{2} = sparse(sum(indnonneg),1);
cpuDNN0 = toc(cpusdpstart); % time to construct input to SDPT3

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

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

optvaldnng = optvaldnng/2; 

cpuupperstart = tic;
Xoptsdpt3A = reshape(Yoptt3A(2:end,1),n,k); % use first column
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);

upperDNN1 = trace(A*pX1*B*pX1')/2;
upperDNN2 = trace(A*pX2*B*pX2')/2;

if upperDNN1 > upperDNN2
  upperDNN = upperDNN2; pX = pX2;
else upperDNN = upperDNN1; pX = pX1;
end
