%%*******************************************************************
%% HKMpred: Compute (dX,dy,dZ) for the H..K..M 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,hRd,schur,solve_ok,UU,VV,Zinv,dd,ee,gamz2] = ...
           HKMpred(blk,Avec,rp,Rd,sigmu,permA,isspA,isspAy,nzlistA,...
           nzlistAsum,nzlistAy,X,Z,invZchol);

    global spdensity iter chol_options nnzschurold; 
    global NNZSCHUR; 
%%
%%  
%%
    m = length(rp);
    schur = sparse(m,m); UU = []; VV = []; 
    dd = cell(size(blk,1),1);   ee = cell(size(blk,1),1); gamz2 = cell(size(blk,1),1); 
    Zinv = cell(size(blk,1),1); EinvRc = cell(size(blk,1),1);
    dX = cell(size(blk,1),1); dy = zeros(m,1); dZ = cell(size(blk,1),1); 
    hRd = zeros(m,1);  hEinvRc = zeros(m,1); 
    existsdpblk = 0; 
%%
%% compute schur matrix
%%
    for p = 1:size(blk,1)
        pblk = blk(p,:); 
        n = sum(pblk{2}); 
        if strcmp(pblk{1},'l')
           dd1 = X{p}./Z{p};
           tmp = Rd{p}.*dd1;
           hRd = hRd + (tmp'*Avec{p})'; 
           EinvRc{p} = sigmu./Z{p}-X{p}; 
           hEinvRc = hEinvRc + (EinvRc{p}'*Avec{p})';  
           [decolidx,spcolidx] = checkdense(Avec{p}'); 
           if ~isempty(decolidx); 
              len = length(decolidx); 
              tmp = (spdiags(sqrt(dd1(decolidx)),0,len,len)*Avec{p}(decolidx,:))'; 
              UU = [UU tmp]; 
              VV = [VV tmp]; 
              dd2 = zeros(n,1);  
              dd2(spcolidx) = dd1(spcolidx);   
              schur = schur + Avec{p}' *spdiags(dd2,0,n,n)*Avec{p};
           else
              schur = schur + Avec{p}' *spdiags(dd1,0,n,n) *Avec{p}; 
           end
        elseif strcmp(pblk{1},'q');       
           s = 1 + [0 cumsum(pblk{2})]; 
           gaptmp = mexqops(pblk,X{p},Z{p},1); 
           gamz2{p}  = mexqops(pblk,Z{p},Z{p},2) + 1e-14; 
           ee{p} = mexqops(pblk,-1./gamz2{p},Z{p},4);  
           dd{p} = mexqops(pblk,gaptmp./gamz2{p},ones(n,1),4); 
           tmp = dd{p}.*Rd{p} + mexqops(pblk,mexqops(pblk,Rd{p},ee{p},1),X{p},3) ...
                 + mexqops(pblk,mexqops(pblk,Rd{p},X{p},1),ee{p},3); 
           hRd = hRd + (tmp'*Avec{p})';
           EinvRc{p} = sigmu*ee{p} - X{p};  
           hEinvRc = hEinvRc + (EinvRc{p}'*Avec{p})';  
           Ax = mexqprod(pblk,Avec{p}',X{p}); 
           Ae = mexqprod(pblk,Avec{p}',ee{p}); 
           [decolidx,spcolidx] = checkdense(Ax);
           ddsch = dd{p};  
           if ~isempty(decolidx); 
              idx = s(decolidx); 
              tmp = zeros(n,1); 
              tmp(idx) = sqrt(2*abs(ddsch(idx))); 
              Ad = mexqprod(pblk,Avec{p}',tmp); 
              ddsch(idx) = abs(ddsch(idx)); 
              UU = [UU Ax(:,decolidx) Ae(:,decolidx)  Ad]; 
              VV = [VV Ae(:,decolidx) Ax(:,decolidx) -Ad]; 
              Ax = Ax(:,spcolidx); Ae = Ae(:,spcolidx); 
              schur = schur + Ax*Ae'+ Ae*Ax';
           else
              schur = schur + Ax*Ae'+ Ae*Ax';
           end
           [decolidx,spcolidx] = checkdense(Avec{p}'); 
           if ~isempty(decolidx); 
              len = length(decolidx);               
              tmp = (spdiags(sqrt(abs(ddsch(decolidx))),0,len,len)*Avec{p}(decolidx,:))'; 
              UU = [UU tmp]; 
              VV = [VV tmp*spdiags(sign(ddsch(decolidx)),0,len,len)]; 
              dd2 = zeros(n,1);  
              dd2(spcolidx) = ddsch(spcolidx); 
              schur = schur + Avec{p}' *spdiags(dd2,0,n,n) *Avec{p};
           else
              schur = schur + Avec{p}' *spdiags(ddsch,0,n,n) *Avec{p}; 
           end  
        elseif strcmp(pblk{1},'s')
           existsdpblk = 1; 
           numblk = length(pblk{2});  n2 = sum(pblk{2}.*pblk{2}); 
           Zinv{p} = Prod2(pblk,invZchol{p},invZchol{p}',1); 
           %%
           %% compute schur for matrices that are very sparse. 
           %% 
           if issparse(schur); schur = full(schur); end; 
           J = max(find(nzlistA{p,1} < inf)) -1; 
           if (J > 0)
              if issparse(X{p}) & ~issparse(Zinv{p}); X{p} = full(X{p}); end
              if ~issparse(X{p}) & issparse(Zinv{p}); Zinv{p} = full(Zinv{p}); end
              mexschur(pblk,Avec{p},nzlistA{p,1},nzlistA{p,2},permA(p,:),...
                       Zinv{p},X{p},J,0,schur);
           end
           %%
           %% compute schur for matrices that are not so sparse or dense.
           %% 
           L = max(find(nzlistAsum{p,1} < inf)) -1;  
           if (J < L)
              len = nzlistAsum{p,1}(J+1); list = nzlistAsum{p,2}(1:len,:); 
           end
           for k = J+1:m 
    	       isspAk = isspA(p,k);
               Ak = mexsmat(blk,Avec,isspAk,p,k); 
    	       if (k <= L)
                  idx1 = nzlistAsum{p,1}(k)+1; idx2 = nzlistAsum{p,1}(k+1);
                  list = [list; nzlistAsum{p,2}(idx1:idx2,:)]; 
                  list = sortrows(list,[2 1]); 
                  tmp = Prod2(pblk,Ak,Zinv{p}); 
                  tmp = mexProd2nz(pblk,X{p},tmp,list); 
               else
                  tmp = Prod3(pblk,X{p},Ak,Zinv{p},0);
               end
               tmp = 0.5*(mexsvec(pblk,tmp) + mexsvec(pblk,tmp,[],1));
               permk = permA(p,k);   
               idx = permA(p,1:k);  
               tmp2 = schur(idx,permk) + mexinprod(blk,Avec,tmp,k,p); 
               schur(idx,permk) = tmp2; 
               schur(permk,idx) = tmp2';
           end;
           clear tmp tmp2;
           tmp = Prod3(pblk,X{p},Rd{p},Zinv{p},0,nzlistAy{p}); 
           tmp = 0.5*(tmp+tmp'); 
           tmp2(permA(p,:),1) = mexinprod(blk,Avec,svec(pblk,tmp),m,p);
           hRd = hRd + tmp2;
           tmp = sigmu*Zinv{p} - X{p}; 
           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; 
%%
%% Cholesky factorization
%% 
    diagschur = diag(schur); 
    if (existsdpblk); 
       pertdiag = 1e-16*sqrt(m)*norm(diagschur);  
    else
       pertdiag = 1e-16*max(abs(diagschur));  
    end
    for k=1:m; 
        schur(k,k) = abs(schur(k,k)) + pertdiag; 
    end;  
    if (iter <= 3); 
       NNZSCHUR = nnz(schur);
    end;
    if (NNZSCHUR > spdensity*m^2); 
       if issparse(schur); schur = full(schur); end; 
       chol_options = 'matlab';
    else
       if ~issparse(schur), schur = sparse(schur); end;
       isdiagschur = (norm(schur-spdiags(diag(schur),0,m,m),inf) == 0);  
       if ~isdiagschur 
          chol_options = 'lipsol'; 
       else
          chol_options = 'matlab'; 
       end
    end
    if (chol_options == 'matlab');
       [schur,indef] = chol(schur);  
       if indef; solve_ok = -2; break; end;
       [dy,solve_ok] = linsysolve(schur,UU,VV,rhs,1); 
    elseif (chol_options == 'lipsol')
       if (iter == 1 | NNZSCHUR ~= nnzschurold);
          symbfct(schur);           
       end; 
       [dy,solve_ok] = linsysolve(schur,UU,VV,rhs,1); 
    end    
    nnzschurold = NNZSCHUR; 
%%
%% compute (dX,dZ)
%%  
    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      
%%*******************************************************************



