function [Aprojeigbnd,upperAbnd,n,Ecard,cputimeprojeigbnd,cputimeupper]=...
                  newV0sparsetest(imax,k,densn)
%%% Generate random problem - adjacency matrix A
%%% Find lower projected eigenvalue bound 
%%%% Using attainment from lower eig part, find upper bound
%%% Uses V0 in the writeup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%INPUT:   k     number of sets
%%%       imax    maximum size of each set
%%%%                    sizes of k sets are randomly generated
%%%       densn   density=densn*n^2  for # nonzeros in random graph
%%%OUTPUT: Aprojeigbnd    lower projected eigenvalue bound using A
%          upperAbnd      upper bound
%          n              number of nodes in graph
%          Ecard          number of edges in graph
%          cputimeprojeigbnd   cpu time for computing lower bound
%          cputimeupper        cpu time for computing upper bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Two nested functions at bottom:
%%  function y=Ahateval(x,A,Vlogical,Vreal,dscale,numl)  % for eigs
%%  function alpha=minscalarprodinnew(u,v)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initialization
disp(' ')
disp('start initializ. for new problem')
cputimeprobgen=tic;   % start time for eigbnd
m=randi(imax,k,1);    % this might yield sets with 'same size'
if sum(m==1)>0,
     m=m+1;
end
n=sum(m);       % number of nodes
A=sprandsym(n,densn);
A(1:n+1:end)=0;
A(abs(A)>0)=1;   % adjacency matrix 0,1 so can be logical
if sum(sum(A)==0)>0,              % zero rows exist - can be deleted???
  disp('Warning!! generated A has zero rows - adding nonzeros')
  inds=find(sum(A)==0);
  for ii=inds,
    temp=sprand(n,1,max(densn,(3/n)));   % at least 3 new els
    temp(abs(temp)>0)=1;
    A(ii,:)=temp';
    A(:,ii)=temp;
    A(ii,ii)=0;
  end
end
A=logical(A);
Ecard=nnz(A)/2;
%%%%%%%%%%%%%%
disp(['cputime for problem generation ',num2str(toc(cputimeprobgen))])
%%
disp(['new problem k,n,density ',num2str([k n (Ecard/(n*(n-1)/2))])]);
cputimestarteig=tic;   % start time for eigbnd
mtilde=sqrt(m);
Mtilde=diag(mtilde);
%%en=ones(n,1);
%%%%%%%%%%%%%%%?????need sparse V!!!!!??????????
%%[qen,~]=qr(en);
%%V=qen(:,2:end);
V=triu(ones(n,n-1));
t1n1=(1:n-1);
V(2:n+1:end)=-t1n1;
%V=sparse(V);
dscale=sqrt(diag(V'*V));
%%%%%%  proper V for  V'AV is ....   Vs=V*inv(diag(dscale))
[qmtilde,~]=qr(mtilde);
W=qmtilde(:,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Xhat=(1/n)*en*m';    % feasible point in EE (linear constraints)
J=ones(k-1);
B=[ J-diag(diag(J))   zeros(k-1,1)
    zeros(1,k-1)         0        ] ;
d=sum(A,2);   % A*Xhat=(1/n)*sum(A,2)*m'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Btilde=Mtilde*B*Mtilde;
%Bhat=W'*Mtilde*B*Mtilde*W;
Bhat=W'*Btilde*W;
Bhat=(Bhat+Bhat')/2;
%%%%%%%%%%
opts.issym=1;
opts.isreal=1;
[Ahatv,eigvav,flagA]=eigs(@(x)Ahateval(x,A,V,dscale),n-1,k-2,'la',opts); 
if flagA~=0,
	disp('ERROR!! not all of largest k-2 eigs converged for Ahat')
end
[eigvav,indseiga]=sort(diag(eigvav),'descend');  %need k-2 largest and smallest
Ahatv=Ahatv(:,indseiga);    % sorted eigenvectors
[Ahatvs,eigvavs,flagAs]=eigs(@(x)Ahateval(x,A,V,dscale),n-1,1,'sa',opts); 
if flagAs~=0,
	disp('ERROR!! smallest eig did not converge for Ahat')
end
eigvav=[eigvav; eigvavs];   % all k-1 sorted eigenvalues
Ahatv=[Ahatv Ahatvs];  % all k-1 sorted eigenvectors
[Bhatv,eigbhat]=eig(Bhat);  
[eigbhat,indsbhat]=sort(diag(eigbhat),'ascend');  %need k-2 smallest
Bhatv=Bhatv(:,indsbhat);
if norm(Bhatv'*Bhatv - eye(size(Bhat,1)),'fro') > 1e-6
  fprintf(' nonorthogonality of eigenvectors of B. Reorthogonalize\n')
  [Bhatv,~] = qr(Bhatv); % reorthogonalize
end
%%%%%Aconst=-trace(A*Xhat*B*Xhat');    % -alpha in Theorem
%Aconstbndpart=(-1/(2*n^2)) * (sum(sum(A))) * (m'*B*m);
Aconstbndpart=(-1/(2*n^2)) * (sum(d)) * (m'*B*m);
Aeigbndpart=eigvav'*eigbhat/2;
%%%% find linear part
v0=(n-m(k))*ones(n-m(k),1);
iii=0;
for ii=1:k-1,
	v0(iii+1:iii+m(ii))= v0(iii+1:iii+m(ii))-m(ii)*ones(m(ii),1);
	iii=iii+m(ii);
end
Alinpart= (1/n)*minscalarprodinnew(d,v0);
Aprojeigbnd= Aconstbndpart+ Aeigbndpart+Alinpart;
cputimeprojeigbnd=toc(cputimestarteig);
disp(['cputime for projeigbnd is ',num2str(cputimeprojeigbnd)]);
disp(' ')
disp(['projected eig lower bound is ',num2str(Aprojeigbnd)])
%if Aprojeigbnd<0,
%	disp('proj eig bnd < 0')
%	spy(A)
%	keyboard
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%Now find closest feasible solutions for upper bounds using eigpart
ZoptpqA=Ahatv*Bhatv';
Xhat=(1/n)*ones(n,1)*m';    % feasible point in EE (linear constraints)
XA=Xhat+(V*ZoptpqA)*(W'*Mtilde);   % attainment of projeig A bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% solve for nearest partition matrix
%options=optimset('largescale','off','simplex','on');
%options=optimset('simplex','on');
Asimplex=[  kron(ones(1,k),speye(n));  kron(speye(k),ones(1,n)) ];
bsimplex=[ones(n,1);m];
Asimplex=Asimplex(1:end-1,:);   % need full row rank
bsimplex=bsimplex(1:end-1,:);   % need full row rank
%% fsimplex from nearest matrix problem min ||XA-X||^2_F
%% equiv to min trace (-XA) X, X partition matrix
fsimplex=-XA(:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[xaoptsimplex,~,flagsimplex,output]=...
     cplexlp(fsimplex,[],[],Asimplex,bsimplex,zeros(n*k,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%options=optimoptions(@linprog,'algorithm','simplex');  % new way
%[xaoptsimplex,~,flagsimplex]=...
%      linprog(fsimplex,[],[],Asimplex,bsimplex,zeros(n*k,1),[],[],options);
cputimeupper=output.time;%cputime-cputimestartup;   % total cpu time for upper LP bnd
if flagsimplex~=1,
	disp('WARNING!! simplex for upper bound did NOT converge')
end
disp(['simplex cputime for upper bnd is ',num2str(cputimeupper)]);
XAoptsimplex=reshape(xaoptsimplex,n,k);
%%%rounding can cause errors close to 0,1
if sum(sum(XAoptsimplex>.0001 & XAoptsimplex<.999)) == 0,
          XAoptsimplex(XAoptsimplex<=.0001)=0;           % round small nos to 0
          XAoptsimplex(XAoptsimplex>=.999)=1;           % round large nos to 1
end
%%%double check for partition matrix
if norm(XAoptsimplex-logical(XAoptsimplex))+...
       norm(sum(XAoptsimplex)-m')+norm(sum(XAoptsimplex,2)-ones(n,1)) > .000001,
   disp('warning!!! - XAoptsimplex is not a partition matrix')
   disp('    return        if you want to continue')
   keyboard
end

upperAbnd=trace(A*XAoptsimplex*B*XAoptsimplex')/2;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=Ahateval(x,A,V,dscale)
%%%%%%  proper V for  V'AV is ....   =V*inv(diag(dscale))
%%      y=(V'*(A*(V*x))); 
      %x=x./dscale; y=(V'*(A*(V*x))); y=y./dscale;
      y=(V'*(A*   (V*   (x./dscale)   )  )  ./dscale); 

end   % of main testsp.... function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function alpha=minscalarprodinnew(u,v)
%function alpha=minscalarprodinnew(u,v)
%%%ensure both are column vectors
u=reshape(u,length(u),1);
v=reshape(v,length(v),1);
t=abs(length(u)-length(v));
if length(u)>=length(v),
	v=[v; zeros(t,1)];
else
	u=[u; zeros(t,1)];
end
u=sort(u);
v=-sort(-v);
alpha=u'*v;
end     % of function minscalarproductinnew



end   % of main testsp.... function
