function [N,flag,bound] = nulls(A,k,tol)
%NULLS  Null space of a (possibly sparse) matrix.
%
%  [N,flag,bound] = nulls(A,k,tol) attempts to compute the null space  of A. 
%
%  k      An initial estimate as to the dimension of the null space.
%         Defaults to 1 if ommitted. It is only a hint, which if
%         accurate, can speed up the computation (but not dramatically). 
%
%  tol    A tolerance: a vector v is considered to be a null
%         vector if norm(A*v) < tol*norm(A,1).
%         Defaults to 1e-14 if ommitted.
% 
%  N      A matrix whose columns are orthonormal and are null vectors of A.
%
%  flag   0 if the null space was computed successfully.
%         1 if the null space might be larger than the number of columns in N.
%
%  bound  If flag==0, this is the dimension of the null space of A.
%         If flag==1, this is an upper bound on the dimensio of the null space.
%
%  The algorithm uses normalized inverse U iteration with condition
%  estimation and normalized inverse L1*U iteration. For futher details,
%  see Gotsman and Toledo, "On the computation of null spaces of sparse
%  rectangular matrices (manuscript, 2005). 
% 
%  The main computational cost in this algorithm is the sparse LU 
%  factorization of A. It is suitable for any matrix that Matlab can
%  factor in sparse mode.
% 
%  Author: Sivan Toledo
%  Version 1.0, July 2005

  t0 = cputime;

  if ( issparse(A) == 0)
    A = sparse(A);
  end;

  [m n] = size(A);
  normA = norm(A,1);

  warning off MATLAB:nearlySingularMatrix

  if (nargin < 1)
    error('Not enough input arguments.');
  end

  if (nargin < 2) 
    k = 1;
  end

  if (nargin < 3)
    tol = 1e-14;
  end

  t1 = cputime;
    [L,U,P,COLPERM] = lu(A,1.0);		
    for j=1:n
      INVPERM(find(COLPERM(:,j))) = j;
    end
  t2 = cputime;
  fprintf(1,'Factorization time in nulls is %.2f seconds\n',t2-t1);
  fprintf(1,'nulls: nnz(L) = %.01e, nnz(U) = %.01e\n',nnz(L),nnz(U));

	% for the accuracy_u experiment
	%adu = abs(diag(U));
	%global acc_u;
	%global t;
	%[acc_u(t,2) acc_u(t,3) ] = min(adu);
	
  % factorization_error_norm = norm( A - P'*L*U*COLPERM' , 1) / norm( A, 1 )

  for i=1:n
    if (U(i,i)==0)
      % U(i,i)=eps*eps*normA;
      U(i,i)=1e-100;
    end
  end
  if (m > n)
    L1 = L(1:n,:);
  else
    L1 = L;
  end

  % factorization_error_norm = norm( A - P'*L*U*COLPERM' , 1) / norm( A, 1 )

  %%%
  %%% first, iterate to find the null space of U
  %%%

  kin  = k;
  kout = 0;
  while ( (kout < n) )
    % disp('U iteration')
    Nu = normalizedinversetriangulariteration(U,kin);
    [Nu,kout] = filter(A,Nu,tol*normA);
    if (kout < kin) 
      break;
    else
      kin = kin*2;
    end;
  end;

  dim_null_U = kout;

  %%%
  %%% now check whether L1 is ill conditioned
  %%%

  % disp('L1 iteration (for condition estimation)')
  Nl = normalizedinversetriangulariteration(L1,1);
  [Nl,kout] = filter(L1,Nl,tol*1.0);

  if (kout == 0)
    
    %%%
    %%% L1 is not ill conditioned, so we return the null space of U
    %%%

    N     = Nu;
    flag  = 0;
    bound = dim_null_U;

    fprintf(1,'nulls: case 1, null(A)==null(U), dim(null(A))==%d\n',bound);

  else 
    
    %%%
    %%% L1 is ill conditioned, so we run inverse iteration on L1*U
    %%%

    fprintf(1,'nulls: L1 is ill conditioned\n');

    kin  = dim_null_U + 1;
    kout = 0;
    while ( (kout < n) )
      % disp('LU iteration')
      Nlu = normalizedinverseLUiteration(L1,U,kin);
      [Nlu,kout] = filter2(L1,U,Nlu,tol*normA);
      if (kout < kin) 
        break;
      else
        kin = kin*2;
      end;
    end;

    bound = kout

    %%%
    %%% There are two special cases were we can still determine the null space
    %%%

    if (bound == dim_null_U)
      % We got lucky
      % L1 was ill conditioned, but the null space of A is that of U
      N     = Nu;
      flag  = 0;
      fprintf(1,'nulls: case 2, null(A)==null(U), dim(null(A))==%d\n',bound);
    elseif (bound == 1 & norm( A*Nlu(:,1) ) > tol*normA)
      % if L1*U has one null vector, and it's not a null vector
      % of A, then A has no null vectors
      N     = zeros(m,0);
      flag  = 0;
      bound = 0;
      fprintf(1,'nulls: case 3, null(A)==null(L1*U), dim(null(A))==0\n');
    else
      [N,R] = qr([ Nu Nlu ],0);
      [N,kout] = filter(A,N,tol*normA);
      kout_combined = kout

      if ( kout == bound )
        flag = 0;
        fprintf(1,'nulls: case 4, found null(A), dim(null(A))==%d\n',bound);
      else
        flag = 1;
        fprintf(1,'nulls: case 5, %d <= dim(null(A)) <= %d\n',kout,bound);
      end;

    end;
  end;

  t3 = cputime;
  fprintf(1,'Total time in nulls is %.2f seconds\n',t3-t0);

  warning on MATLAB:nearlySingularMatrix

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % normalized inverse LU iteration                    %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function W = normalizedinverseLUiteration(L,U,k)

    function X = invLU(W)
      Tt = U' \ W;
      Q = L' \ Tt;
      Y = L\Q;
      X = U \ Y;
    end;

    if (0) 

      opts.issymm = true;
      opts.isreal = true;

      [W,D] = eigs( @invLU, size(L,1), k, 'LM', opts);

    else % of if (1)

    W = rand(n,k);
    [W,R] = qr(W,0);

    for i=1:3;
      % iter = i
      % now we have an approximation W for the right singular vectors
      Tt = U' \ W;
      Q = L' \ Tt;
      [Q,R] = qr(Q,0);

      % we start with an approximation Q for the left singular vectors
      Y = L\Q;
      X = U \ Y;
      [W,R] = qr(X,0);
    end

    end; % of if (1)

    W = W(INVPERM,:);

    if (norm( isnan(W) + isinf(W) , 1 ) > 0)
      disp('nulls: ****************************************')
      disp('nulls: *** infs or nans in iverse iteration ***')
      disp('nulls: ****************************************')
      disp('nulls: *** infs or nans in iverse iteration ***')
    end;

  end; % of nested function

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % normalized inverse iteration                       %
  % for a triangular matrix Z                          %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function W = normalizedinversetriangulariteration(Z,k)

    function X = invZ(W)
      Q = Z' \ W;
      X = Z \ Q;
    end;

    if (0) 

      opts.issymm = true;
      opts.isreal = true;

      [W,D] = eigs( @invZ, size(Z,1), k, 'LM', opts);

    else % of if (1)

    W = rand(n,k);
    [W,R] = qr(W,0);

    for i=1:3;
      % now we have an approximation W for the right singular vectors
      Q = Z' \ W;
      [Q,R] = qr(Q,0);

      % we start with an approximation Q for the left singular vectors
      X = Z \ Q;
      [W,R] = qr(X,0);
    end

    end; % of if (1)

    W = W(INVPERM,:);

    if (norm( isnan(W) + isinf(W) , 1 ) > 0)
      disp('nulls: ****************************************')
      disp('nulls: *** infs or nans in iverse iteration ***')
      disp('nulls: ****************************************')
    end;

  end; % of nested function

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % filter only the prefix columns of W that satisfy   %
  % norm( Z*W(:,k) ) <= thresh                         %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function [Wout,kout] = filter(Z,W,thresh)

    [dummy k] = size(W);
    kout=0;
    while ((kout < k) & norm(Z*W(:,kout+1)) < thresh)
      kout = kout+1;
    end;
    Wout = W(:,1:kout);

  end; % of nested function

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % filter only the prefix columns of W that satisfy   %
  % norm( Z*T**W(:,k) ) <= thresh                      %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function [Wout,kout] = filter2(Z,T,W,thresh)

    [dummy k] = size(W);
    kout=0;
    while ((kout < k) & norm(Z*(T*W(:,kout+1))) < thresh)
      kout = kout+1;
    end;
    Wout = W(:,1:kout);

  end; % of nested function

end % of main function