function [L, p] = cholincsp(X, tol)
%
% CHOLINCSP - Incomplete Cholesky factorization with symmetric pivoting
%
%    R = CHOLINCSP(X) computes an inclomplete Cholesky factorization of a
%    symmetrix matrix X, such that R*R' is a permutation of the rows and
%    columns of X.  The columns of Cholesky factor R are computed one at
%    a time, with symmetric pivoting, as described in [1]. 
%
%    [R,I] = CHOLINC(X) also gives the permutation in I, such that
%    R*R' = X(I,I) .  If R contains n columns, then the first n elements of I
%    specify a set of linearly independent columns forming a basis for the
%    remaining columns of X.
%
%    R = CHOLINCSP(X, TOL) allows the termination criterion (TOL) to be
%    adjusted (a default value of 1e-12 is used if only one argument is
%    given.
%
%    References :
%
%       [1] Fine, S. and Scheinberg, K., "Efficient SVM training using
%           low-rank kernel representations", Journal of Machine Learning
%           Research, vol. 2, pp 243-264, December 2001.
%
%    See also : CHOLINCSP_TEST

%
% File        : cholincsp.m
%
% Date        : Wednesday 11th August 2004
%
% Author      : Gavin C. Cawley (with help from Katya Scheinberg)
%
% Description : Vectorised MATLAB implementation of Fine and Scheinberg's
%               incomplete Cholesky factorisation with symmetric pivoting [1].
%
% References  : [1] Fine, S. and Scheinberg, K., "Efficient SVM training
%                   using low-rank kernel representations", Journal of Machine
%                   Learning Research, vol. 2, pp 243-264, December 2001.
%
% History     : 10/08/2004 - v1.00
%
% To do       : Write an even faster MEX version using BLAS etc!
%
% Copyright   : (c) Dr Gavin C. Cawley, August 2004.
%
%    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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

if nargin == 1

   tol = 1e-12;

end

m         = size(X,1);
p         = 1:m;
trc       = zeros(1,m);
L         = zeros(m, m);
L(:,1)    = X(:,1);
dL        = zeros(1, m);
diagonal  = sub2ind([m,m], 1:m, 1:m); 

for i=1:m

   L(diagonal(i:m)) = X(diagonal(p(i:m))) - sum(L(i:m,1:i-1).^2,2)';

   dL = diag(L);

   trc(i)=sum(dL(i:m));

   if  trc(i) <= tol

     break

   else

      [v q] = max(dL(i:m));

      q = q+i-1;

      t    = p(i);
      p(i) = p(q);
      p(q) = t;

      tt=L(i, 1:i-1);
      L(i, 1:i-1)=L(q, 1:i-1);
      L(q, 1:i-1)=tt;

      v = sqrt(v);

      L(i,i)=v; 

      L(i+1:m,i) = (X(p(i), p(i+1:m))' - L(i+1:m,1:i-1)*L(i,1:i-1)')/v;

   end

end

if i ~= m

   L = L(:,1:i-1);

end

% bye bye...

