%%**********************************************************************
%% NTpred: Compute (dX,dy,dZ) for NT direction. 
%%                       
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'. 
%% compute W satisfying W*Z*W = X. 
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%
%% 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,W,G,sv,dd,ee,rr,vv,gamx,gamz] = ...
          NTpred(blk,Avec,rp,Rd,sigmu,permA,isspA,isspAy,nzlistA,nzlistAsum,...
          nzlistAy,X,Z,Zchol,invZchol);

    global spdensity iter chol_options nnzschurold; 
    global NNZSCHUR; 
%%
%% 
%%
    m = length(rp); 
    schur = sparse(m,m); 
    UU = []; VV = [];
    hRd = zeros(m,1);  hEinvRc = zeros(m,1); 
    W = cell(size(blk,1),1);  G = cell(size(blk,1),1);   
    EinvRc = cell(size(blk,1),1);  sv = cell(size(blk,1),1);  
    existsdpblk = 0; 
%%
%% compute schur matrix
%%
    for p = 1:size(blk,1)
        pblk = blk(p,:); 
        numblk = length(pblk{2});  
        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})]; idx1 = s(1:numblk); 
           gaptmp = mexqops(pblk,X{p},Z{p},1); 
           gamx{p} = sqrt(mexqops(pblk,X{p},X{p},2)); 
           gamz2 = mexqops(pblk,Z{p},Z{p},2); 
           gamz{p} = sqrt(gamz2); 
           w2 = gamz{p}./gamx{p};  w = sqrt(w2); 
           dd{p} = mexqops(pblk,1./w2,ones(n,1),4);       
           tmp = mexqops(pblk,1./w,Z{p},3) - mexqops(pblk,w,X{p},4);
           ee{p} = mexqops(pblk,1./(w.*sqrt(gamx{p}.*gamz{p}+gaptmp)),tmp,4);
           rr{p} = mexqops(pblk,1./sqrt(2*(gamx{p}.*gamz{p}+gaptmp)),tmp,3); 
           inprod = mexqops(pblk,rr{p},X{p},1); 
           vvtmp = X{p} + mexqops(pblk,(X{p}(idx1)+inprod)./(1+rr{p}(idx1)),rr{p},3); 
           vvtmp(idx1) = inprod;       
           vv{p} = mexqops(pblk,w,vvtmp,3); 
           %%
           tmp = dd{p}.*Rd{p} + mexqops(pblk,mexqops(pblk,Rd{p},ee{p},1),ee{p},3); 
           hRd = hRd + (tmp'*Avec{p})';
           EinvRc{p} = mexqops(pblk,-sigmu./gamz2,Z{p},4)-X{p};
           hEinvRc = hEinvRc + (EinvRc{p}'*Avec{p})';         
           Ae = mexqprod(pblk,Avec{p}',ee{p}); 
           [decolidx,spcolidx] = checkdense(Ae);
           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 Ae(:,decolidx)  Ad]; 
              VV = [VV Ae(:,decolidx) -Ad]; 
              Ae = Ae(:,spcolidx); 
              schur = schur + Ae*Ae';
           else
              schur = schur + Ae*Ae';
           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}); 
           %%
           %% compute NT scaling matrix W = G'*G. 
           %% important to keep W{p} symmertic.
           %%
           tmp = Prod2(pblk,Zchol{p},X{p},0); 
           tmp = Prod2(pblk,tmp,Zchol{p}',1); 
           [sv2,V] = blkeig(pblk,tmp); 
           sv{p} = sqrt(sv2) + eps*max(sqrt(sv2));        
           tmp  = Prod2(pblk,invZchol{p},V); 
           clear V;
           G{p} = Prod2(pblk,spdiags(sqrt(sv{p}),0,n,n),tmp');
           clear tmp;
           W{p} = Prod2(pblk,G{p}',G{p},1);                    
           %%
           %% compute columns of 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)
              mexschur(pblk,Avec{p},nzlistA{p,1},nzlistA{p,2},permA(p,:),...
              W{p},W{p},J,1,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 = Prod3(pblk,W{p},Ak,W{p},1,list); 
              else
                 tmp = Prod3(pblk,W{p},Ak,W{p},1); 
              end
              permk = permA(p,k); 
              idx = permA(p,1:k);   
              tmp2 = schur(idx,permk) + mexinprod(blk,Avec,svec(pblk,tmp),k,p);
              schur(idx,permk) = tmp2; 
              schur(permk,idx) = tmp2'; 
           end
           clear tmp tmp2;
           tmp = Prod3(pblk,W{p},Rd{p},W{p},1,nzlistAy{p}); 
           tmp2(permA(p,:),1) = mexinprod(blk,Avec,svec(pblk,tmp),m,p);
           hRd = hRd + tmp2;  
           tmp = spdiags(sigmu./sv{p} -sv{p},0,n,n);
           tmp = Prod3(pblk,G{p}',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; 
%%
%% 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)
%%
    dX = cell(size(blk,1),1); dZ = cell(size(blk,1),1); 
    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      
%%**********************************************************************

