function [Poptfull,rankPs,Fstep]=rankredTPCP(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres,rankPtr);
% [Poptfull,rankPs,Fstep]=rankredTPCP(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres,rankPtr);
%INPUT: 
%  n,m  input matrix sizes of A,B
%  A,B  cells of size k containing input of Hermitian matrices
%        for constraints   T(Ai)=Bi  i=1...k
%        We assume the last of these constraint is  unital T(In)=Im
%        We are including the trace preserving constraint
%  toler desired stopping tolerance   default .01
% iterlimit   iteration limit    default 30
% P0          initial estimate for compl pos Choi matrix P, 
%                    use P0=[] if initial estimate not available
% rankP     desired rank for final P - means projection on faces of rank
%                                          rankP are used. (default mn)
% unital    is 1 (default) if unital constraint is included and 0 otherwise
%           If 1, then it is assumed that the last constraint in A is the
%           unital constraint A{k}=eye(n)/n;B{k}=eye(m)/m;
% tracepres is 1 (default) if trace preserving constraint 
%                    should be included and 0 otherwise.
%               If 0, then it is added in matrix representation in
%               matrixrep.m
% rankPtr    rank of initial TPCP used to generate data
%%%%%%
%OUTPUT: Poptful      Hermitian size m*n for compl pos. T optimal Ps
%        rankPs       ranks for the optimal Ps
%        Fstep        number of facial reductions + 2



%%% solve for first optimal TPCP P
mn=m*n;
rankPs{1}=rankPtr;
starttotal=cputime;
[P,L,R]=altproj(n,m,A,B,toler,iterlimit,P0,rankP,unital,tracepres);
endtotaltime=cputime-starttotal;
fprintf('total time from rankrunTPCP for first altproj steps is  %g   \n',endtotaltime)

%symmetrize optimal P to avoid roundoff error problems with eigs
P=(P+P')/2;
Poptful{1}=[];
fprintf('First optimal P=Poptfull{2}=Popt{2} is found in rankrunTPCP.m \n')
Popt{2}=P;    % first optimal solution from altproj
Us{1}=eye(mn);   % eigenvectors for Popt{1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Optimal initial P=Popt{2} has been found - assuming 'good' accuracy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%P=Popt{2} is opt from altproj using data LsHvec(P)=R,P>=0
nonunique=true;
Fstep=2; % start with P=Popt{2} from altproj
Lorig=L;Rorig=R;  % for error check at end
rankL=rank(full(L));
Ulows=eye(mn);
while nonunique,    % while optimal solution P is not unique
	fprintf('\n\nstarting Fstep %i in while loop \n',Fstep);
%Check for unique solution of linear system first  L sHvec(P)=R
% i.e. full column rank for L. ... 
fprintf('[ norm(L*sHvec(P) -R) and  min(eig(P))  are  %g %g \n',...
               norm(L*sHvec(P) -R),min(eig(P)));
[U,DP]=eig(P);
[lamsP,inds]=sort(diag(DP),'descend');
rankPs{Fstep}=rank(DP); % already diagonalized
DP=diag(lamsP(1:rankPs{Fstep}));   % feasible for facially reduced problem
Us{Fstep}=U(:,inds(1:rankPs{Fstep}));    % current U
if norm(P-Us{Fstep}*DP*Us{Fstep}')>1e-10,
	fprintf('norm error reconstr P=UDU* from *compact* eig decomp %g \n',...
	        norm(P-Us{Fstep}*DP*Us{Fstep}'));
	fprintf('enter return to continue \n')
	keyboard
else
	fprintf('P reconstr. using compact spectral decomp verified \n')
end
Ulows=Ulows*Us{Fstep};
Plowrank=Ulows*DP*Ulows';
Plowrank=(Plowrank+Plowrank')/2;
if norm(Lorig*sHvec(Plowrank)-Rorig)>1e-10,
	fprintf('norm error in reconstr Plowrank using Ulows %g \n',...
             norm(Lorig*sHvec(Plowrank)-Rorig));
	fprintf('enter return to continue \n')
	keyboard
else
	fprintf('facially reduced P feas. verified for original probl.  at Fstep = %i \n',Fstep)
end
Poptfull{Fstep}=Plowrank;  % save full matrices in cell 
%%%%%%%%%%%%%%%
%% update constraint matrix L
Lnew=zeros(size(L,1),rankPs{Fstep}^2);
for ii=1:size(L,1),
	Lnew(ii,:)=real(sHvec(Us{Fstep}'*sHMat(L(ii,:)')*Us{Fstep}));
end
Rnew=R;
residualface=Rnew-Lnew*sHvec(DP);
if norm(residualface)>1e-10,
	fprintf('norm error in feas. check on face %g \n',...
                 norm(residualface));
	fprintf('enter return to continue \n')
	keyboard
else
	fprintf('diagonal matrix(largest eigs P) is feas on face/verified \n')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rankLnew=rank(Lnew);
if rankLnew==size(Lnew,2),
	fprintf('Lnew full col. rank; current P unique; no further rank reduction possible \n')
	nonunique=false;
else
%%Set P0 indefinite to ensure starting outside cone
  nn=rankPs{Fstep};
  mm=size(Lnew,1);
  tt=max(1,floor(nn/2));
  ttn=nn-tt;
  P0=blkdiag(-eye(tt),eye(ttn));  % ensure indef start
  if rankLnew < rankL,
	%%%%%%%%%%%?????following heuristic is NOT!!! rigorous yet!!!
	fprintf('rank of Lnew %i < rank L %i - so decreasing size L \n',...
	            rankLnew,rankL);
	[ll,uu,pp]=lu(Lnew,'vector');
	Lnew=Lnew(1:rankLnew,:); % ??temp fix?? cannot always discard thse rows
	Rnew=R(1:rankLnew);
'before < rankLnew ???? for debugging???? not done yet????'
	keyboard
  end
starttotal=cputime;
P=rankaltproj(nn,toler,iterlimit,P0,rankP,Lnew,Rnew);
fprintf('total time from rankrunTPCP for altproj steps is  %g   \n',endtotaltime)
%%%%%%%%%%%updates
Fstep=Fstep+1;
Popt{Fstep}=P;
L=Lnew;
rankL=rank(L);
R=Rnew;
rankPs{Fstep}=rank(Popt{Fstep});
fprintf('full mn by mn optimal matrices are in cell Poptfull \n')
fprintf('\nranks for initial Ptr, first opt P and %i optimal Ps from facial reduction steps are  \n', (Fstep-2))
for ii=1:Fstep,
	fprintf('rank is %i \n',rankPs{ii});
end

end   % check on full column rank

end   % while nonunique


%%%%%%%%%%%%%%%check on equations for last full size facially reduced P
P=Poptfull{end};
k=length(A);
for kk=1:k,
	if norm(B{kk}-TP(m,n,P,A{kk}))>(max(m,n)*toler),
		if kk==k && unital,
		   fprintf('error for unital k= %i for ||TP(A_k)-B_k||,  is  %g \n',...
		   kk,norm(B{kk}-TP(m,n,P,A{kk})));
		else
		   fprintf('error for k= %i for ||TP(A_k)-B_k||,  is  %g \n',...
		   kk,norm(B{kk}-TP(m,n,P,A{kk})));
	        end
	end
end
if tracepres && norm(eye(n)-partialtrace(m,n,P))>((mn)*toler),
          fprintf('error for partial trace ||parttr(P)-I_n||,  is  %g \n',...
		   norm(partialtrace(m,n,P)-eye(n)));
end
fprintf('error check opt low rank P=Poptfull{%i} for orig. constrs. completed in rankredTPCP \n',length(Poptfull))
fprintf('%i optimal Ps are in cell Poptfull starting at Poptfull{2} \n',...
           (Fstep-1))
