%                                      [d,vfrm] = updtransfo(x,z,w, dIN,K)
% UPDTRANSFO  Updated the Nesterov-Todd transformation using a
%  numerically stable method.
%
% SEE ALSO sedumi
% **********  INTERNAL FUNCTION OF SEDUMI **********
function  [d,vfrm] = updtransfo(x,z,w, dIN,K)
 %  
 %   This file is part of SeDuMi 1.05
 %   Copyright (C) 2001 Jos F. Sturm
 %     Dept. Econometrics & O.R., Tilburg University, the Netherlands.
 %     Supported by the Netherlands Organization for Scientific Research (NWO).
 %   Affiliation SeDuMi 1.03 and 1.04Beta (2000):
 %     Dept. Quantitative Economics, Maastricht University, the Netherlands.
 %   Affiliations up to SeDuMi 1.02 (AUG1998):
 %     CRL, McMaster University, Canada.
 %     Supported by the Netherlands Organization for Scientific Research (NWO).
 % 
 %   This program is free software; you can redistribute it and/or modify
 %   it under the terms of the GNU General Public License as published by
 %   the Free Software Foundation; either version 2 of the License, or
 %   (at your option) any later version.
 % 
 %   This program is distributed in the hope that it will be useful,
 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 %   GNU General Public License for more details.
 % 
 %   You should have received a copy of the GNU General Public License
 %   along with this program; if not, write to the Free Software
 %   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 %

% ------------------------------------------------------------
% PSD:
% Given w = D(Xscl)Zscl, compute spec-factor Q*WLAB*Q' = W
% ------------------------------------------------------------
 [wlab,q] = psdeig(w.s,K);
 w.lab(K.l+2*length(K.q) + 1:end) = wlab;
% ------------------------------------------------------------
% lambda(v) = sqrt(lambda(w))
% ------------------------------------------------------------
 vfrm.lab = sqrt(w.lab);
% ------------------------------------------------------------
% LP : d.l = dIN.l .* (x ./ z)
% ------------------------------------------------------------
 d.l = dIN.l .* (x(1:K.l) ./ z(1:K.l));
% ------------------------------------------------------------
% Lorentz:
% Auxiliary: s := sqrt(det(x)./det(z))
%   chi = (x[k] + s(k)*Jz[k]) / (lab1(v)+lab2(v))
%   psi = (x[k] - s(k)*Jz[k]) / (lab2(v)-lab1(v))
% Scale: detd = detdIN .* s and d = D(dIN)*chi
% ------------------------------------------------------------
 if length(K.q) == 0
   d.det = zeros(0,1); d.q1 = zeros(0,1); d.q2 = zeros(0,1);
   d.auxdet = zeros(0,1); d.auxtr = zeros(0,1);
   vfrm.q = zeros(0,1);
 else
   i1 = K.mainblks(1); i2 = K.mainblks(2); nq = i2 - i1; j3 = i2+nq-1;
   s = sqrt(w.tdetx ./ w.tdetz);
   d.det = dIN.det .* s;
   psi1 = s.*z(i1:i2-1); psi2 = qblkmul(s,z,K.qblkstart);       %s * z
   tmp = vfrm.lab(i1:i2-1) + vfrm.lab(i2:j3);
   chi1 = (x(i1:i2-1)+psi1)./tmp;
   chi2 = qblkmul(1./tmp, x(i2:K.lq)-psi2,K.qblkstart);
   psi1 = x(i1:i2-1)-psi1;
   psi2 = x(i2:K.lq)+psi2;
   dq = asmDxq(dIN,[chi1;chi2],K);                     % d = D(dIN)*chi
   d.q1 = dq(1:nq);
   d.q2 = dq(nq+1:end);
   d.auxdet = sqrt(2*d.det);
   d.auxtr = sqrt(2)*(d.q1 + d.auxdet);
   alpha = (dIN.q1 .* psi1 + ddot(dIN.q2,psi2,K.qblkstart)) ./ d.auxtr;
   tmp = 2*sqrt(s);
   psi1 = (psi1 - alpha .* chi1)./tmp;
   psi2 = psi2 - qblkmul(alpha,chi2,K.qblkstart);
   psi2 = qblkmul(1./tmp, psi2,K.qblkstart);
   gamma = (sqrt(2)*psi1+alpha) ./ dIN.auxtr;
   tmp = vfrm.lab(i2:j3) - vfrm.lab(i1:i2-1);
   tmp(tmp == 0) = 1;                                %avoid division by zero
   psi2 = psi2 + qblkmul(gamma,dIN.q2,K.qblkstart);
   vfrm.q = qblkmul(1./tmp, psi2, K.qblkstart);
 end
% ------------------------------------------------------------
% D= QUD'*QUD, where QUD(:,udIN.perm) = diag(1./sqrt(vlab))*Q'*ux*ud.
% Let vinv = diag(1./sqrt(vlab))*Q'.
% ------------------------------------------------------------'
 d.u = triumtriu(w.ux, dIN.u, K);             % ux * ud
 [d.u,d.perm,gjc,g] = urotorder(d.u,K, 1.1, dIN.perm);  % stable reordering
 q = givensrot(gjc,g,q,K);    % ROTATE Q accordingly: (G*Q)'*(G*ux*ud).'
 vinv = sqrtinv(q,vfrm.lab,K);
% ------------------------------------------------------------
% QR-FACTORIZE: Qv * VINV = R
% Then the new ud is simply R*ux*ud, which is upper triangular.
% ------------------------------------------------------------
 [vfrm.s, r] = qrK(vinv,K);
 d.u = triumtriu(r, d.u, K);
