%%*******************************************************************
%% HSDHKMrhsfun: compute the right-hand side vector of the
%%               Schur complement equation for the  HKM direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [rhs,EinvRc,hRd] = HSDHKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ)

m = par.m;
if (nargin > 8)
    corrector = 1;
else
    corrector = 0;
    hRd = zeros(m+2,1);
end
hEinvRc = zeros(m+2,1);
EinvRc = cell(size(blk,1),1);
if length(sigmu)==1; sigmu = sigmu*ones(1,size(blk,1)); end
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});
    if strcmp(pblk{1},'l')
        if (corrector)
            Rq = dX{p}.*dZ{p};
        else
            Rq = sparse(n,1);
            tmp  = par.dd{p}.*Rd{p};
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)./Z{p}-X{p} -Rq./Z{p};
        tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'q')
        if (corrector)
            ff{p} = qops(pblk,1./par.gamz{p},Z{p},3); %#ok
            hdx = qops(pblk,par.gamz{p},ff{p},5,dX{p});
            hdz = qops(pblk,par.gamz{p},ff{p},6,dZ{p});
            hdxdz = Arrow(pblk,hdx,hdz);
            Rq = qops(pblk,par.gamz{p},ff{p},6,hdxdz);
        else
            Rq  = sparse(n,1);
            tmp = par.dd{p}.*Rd{p} ...
                + qops(pblk,qops(pblk,Rd{p},par.Zinv{p},1),X{p},3) ...
                + qops(pblk,qops(pblk,Rd{p},X{p},1),par.Zinv{p},3);
            tmp2 = mexMatvec(At{p,1},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)*par.Zinv{p}-X{p} -Rq;
        tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'s')
        if (corrector)
            Rq = Prod3(pblk,dX{p},dZ{p},par.Zinv{p},0);
            Rq = 0.5*(Rq+Rq');
        else
            Rq = sparse(n,n);
            tmp = Prod3(pblk,X{p},Rd{p},par.Zinv{p},0,par.nzlistAy{p});
            EinvRc{p} = 0.5*(tmp+tmp');
            tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)*par.Zinv{p}-X{p}- Rq;
        tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
        hEinvRc = hEinvRc + tmp2;
    end
end
%%
rhs = rp + hRd - hEinvRc;
rhs(m+1) = rhs(m+1) + (par.mu/par.tau - par.kap);
if (corrector)
    rhs(m+1) = rhs(m+1) - par.dtau*par.dkap/par.tau;
end
%%*******************************************************************
