
function [flag,Lappr,indsnzoR,indsnzoC]= rpcasolve(Yfinal,ZsClique,Ztemp,ZClique,LClique,Zpart,Lorig,m,n,opts,params,itergrow)
%%%% checking if Yfinal has 0 on the diagonal and solve the 
%%%% small problem by removing the zero cols and rows
% flag = 0 -solved
% flag = 1 - unsolved

start_time = tic;
tolerrank=params.tolerrank;
tolrankp = opts.tolrankp;  % for computing rank of exposing vector
tolrankq = opts.tolrankq;  % for computing rank of exposing vector
sqrt2 = sqrt(2);
  indszeroY=find((full(diag(Yfinal)==0)));
   numzerosYfinal=length(indszeroY);
   indszeroR = find((full(diag(Yfinal(1:m,1:m))==0)));
   indsnzoR = setdiff(1:m,indszeroR);
   indszeroC = find((full(diag(Yfinal(m+1:m+n,m+1:m+n))==0)));
   indsnzoC = setdiff(1:n,indszeroC);
   numzerosR=length(indszeroR);
   numzerosC=length(indszeroC);
   mm=length(indsnzoR);
   nn=length(indsnzoC);
  if numzerosYfinal > 0
     %indsnonzero=setdiff(1:length(Yfinal),indszeroY);
     %spy(Yfinal( [indsnonzero indszeroY'],[indsnonzero indszeroY']))
      fprintf('WARNING:  number of total zero  %i >0 \n', numzerosYfinal)
       fprintf('WARNING:  number of zero rows  %i >0 \n', numzerosR)
         fprintf('WARNING:  number of zero column  %i >0 \n', numzerosC)
      fprintf('WARNING:  shifting out diagonal \n')
    
  end
  Znzo = Ztemp(indsnzoR,indsnzoC);
  ZpartN = Zpart(indsnzoR,indsnzoC);
  LCliqueN= LClique(indsnzoR,indsnzoC);
  ZCliqueN = ZClique(indsnzoR,indsnzoC);
  ZsCliqueN = ZsClique(indsnzoR,indsnzoC);
  LorigN = Lorig(indsnzoR,indsnzoC);
  [YfinalN,ZCliqueN,ZsCliqueN,LCliqueN]= YfindExp(Znzo,ZpartN,ZCliqueN,ZsCliqueN,LCliqueN,LorigN,mm,nn,params,opts,itergrow);
  indsZN = find(Znzo>0);
  indsZcN = find(ZCliqueN>0);
  indsZcNL = find(ZsCliqueN>0);
  indSCN = setdiff(indsZN,indsZcNL);
  
%Yfinalnz = zeros(numnzoR+numnzoC,numnzoR+numnzoC);
%Yfinal = Yfinal ;
Ypn= Yfinal(indsnzoR,indsnzoR) + YfinalN(1:mm,1:mm);
Yqn= Yfinal(indsnzoC+m,indsnzoC+m) + YfinalN(mm+1:end,mm+1:end);
if opts.spn ==0
Vpn=null(full(Ypn));                 %%?????change to sparse???????
Vqn=null(full(Yqn));
else
   % Vpn1=null(full(Ypn));                 %%?????change to sparse???????
% Vqn1=null(full(Yqn));
%Vpn=nulls(Ypn);                 %%?????change to sparse???????
%Vqn=nulls(Yqn);
%Vpn = full(spspaces(Ypn,2,tolrankp));
%Vqn = full(spspaces(Yqn,2,tolrankq));
% 
%   eigsOpts.issym=1;
%   eigsOpts.isreal=1;
%   %eigsOpts.tol=1e-29;
%   eigsOpts.v0=randn(m,1);
%   eigsOpts.p=100;
 try 
     rd =1;
     np = 40;
     nq = 40;
     while rd == 1 
  numbereigsp=min(mm-1,np);
  numbereigsq=min(nn-1,nq);
 warning('off','all');
 try
  [UYp,eYp,~] = eigs(Ypn,numbereigsp,'sm');%,eigsOpts); 
 catch
       [UYp,eYp,~] = eigs(Ypn+eps*eye(mm),numbereigsp,'sm');%,eigsOpts);
 end
 try
  [UYq,eYq,~] = eigs(Yqn,numbereigsq,'sm');%,eigsOpts); 
 catch
     [UYq,eYq,~] = eigs(Yqn+eps*eye(nn),numbereigsq,'sm');
 end
  deYp=diag(eYp);
  deYq=diag(eYq);
  %deY=deY(1:end);
  % rp = find(abs(deYp) >= tolrankp, 1, 'first');
   indrp = abs(deYp) < tolrankp;
   rp = find(abs(deYp) >= tolrankp, 1);
   %rq = find(abs(deYq) >= tolrankq, 1, 'first');
   indrq = abs(deYq) < tolrankq;
   rq = find(abs(deYq) >= tolrankq, 1);
  % Vpn=UYp(:,1:rp-1);   % for eigs
  % Vqn=UYq(:,1:rq-1);
   Vpn=UYp(:,indrp);   % for eigs
   Vqn=UYq(:,indrq);
   warning('on', 'all');
   rd = 0;
%%%%%% End: Aug 11 changed %%%%%%
 if isempty(rp) && np < mm
     np = np + 40;
     rd = 1;
 end
 if  isempty(rq) && nq < nn
     nq = nq + 40;
     rd = 1;
 end
  %%%% finding the final exposing vector %%%%
  
  end
 catch
     Vpn=[];                 %%?????change to sparse???????
     Vqn=[] ;
    % keyboard;
 end
%{
[Vpn,Lpn] = eig(full(Ypn));   
[Vqn,Lqn] = eig(full(Yqn));   
indspn=find(diag(Lpn)>tolrankp);
indsqn=find(diag(Lqn)>tolrankq);
rpn=mm-indspn(1)+1;
rqn=nn-indsqn(1)+1;

rpn = rpn -opts.redrank;
rqn = rqn -opts.redrank;


Vpn = Vpn(:,1:mm-rpn);          % nullspace eigenspace
Vqn = Vqn(:,1:nn-rqn);  
%}
end

  Vn = [Vpn   zeros(mm,size(Vqn,2))
       zeros(nn,size(Vpn,2))   Vqn];
  rnn = size(Vn,2);
  tvn = rnn*(rnn+1)/2;
  
   
 % Imn=eye(mm+nn);
%Zs = spalloc(mm,nn,length(indSCN));
%Zs = false(mm,nn);
%Zs(indSCN) = 1;
%indCN=find([Zs;sparse(nn,nn)])+(mm+nn)*mm;  % indC not bicliques
%A=[(true(mm))   Zs
 %          Zs'     (true(nn))];
 %A(Imn==1)=0;  
 %[indci,indcj]=ind2sub([mm+nn,mm+nn],indCN);
       %indsZc = find(abs(Zorig)>0);

%Zc = spalloc(mm,nn,length(indsZcN));
%Zc = false(mm,nn);
 %   Zc(indsZcN) = 1;
  %  indscN=find([Zc;sparse(nn,nn)])+(mm+nn)*mm; % for matrix rep: ALL bicliques
  %  A=[(true(mm))   Zc
    %       Zc'     (true(nn))];
   %     A(Imn==1)=0;  
 %[indsi,indsj]=ind2sub([mm+nn,mm+nn],indscN);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %indsdiag=(1:rnn)';
  %indsdiag=indsdiag.*(indsdiag+1)/2;
  %indsoffdiag=setdiff(1:tvn,indsdiag);
  
  %E=ones(rnn);
  %[mi,mj]=find(triu(E,1));
  %indsu=find(triu(E,0));
  %indsuu=find(triu(E,1));

  %%% find matrix representations
  %%%  matreps in bicliques; matrrepsc NOT in bicliques
  % tempYc=[sparse(mm,mm) ZpartN;sparse(nn,mm+nn)];
  %tempY=[sparse(mm,mm) LCliqueN;sparse(nn,mm+nn)];
  % allzinc=tempY(indscN); % vector of known entries of Z in bicliques
  % allznotc = tempYc(indCN); % vector of known entries of Z not incliques
   allzinc = LCliqueN(indsZcN); % vector of known entries of Z in bicliques
   allznotc = ZpartN(indSCN); % vector of known entries of Z not incliques
  
  % clear tempY tempYc
  mp = size(Vpn,2);
  mq = size(Vqn,2);
  %M = [zeros(mp,mp),ones(mp,mq)
   %    ones(mq,mp), zeros(mq,mq)];
  %indcN= find(localsvec(M,indsu,indsuu))';
  % clear M
  %VT=Vn';
  %%% matrix representation for points in the  bicliques (4.5a) in paper
  %matreps=zeros(length(indscN),tvn);
  % for zz=1:length(indscN)   % all elements of bicliques are sampled
	%  Ri=VT(:,indsi(zz));
	 % Rj=VT(:,indsj(zz));
	 % Tdiag=Ri.*Rj; 
	 % Toffdiag = (Ri(mj).*Rj(mi)+Ri(mi).*Rj(mj))/sqrt2; 
	 % matreps(zz,indsdiag)=Tdiag;   
	 % matreps(zz,indsoffdiag)=Toffdiag;
  % end
  %matreps = matreps(:,indcN);    %  (4.5a) in bicliques in paper
  % indsZcN = indsZcN(1:mp*mq);
  % indSCN = indSCN(1:mp*mq);
  try
  matreps = zeros(length(indsZcN),mp*mq);
  matrepsc = zeros(length(indSCN),mp*mq);
  catch
      flag = 1;
      Lappr = [];
      return
  end
  IV = eye(mp*mq);
  %matreps = sparse(length(indsZcN),mp*mq);
  %matrepsc = sparse(length(indSCN),mp*mq);
  %IV = speye(mp*mq);
  for i = 1:mp*mq
      IVC = IV(:,i);
   Temp = Vpn*reshape(IVC,mp,mq)*Vqn';
   matreps(:,i) = Temp(indsZcN);
   matrepsc(:,i) = Temp(indSCN);
  end
  
     [matrepssub,indssub]=licols(matreps',tolerrank); %in bicliques
      matrepssub = matrepssub';
      allzsub =  allzinc(indssub);  
  %%%% matrix rerpesntation for points not in the bicliques (4.5b) in paper
 % %{
  if opts.model == 1 &&  size(matrepssub,1) ~= size(matrepssub,2)
 % matrepsc=zeros(length(indCN),tvn);
  %for zz=1:length(indCN)
	%  Ri=VT(:,indci(zz));
	 % Rj=VT(:,indcj(zz));
	 % Tdiag=Ri.*Rj; 
	 % Toffdiag = (Ri(mj).*Rj(mi)+Ri(mi).*Rj(mj))/sqrt2; 
	 % matrepsc(zz,indsdiag)=Tdiag;   
	 % matrepsc(zz,indsoffdiag)=Toffdiag;
  % end
  %matrepsc = matrepsc(:,indcN);    %  (4.5b) not in bicliques in paper
  
  %%%%%%%%%%%%
  %%%???? why is matrepsc not full col rank?????? underdeterm????
  %%%%%% need to look at underdetermined sketch matrix applic????
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % recover the big matrix Z by nuclear norm minimization
   
   
   
      %%%%%%%%%%%%%
      tic
      rc = rank(matrepssub);
      [L,U,P ]= lu([matrepssub',matrepsc']);
      U1 = U(rc+1:end,:);
      U2 = U1(:,size(matrepssub,1)+1:end);
      indsc = find(vecnorm1(U2)>tolerrank);
      toc
%       tic
%       [L0,U0,P0] = lu(matrepssub');
%        rd = rank(U0);
%       Q = (L0)\(P0*(matrepsc'));
%       indsc2 = find(vecnorm(Q(rd+1:end,:))>tolerrank);
%       toc
%       setdiff(indsc1,indsc2);
     %matrepsc = matrepsc(indsc,:);
     % keyboard
      matrepscsub =  matrepsc(indsc,:);
      allzcsub = allznotc(indsc,:);
      %%%%%%%%%%%%%factor out nullspace for underdetermined system
      [~,smtran,vmtran]= svd(matrepscsub,0);  
      rvmtran= rank(smtran);
      rgmattran=vmtran(:,1:rvmtran);   % basis range of transpose
      nullmattran=vmtran(:,rvmtran+1:end);   % basis of nullspace
      rgmatrepscsub = matrepscsub*rgmattran; % for obj function l1 norm
      %  muu = 1/(sqrt(mq*mp)*params.p)%1*params.muu;
        muu = params.muu
	%%%%!!!??? uncomment if sketch used?????
	%Sk=randn(8*((mp*mq)),size(matrepscsub,1));  %% ??? what size???
        %Skallzcsub = Sk*allzcsub;
        %Skrgmatrepscsub=Sk*rgmatrepscsub;
   %%%%%%%%%%%
      matrepssrg= matrepssub*rgmattran;
      indrg=find(vecnorm1(matrepssrg)>1e-10);  % nonzero cols indices
      matrepssrg=matrepssrg(:,indrg);
   %%%%%%%%%%%
      matrepssrn=matrepssub*nullmattran;
      indrn=find(vecnorm1(matrepssrn)>1e-10);  % nonzero cols indices
      if isempty(indrn)
	      fprintf('size matrepssrn is %i and  set indrn is empty for nullspace!! \n',size(matrepssrn,2))
      end
      matrepssrn=matrepssrn(:,indrn);
   %%%%%%%%%%%

   fprintf('rng bicliques matrepssrg %i by %i ; nullsp bicliques matrepssrn %i by %i\n',...
   size(matrepssrg),size(matrepssrn))
   fprintf('objective fn rgmatrepscsub %i by %i \n',...
           size(rgmatrepscsub))
  end
    %for iii=1:1
    %}
     if size(matrepssub,1)== size(matrepssub,2)
	        Rv = matrepssub\allzsub;
            R = reshape(Rv,mp,mq);
            fprintf('CVX is not used!! \n')
        
      else
    fprintf('starting cvx loops with R %i by %i \n',mp,mq)
      startcvx=tic;
      if opts.model == 1
      try
      cvx_clear
      cvx_quiet true
      %cvx_solver sedumi
      cvx_begin sdp
      %cvx_solver sdpt3
      cvx_solver mosek
      cvx_precision best 
         variable R(mp,mq)
         variable Rrng(size(rgmattran,2),1)  % basis for range of transp.
         variable Rnull(size(nullmattran,2),1) % basis for nullspace
   %%% reshape for R and get rid of R???? %%% use sketch  Sk or not????
   %minimize (norm_nuc(rgmattran*Rrng+nullmattran*Rnull) + muu*norm(allzcsub-rgmatrepscsub*Rrng,1))
   minimize (norm_nuc(R) + muu*norm(allzcsub-rgmatrepscsub*Rrng,1))
   %minimize (norm_nuc(R) + muu*norm(Skallzcsub-Skrgmatrepscsub*Rrng,1))
         subject to
	    %norm_nuc(R) <= 2*muu*1e3; 
    	R(:) == rgmattran*Rrng+nullmattran*Rnull;
       if isempty(indrn)
                      matrepssrg*Rrng(indrg) == allzsub;
	    else
	              (matrepssrn*Rnull(indrn))+...
                         (matrepssrg*Rrng(indrg)) == allzsub;
	    end
            %matrepssub*R(:) == allzsub;   % for bicliques
      cvx_end
      catch
          flag = 1;
        keyboard; 
      end     % end of try/catch
      elseif opts.model == 0
          try
          cvx_clear
      cvx_quiet true
      %cvx_solver sedumi
      cvx_begin sdp
      %cvx_solver sdpt3
      cvx_solver mosek
      cvx_precision best 
         variable R(mp,mq)
   %%% reshape for R and get rid of R???? %%% use sketch  Sk or not????
   %minimize (norm_nuc(rgmattran*Rrng+nullmattran*Rnull) + muu*norm(allzcsub-rgmatrepscsub*Rrng,1))
   minimize norm_nuc(R) 
   %minimize (norm_nuc(R) + muu*norm(Skallzcsub-Skrgmatrepscsub*Rrng,1))
         subject to
	    %norm_nuc(R) <= 2*muu*1e3; 
    	matrepssub*R(:) == allzsub;
       
            %matrepssub*R(:) == allzsub;   % for bicliques
      cvx_end
    % keyboard
          catch
            %  keyboard;
          end
      end
      time2=toc(startcvx);
      fprintf('time for cvx is  %g   \n',time2)
      if isnan(R)
             fprintf('R is NAN\n')
            % keyboard
             flag =1;
      end
    end
     % Zorig = Zorig(indsnzoR,indsnzoC);
      R = full(R);
      Lappr = (Vpn*R*Vqn');  % approx for Lorig
      Lappr = round(Lappr);
   %Snoclq=allzcsub-rgmatrepscsub*Rrng;  % error Zorig for sampled/non biclique
   %Lappr = round(Lappr);
   %Sappr=round(Zorig-Lappr);   % new sparsity part 'cheating'
   %Zappr=Sappr+Lappr;
   %% relative error for sampled elements for Lorig
   %resids= norm(Lappr-Lorig,'fro')/norm(Lorig,'fro');
%    resids(2) = norm(Zappr-Zorig,'fro')/norm(Zorig,'fro');
%    fprintf('abs and relat. residual low rank total error %g  %g \n',...
%                  norm(Lappr-Lorig,'fro'),...
%                  resids(1))
 %  fprintf('abs and relat. residual low rank error sampled data %g  %g \n',...
  %              norm(Lappr(indsZcN)-LClique(indsZcN),'fro'),...
   %             norm(Lappr(indsZcN)-LClique(indsZcN),'fro')/norm(LClique(indsZcN),'fro'));
  % fprintf('abs and relat. residual sparse sampled error %g  %g \n',...
   %        norm(Sappr(indsZN)-Sorig(indsZN),'fro'),...
    %       norm(Sappr(indsZN)-Sorig(indsZN),'fro')/max(norm(Sorig(indsZN),'fro'),1))
   %fprintf('abs and relat. residual sparse total error %g  %g \n',...
    %       norm(Sappr-Sorig,'fro'),...
     %      norm(Sappr-Sorig,'fro')/max(norm(Sorig,'fro'),1))
     %%%%%
     %muu=.7*muu+.3*2;
     %end
     %  norm(round(Lappr(indsZcN))-LClique(indsZcN),'fro') == 0 &&...
     if   rank(R/norm(R,2),params.tolerrank) == params.r
     flag = 0;
     else
      flag = 1;
     end
 
time1=toc(start_time);
fprintf('Time for file rpcasolve   %g\n',time1)

%%%%%%%%%%%%%%%temp output for now?????????????
 %keyboard;
 clear matreps matrepsc
 %rn = size(V,2);   % output rn done above for initial rank
   % end of function completeZ.m for exposvct
end
