function [ainorm, y, it] = invnormest_simple(L,U)
%   ainorm = invnormest_simple(L,U)
%   given L,U from LU-decomposition estimate 1-norm of inverse matrix
%   Here L,U were obtained by [L,U,p] = lu(A,'vector') 
%
%   invnormest estimates the norm by solving (usually) 2 linear systems with A and 
%   2 linear systems with the matrix A' (using the given LU-decomposition for both).
% 
%   [ainorm,v] = invnormest(L,U)
%   returns the estimate ainorm, and a vector v such that 
%     norm(A*v,1) = norm(v,1)/ainorm
%   If ainorm is large, then v is an approximate null vector of A.
%   (Note: if a column permutation q was used, v is an approximate null vector of A(:,q) instead of A.)
%
%   [ainorm,v,it] = invnormest(L,U)
%   also returns the number of iterations. The algorithm solves it linear systems with the
%   matrix A and it linear systems with the matrix A'.

% 1-norm power method (by Boyd, Hager)
% see Higham, Tisseur: A Block Algorithm for Matrix 1-Norm Estimation

                  % let M = A(p,:)
n = size(L,1);
x = ones(n,1)/n;
for it=1:10       % usually only 2 iterations needed until termination
  y = U\(L\x);    % solve M*y = x
  xi = sign(y);
  z = L'\(U'\xi); % solve M'*z = xi
  if norm(z,Inf) <= z'*x
      ainorm = norm(y,1);
      break
  end
  [m,j] = max(abs(z));
  x = zeros(n,1); x(j) = 1;
end

