function [qpbd,pX,upperQ,fval,cpuqp,cpuqpupper] = myqpbd(L,m)

% This computes qp bound
% The input is Laplacian L of a graph, and the partition vector m.
% input
% L    - Laplacian
% m    - set size vector
% Output
% qpbd           - quadratic lower bounds A
% pX             - Partition matrix attaining upper bound
% upperQ         - upper bound
% fval           - optimal value of the QP solved
% cpuqp          - cpu time for lower bound
% cpuqpupper     - cpu time for upper bound

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

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

en = ones(n,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

% compute scalar product and constant
Bhat = W'*Mtilde*B*Mtilde*W;
Bhat = (Bhat + Bhat')/2;

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

C = V'*A*Xhat*B*Mtilde*W; % a multiple of 2 is not needed for input to quad prog

[Ahatv,eigvav] = eig(Ahat); 
eigvav = diag(eigvav);
if norm(Ahatv'*Ahatv - eye(size(Ahat,1)),'fro') > 1e-6
  fprintf(' error in orthogonality in Ahat; reorthogonalize\n')
  [Ahatv,~] = qr(Ahatv); % orthogonalize
end

[Bhatv,eigbtilde] = eig(Bhat);
eigbtilde = diag(eigbtilde);
if norm(Bhatv'*Bhatv - eye(size(Bhat,1)),'fro') > 1e-6
  fprintf(' error in orthogonality in Bhat; reorthogonalize\n')
  [Bhatv,~] = qr(Bhatv); % orthogonalize
end

Aminscalarprod = minscalarprod(eigbtilde,eigvav);
Aconst = (sum(sum(A)))*(m'*B*m)/(n*n); % alpha in the writeup
qpconst = (Aconst + Aminscalarprod)/2;

% Find the Hessian for QP bound: first solve a dual LP; needs to be
% optimized !!!

diagprod = eigbtilde*eigvav'; % ij-th entry is eigb_i*eiga_j
diagprod = diagprod(:);

% options = optimset('largescale','off','simplex','on');
% options = optimoptions(@linprog,'algorithm','simplex');
Alinprog = [kron(ones(n-1,1),eye(k-1)), kron(eye(n-1),ones(k-1,1))];
% st = linprog(-ones(2*n-2,1),Alinprog,diagprod,[],[],[],[],[],options); % maximize

% use cvx %% write directly into SDPT3/SeDuMi format?

cvx_solver SDPT3

        cvx_begin
          variable st(n+k-2,1)
          minimize -sum(st)
          subject to
            Alinprog*st <= diagprod;
            st(k:end) <= 0;
        cvx_end
        
if abs(Aminscalarprod - sum(st)) > 1e-3 || sum(max( Alinprog*st - diagprod,0)) > 1e-6 || sum(max(st(k:end),0)) > 1e-6 % make sure the LP is solved reasonably accurately
  fprintf(' LP subproblem in QP bound not solved accurately\n')
  keyboard
end

S = Ahatv*diag(st(k:end))*Ahatv';
T = Bhatv*diag(st(1:k-1))*Bhatv';

Qhat = kron(Bhat,Ahat) - kron(eye(k-1),S) - kron(T,eye(n-1)); 
Qhat = (Qhat + Qhat')/2;
% keyboard

% min(eig(Qhat))

Chat = C(:);
Aineq = kron(Mtilde*W,V);

% use cvx %% write directly into SDPT3/SeDuMi format?
% On Dec 21: much smaller QP

        cvx_begin
          variable x((n-1)*(k-1),1)
          minimize 0.5*quad_form(x,Qhat) + dot(x,Chat)
          subject to
            Aineq*x >= -Xhat(:);
        cvx_end

fval = 0.5*x'*Qhat*x + Chat'*x;
qpbd = fval + qpconst;
cpuqp = toc(cpuqpstart);

% compute feasible solution for upper bound computation
Z = reshape(x,n-1,k-1);

X = Xhat + (V*Z)*(W'*Mtilde);
cpuqpupperstart = tic;
pX = proj_partition(X,m);
cpuqpupper = toc(cpuqpupperstart);

upperQ = trace(A*pX*B*pX')/2;
