function L=matrixrep(m,n,A,tracepres)
%function L=matrixrep(m,n,A,tracepres)
%% find the matrix representation for completely positive T constraints
%%  T maps Hn ---> Hm  with constraints T(Ai)=Bi for all i
%INPUT:
%     m,n     order of Hermitian inputs
%     A cell of size k for input matrices (size is k+1 as unital
%             constraint is added in run.m file.)
%     tracepres    is 1 if trace preserving constraint should be added
%                  is 0 otherwise
%OUTPUT:
%   L        matrix representation (full)
%      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Mar 22/14 ... code can be sped up by vectorizing the for kk= ...  loop
%%%      using a kronecker product ... some details are in previous versions
%i=sqrt(-1);     % better to use   1i
%%Apr 5/14 ... completed the sparse representation of the input
%%             i.e. still looping for each column but just saving
%               nonzeros as a sparse representation
sqrt2=sqrt(2);
k=size(A,1);  % number of Ai,Bi constraints
mn=m*n;
tmm=m*(m-1)/2;
tnm=n*(n-1)/2;
if tracepres
      nconstr=k*m^2+n^2;
      else
      nconstr=k*m^2;
end
%%%%%%%%%%%%%%%%%%%%testing
M=zeros(k,m^2);
for ii=1:k,
	M(ii,:)=sHvec(A{ii})';
end
%%%%%%%%%%%%%%%%%%%%testing

totaltemptd=zeros(n^2,3);
totaltemptr=zeros(k*mn*(mn-1)/2,3);
totaltemptri=zeros(2*k*(mn*(mn-1)/2)-k*(m*n*(n-1)/2),3);   %/????need better counting????
trstartd=0;
trstartr=0;
trstarti=0;
			ttempd=[];
			ttempr=[];
			ttempri=[];
for qq=1:mn,    % for each basis element - off diagonal done 'twice'
	for pp=1:qq,  % by column and down the columns
		%%find the block (ss,tt) for the nonzero in basis elem. (pp,qq) 
		ss=ceil(pp/m);
		tt=ceil(qq/m);
		pst=pp-(ss-1)*m;   % position of elem in block
		qst=qq-(tt-1)*m;   % position of elem in block
		%if ss==tt,   % doing this for ss==tt only for now????
		if pst<qst,  % position for SHvec
			pqupper=(qst-1)*(qst-2)/2+pst;
			pqupperi=tmm+(qst-1)*(qst-2)/2+pst;
		elseif qst<pst  % reversed roles  ????
			pqupper=(pst-1)*(pst-2)/2+qst;   %%???
			pqupperi=tmm+(pst-1)*(pst-2)/2+qst;   %%???
		else
			pqupper=m*(m-1)+pst;
			pqupperi=m*(m-1)+pst;
		end
		%end
		if pp<qq,
			offdiagcols=((qq-2)*(qq-1)/2)+pp;
		end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		if ss==tt,   % this includes pp=qq case
			stupper=m*(m-1)+ss;  % position in sHvec(A^i)
			stupperi=m*(m-1)+ss;
		      if pp<qq,
				%%%following is for imag
			      irows=(pqupperi:m^2:k*m^2)';
			      jcols=offdiagcols*ones(k,1);
			      svals=M(:,stupperi);
		              ttempri=[irows jcols svals];
%%				%%%following is for real ... move up eventually
			      irows=(pqupper:m^2:k*m^2)';
			      jcols=offdiagcols*ones(k,1);
			      svals=M(:,stupper);
		              ttempr=[irows jcols svals];
		      else   % pp==qq
			      irows=(pqupper:m^2:k*m^2)';
			      jcols=pp*ones(k,1);
			      svals=M(:,stupper);
		              ttempd=[irows jcols svals];
		      end
	       else   %%%%%%%%%%%%%%%%ss<tt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ss==tt
			stupper=(tt-1)*(tt-2)/2+ss;
			stupperi=tmm+(tt-1)*(tt-2)/2+ss;
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			if pst == qst,    % entry is real 
			      jcols=offdiagcols*ones(k,1);
			      irows=(pqupper:m^2:k*m^2)';
			      svals=M(:,stupper);
		              ttempr=[irows jcols svals];
			      irowsi=(pqupperi:m^2:k*m^2)';
			      svalsi=-M(:,stupperi);
		              ttempri=[irowsi jcols svalsi];
		        elseif pst<qst,  
			     irows=[(pqupper:m^2:k*m^2)';(pqupperi:m^2:k*m^2)'];
			      jcols=offdiagcols*ones(2*k,1);
			      svals=[M(:,stupper);M(:,stupperi)]/sqrt2;
		              ttempr=[irows jcols svals];
			      irowsi=irows;
			      svalsi=[-M(:,stupperi);M(:,stupper)]/sqrt2;
		              ttempri=[irowsi jcols svalsi];
		        elseif qst<pst,  
			      jcols=offdiagcols*ones(2*k,1);
			     irows=[(pqupper:m^2:k*m^2)';(pqupperi:m^2:k*m^2)'];
			      svals=[M(:,stupper);-M(:,stupperi)]/sqrt2;
		              ttempr=[irows jcols svals];
			      irowsi1=(pqupper:m^2:k*m^2);
			      irowsi2=(pqupperi:m^2:k*m^2);
			      irowsi=[irowsi1';irowsi2'];
			      svalsi=[-M(:,stupperi);-M(:,stupper)]/sqrt2;
		              ttempri=[irowsi jcols svalsi];
		        end
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	       end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ss==tt
		%%% add the last constraint - trace preserving
	        if tracepres,
	           if pst==qst,    % diagonal elem of block is not zero
			if pp<qq,
				%need to still add pp==qq for ttempr
			      ttempr=[ttempr;k*m^2+stupper offdiagcols 1];
			      ttempri=[ttempri;k*m^2+stupperi offdiagcols 1];
		        else
			      ttempd=[ttempd;k*m^2+stupper pp 1];
		        end
		   end
	        end    % end of if trace preserving
		if pp==qq,
		        trendd=trstartd+length(ttempd);
			totaltemptd(trstartd+1:trendd,:)=ttempd;
			trstartd=trendd;
		else
		        trendr=trstartr+length(ttempr);
		        trendi=trstarti+length(ttempri);
			totaltemptr(trstartr+1:trendr,:)=ttempr;
			totaltemptri(trstarti+1:trendi,:)=ttempri;
			trstartr=trendr;
			trstarti=trendi;
		end

	end
end
if size(totaltemptri,1)~=trendi,
fprintf('at end size error in totaltemptri - still to fix');     %%?????
totaltemptri=totaltemptri(1:trendi,:);
end
if size(totaltemptd,1)~=trendd,
fprintf('at end size error td')
keyboard
end
if size(totaltemptr,1)~=trendr,
fprintf('at end size error tr')
keyboard
end
finalsparsematd=...
sparse(totaltemptd(:,1),totaltemptd(:,2),totaltemptd(:,3),...
     nconstr,mn,length(totaltemptd));
finalsparsematr=...
sparse(totaltemptr(:,1),totaltemptr(:,2),totaltemptr(:,3),...
     nconstr,(mn^2-mn)/2,length(totaltemptr));
finalsparsematri=...
sparse(totaltemptri(:,1),totaltemptri(:,2),totaltemptri(:,3),...
     nconstr,(mn^2-mn)/2,length(totaltemptri));
     %figure(1)
     %normest(finalsparsematri-Lcolsi)
     %spy(finalsparsematri-Lcolsi)
L=[finalsparsematr finalsparsematri finalsparsematd];
%L=[finalsparsematr Lcolsi finalsparsematd];
%L=[Lcolsr Lcolsi Lcolsd];
%% condition number check is too expensive ... skip for now???
%fprintf(' condition number of full matrix repres. L is %g \n',cond(full(L)));
