function P=rankaltproj(nn,toler,iterlimit,P0,rankP,L,R)
%function P=rankaltproj(nn,toler,iterlimit,P0,rankP,L,R)
%INPUT: 
%  nn,mm  input size of unknown matrix P (nn by nn) and number of constrs
%  toler desired stopping tolerance   default .01
% iterlimit   iteration limit    default 30
% P0          initial estimate for facial reduced 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 nn)
% L,R        current matrix representation and RHS with mm rows
%             This allows a matrix represt to be passed that reduces the
%             rank of a previous obtained optimal solution
%%%%%%
%OUTPUT: P      Hermitian size m*n for compl pos. T
%                redundant constraints but before changing column order
%                This is needed when trying to reduce rank
%%%%%%%%%%%%%%%%%%%
%%



fromstart=cputime;   % time complete algorithm
%%%%%%%%%%%%%Initialize
tnnm=nn*(nn-1)/2;
svindu=logical(triu(ones(nn),1));
sqrt2=sqrt(2);
if ~exist('toler','var')
	toler=.01;
end
if ~exist('iterlimit','var')
        iterlimit=30;
end
if ~exist('rankP','var'),
	rankP=nn;
	fprintf('setting rankP to nn =  %i  \n',nn);
end
if ~exist('P0','var')   ||   size(P0,1)==0,
        P0=nn*blkdiag(eye(rankP),zeros(nn-rankP));  % initial starting point
	fprintf('using [nn*I_r 0] as initial starting point for P \n')
end
fprintf('solving using alt. proj. with fixed rank <=   %i \n',...
     rankP)
P=P0;
clear P0;
%find the matrix representation of the constraints and RHS
%startLtime=cputime;

residual=R-(L*sHvec(P));
noiter=0;
norms=[];
cosines=[];
normPdiffold=norm((P));    % using P_-1 = 0
%%%%%%%%%%%%%%%%%%%test permutation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%permute rows first%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
Lp=L;    % no sparsity expected for noww/?????
%pcols=colamd(L);    % permutation of cols for sparsity
%%Lp=L(:,pcols);
%%figure(2)
%%spy(Lp)
%%xlabel(['matrix repr. L after row permut. and colamd, mm,nn is  ' num2str([mm nn])])
%%[~,pcolsinv]=sort(pcols);   % inverse perm for cols
%%%%%%%%%%%%%%%%%%%testing  permute rows first%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% temporary fix!!!!! for roundoff ... get/find structure!!!
%%fprintf('temporarily!!!! permute cols using colamd for sparsity!!!! \n')
%%%pinvL=(pinv(full(L(:,pcols))));  % what if L is sparse ????? and it is!!!!
%%%pinvL=Lp'*inv(Lp*Lp');
pinvL=(inv(Lp*Lp')*Lp)'; % ' at end appears faster
%fprintf('time for colamd and pinv(L) is  %g   \n',toc);
indstemp=find((abs(pinvL)<1e-10)& (abs(pinvL)>0));
fprintf('temporarily!!! removing %i small elements in pinvL - roundoff \n',...
           length(indstemp));
pinvL(indstemp)=0;
%figure(3)
%spy(pinvL)
%xlabel(['gen inv after removing roundoff, mm,nn is  ' num2str([mm nn])]);
%%Lt=L*L';
%%figure(4)
%%spy(Lt)
%%title('L L'' gives ')
%if ~issparse(pinvL),
%     fprintf('changing pinvL to sparse(pinvL)')
%     pinvL=sparse(pinvL);
%end
%%%% temporary fix!!!!! for roundoff ... get/find structure!!! of L!!!
fprintf('nnz of L, pinvL and density pinvL are  %i %i %g \n',...
       nnz(L),nnz(pinvL),nnz(pinvL)/numel(L));
%%%possible 'large' sparse case for pinvL???
%LL=sparse(L); [uu,ss,vv]=svds(LL,size(LL,1)); pinvsvd=vv*(ss\uu');
feasdetect=0;
fprintf('\n iter norms-ratio     cos-vectors   norm-residual eig-proj \n')
startwhiletime=cputime;   % start iterations
while noiter < iterlimit &&...
        (norm(residual))>toler && feasdetect<4,
	noiter=noiter+1;
	Pvec=pinvL*residual; % Pl-P;
	%Pvec=Pvec(pcolsinv); % using inverse cols perm. for perm in pinvL
	% construct Pdiffold=sHMat(Pvec) here
	Sympart=zeros(nn);
	Sympart(svindu)=Pvec(1:tnnm)/sqrt2;
	Sympart=Sympart+Sympart';
	Skewpart=zeros(nn);
	Skewpart(svindu)=Pvec(tnnm+1:2*tnnm)/sqrt2;
	Skewpart=Skewpart-Skewpart';
	Pdiffold=Sympart+1i*Skewpart+diag(Pvec(2*tnnm+1:end)); % 1i is sqrt(-1)
	%Pdiffoldtemp2=...
	%     usHMat(Pvec(1:tmnm),mn) +...
	%           i*uskewHMat(Pvec(tmnm+1:2*tmnm),mn) +diag(Pvec(2*tmnm+1:end));
	%Pdiffoldtemp=sHMat(pinvL*residual); % Pl-P
	Pl=P+Pdiffold;  % projection onto linear manifold
	%%%check Pl is on linear manifold
	normPdiffnew=norm(Pl-P,'fro');
	norms=[norms normPdiffnew/normPdiffold];
	%fprintf('norm of residual after proj onto lin manifold  %g \n',...
	%                norm(full(R-(L*sHvec(Pl)))));
	%fprintf('ratio of norms after proj onto lin manifold  %g \n',...
	%                normPdiffnew/normPdiffold);
	[u,d]=eig(Pl);
	d=real(d);  % to avoid roundoff in eig calculation
	[dd,indsall]=sort(diag(d),'descend');   
        indsneg=find(dd<=0);
        % project using largest pos. rankP eigs
	indspos=indsall(1:min(nn-length(indsneg),rankP)); 
        P=u(:,indspos)*d(indspos,indspos)*u(:,indspos)';
	           % projection onto psd
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if size(P,1)==0,           % proj. of nsd is 0
		P=zeros(nn);
	end
	P=(P+P')/2;    
	%P=(nn/trace(P))*P;  % normalize for trace ?????needed? helps???
	Pdiffnew=P-Pl;
	cosvalue=real(trace((Pdiffold/norm(Pdiffold,'fro'))'*...
	                  (Pdiffnew/norm(Pdiffnew,'fro'))));
	cosines=...
	   [cosines cosvalue];
	RP=real(P);
	IP=imag(P);
	vP=[sqrt2*RP(svindu);sqrt2*IP(svindu);diag(P)];
        residual=R-(L*vP); % update residual
        %residual=R-(L*sHvec(P)); % update residual
	normratio=normPdiffnew/normPdiffold;
	%%%?????need to decide theoretically how to detect infeas
	%%%/??? what accuracy to use below
	if abs(1-normratio)+abs(-1-cosvalue)<1e-6, % not converging
	           feasdetect=feasdetect+1;
	end
	if 100*round(noiter/100)==noiter,
	  fprintf('%i  %10.4g  %14.4g  %14.4g  %8.4g \n',...
	  noiter,normPdiffnew/normPdiffold,cosvalue,norm(residual),norm(P-Pl));
	end
	normPdiffold=normPdiffnew;
	%%% check for separating hyperplane after ?? iterations???
	%%%%?????
end
fprintf('%i  %10.4g  %14.4g  %14.4g  %8.4g \n',...
         noiter,normratio,cosvalue,norm(residual),norm(P-Pl));
fprintf('time for iterations/while loop is  %g   \n',cputime-startwhiletime);
if feasdetect>=4,
	fprintf('algorithm shows - infeasibile???\n')
	fprintf('a separating hyperplane is ????\n')
end
%%%%%%%%%%%remove comments to plot ratios of norms and cosines
%figure(1)
%plot(norms)
%title('ratios of norms of consecutive iterates')
%figure(2)
%plot(cosines)
%title('cosines of angles between consecutive iterates')
fprintf('\nfinal norms ratio is  %g ; final cosine is  %g  \n',...
                  norms(end),cosines(end));
%%%%%%%%%%%%%%%%%verify/check constraints
fprintf('total time for for algorithm from rankaltproj.m file is  %g \n',...
cputime-fromstart)
