%%************************************************************************
%% NTcorr: corrector step for the NT direction. 
%%
%% SDPT3: version 3.0 
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%************************************************************************

  function [dX,dy,dZ,solve_ok] = NTcorr(blk,dX,dZ,Avec,rp,Rd,sigmu,hRd,permA,isspAy,...
                        nzlistAy,schur,UU,VV,X,Z,W,G,sv,dd,ee,rr,vv,gamx,gamz); 

    global spdensity chol_options; 
%%
    m = length(rp); 
    EinvRc = cell(size(blk,1),1); hEinvRc = zeros(m,1); 
    for p = 1:size(blk,1)
       pblk = blk(p,:);  n = sum(pblk{2});
       if strcmp(pblk{1},'l');
          tmp = (sigmu - dX{p}.*dZ{p})./Z{p} -X{p};  
          EinvRc{p} = tmp; 
          hEinvRc = hEinvRc + (tmp'*Avec{p})'; 
       elseif strcmp(pblk{1},'q');       
          numblk = length(pblk{2});
          s = 1 + [0 cumsum(pblk{2})]; idx1 = s(1:numblk); 
          aa = mexqops(pblk,rr{p},dX{p},1); 
          bb = mexqops(pblk,rr{p},dZ{p},2); 
          cc = (aa.*dZ{p}(idx1)-bb.*dX{p}(idx1))./(1+rr{p}(idx1)); 
          hdxdz = mexqops(pblk,aa,dZ{p},3) + mexqops(pblk,bb,dX{p},3) ...
                 - mexqops(pblk,cc,rr{p},3); 
          hdxdz(idx1) = mexqops(pblk,dX{p},dZ{p},1); 
          tmp = mexqops(pblk,hdxdz,vv{p},2)./(gamx{p}.*gamz{p}); 
          Vihdxdz = mexqops(pblk,1./vv{p}(idx1),hdxdz,3) ...
                    - mexqops(pblk,tmp./vv{p}(idx1),vv{p},3); 
          Vihdxdz(idx1) = tmp; 
          tmp = mexqops(pblk,rr{p},Vihdxdz,2);
          dxdz = Vihdxdz - mexqops(pblk,(Vihdxdz(idx1)+tmp)./(1+rr{p}(idx1)),rr{p},3); 
          dxdz(idx1) = tmp; 
          dxdz = mexqops(pblk,sqrt(gamx{p}./gamz{p}),dxdz,3); 
          tmp = mexqops(pblk,-sigmu./(gamz{p}.*gamz{p}),Z{p},4)-X{p}-dxdz;
          EinvRc{p} = tmp;  
          hEinvRc = hEinvRc + (tmp'*Avec{p})'; 
       elseif strcmp(pblk{1},'s')
          numblk = length(pblk{2}); n2 = sum(pblk{2}.*pblk{2}); 
          hdZ = Prod3(pblk,G{p},dZ{p},G{p}',1); 
          hdX = spdiags(-sv{p},0,n,n)-hdZ;          
          tmp = Prod2(pblk,hdX,hdZ,0);  clear hdX hdZ; 
	  tmp = 0.5*(tmp+tmp');  
          if (numblk == 1) 
             d = sv{p};
   	     e = ones(pblk{2},1); 
  	     Rq = 2*tmp./(d*e'+e*d'); clear tmp; 
             if (nnz(Rq) <= spdensity*n2); Rq = sparse(Rq); end
          else
             Rq = sparse(n,n);
             s = [0 cumsum(pblk{2})]; 
             for i = 1:numblk
                 pos = [s(i)+1 : s(i+1)]; 
                 d = sv{p}(pos); e = ones(length(pos),1); 
                 Rq(pos,pos) = 2*tmp(pos,pos)./(d*e' + e*d'); 
             end;
          end
          tmp = spdiags(sigmu./sv{p}-sv{p},0,n,n); 
          tmp = Prod2(pblk,G{p}',tmp-Rq,0);  clear Rq; 
          tmp = Prod2(pblk,tmp,G{p},1);  
          EinvRc{p} = tmp;
          tmp2(permA(p,:),1) = mexinprod(blk,Avec,svec(pblk,tmp),m,p); 
          hEinvRc = hEinvRc + tmp2; 
       end
    end
    clear tmp tmp2; 
    rhs = rp + hRd - hEinvRc;
%%
    if (chol_options == 'matlab'); 
       [dy,solve_ok] = linsysolve(schur,UU,VV,rhs,0); 
    elseif (chol_options == 'lipsol')
       [dy,solve_ok] = linsysolve(schur,UU,VV,rhs,0);  
    end
%%
   for p=1:size(blk,1)
       pblk = blk(p,:); 
       if strcmp(pblk{1},'l')
          dZ{p} = Rd{p} - Avec{p}*dy;     
          dX{p} = EinvRc{p} - dZ{p}.*X{p}./Z{p};
       elseif strcmp(pblk{1},'q'); 
          dZ{p} = Rd{p} - Avec{p}*dy;  
          tmp = dd{p}.*dZ{p} + mexqops(pblk,mexqops(pblk,dZ{p},ee{p},1),ee{p},3);
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'s'); 
          numblk = length(pblk{2}); n2 = sum(pblk{2}.*pblk{2}); 
          dZ{p} = Rd{p} - smat(pblk,Avec{p}*dy(permA(p,:)),isspAy(p)); 
          tmp = Prod3(pblk,W{p},dZ{p},W{p},1); 
          dX{p} = EinvRc{p}-tmp;
       end
    end 
%%========================================================================





