function [ p, delta ] = isl( p, nabla, n, r )
%ISL Applies the Iterated Scaling algorithm. 
%   The first input argument is the set of vectorized factor matrices.
%   The second input argument is the perturbation vector.
%   The third input argument is an array of integers containing the
%     dimensions of the tensor represented by the vectorized factor
%     matrices in p.
%   The fourth input argument is the number of terms in the rank
%     decomposition of the tensor represented by p.
%   The first output argument is the rescaled set of vectorized factor
%     matrices, representing the same tensor as p.
%   The second output argument is a vector perpendicular to the kernel
%     of Terracini's matrix, whose norm is approximately the square of the
%     norm of nabla.

% Constants
d = length(n);
cn = cumsum(n);
Sigma = cn(d);

% Construct kernel
K = zeros(r*Sigma, r*(d-1));
li = 1;
for i = 1 : r
   pt = p((i-1)*Sigma+1:i*Sigma);
   for k = 2 : d
       off = (i-1)*Sigma;
       K(off+1:off+cn(1),li) = pt(1:cn(1));
       K(off+cn(k-1)+1:off+cn(k),li) = -pt(cn(k-1)+1:cn(k));
       li = li + 1;
   end
end

% Iterated scaling algorithm
p0 = p;
pk = p;
delta = zeros(size(p));
gam = ones(d-1,r);
[U,S,V] = svd(K,0);
nablanorm = inf;
while nablanorm > 10*eps
    v = V*(S\(U'*nabla));
    gam = gam - reshape(v, [d-1 r]);
    zk = pk + K*v;
    off = 0;
    for i = 1 : r
        pk(off+1:off+cn(1)) = p0(off+1:off+cn(1))/prod(gam(:,i));
        off = off + n(1);
        for k = 2 : d
            pk(off+1:off+n(k)) = p0(off+1:off+n(k))*gam(k-1,i);
            off = off + n(k);
        end
    end
    nabla = zk - pk;
    delta = delta + (nabla - U*(U'*nabla));
    nabla = U*(U'*nabla);
    nablanorm = norm(nabla);
end
p = pk;
end

