%%*****************************************************************
%% HKMcorr: corrector step for the HKM 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] = HKMcorr(blk,dX,dZ,Avec,rp,Rd,sigmu,hRd,...
           permA,isspAy,nzlistAy,schur,UU,VV,X,Z,Zinv,dd,ee,gamz2);

    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,:); 
       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');       
          aa = mexqops(pblk,Z{p},dX{p},1); 
          bb = mexqops(pblk,Z{p},dZ{p},2); 
          cc = mexqops(pblk,dZ{p},dX{p},1);
          dxdz = mexqops(pblk,aa,dZ{p},4) + mexqops(pblk,bb,dX{p},3) ...
                - mexqops(pblk,cc,Z{p},4);     
          dxdz = mexqops(pblk,1./gamz2{p},dxdz,3); 
          tmp = sigmu*ee{p}-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}); 
          tmp = Prod3(pblk,dX{p},dZ{p},Zinv{p},0); 
	  tmp = 0.5*(tmp+tmp'); 
          tmp = sigmu*Zinv{p}-X{p}-tmp;
          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),X{p},3) ...
                + mexqops(pblk,mexqops(pblk,dZ{p},X{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,X{p},dZ{p},Zinv{p},0); 
          tmp = 0.5*(tmp+tmp'); 
          dX{p} = EinvRc{p}-tmp;
       end
    end 
%%=========================================================
