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') 
%     For sparse A use
%     [L,U,p,q] = lu(A,'vector');       % A must be in sparse format
%
%   If A is symmetric positive definite: 
%     C = chol(A) gives Cholesky decomposition A = C'*C, then use
%     ainorm = invnormest(C',C)   
%     For sparse symm. pos. def. A use
%     [C,flag,p] = chol(A,'vector');    % A must be in sparse format
%
%   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)
% improved algorithm:  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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IDEA of ALGORITHM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Let B=inv(M), we want to estimate norm(B,1). We actually don't know the components of B,
% but we can easily compute matrix-vector products B*v and B'*v using the LU-decomposition.
%
% We know that norm(B,1) can be obtained as norm(B*x,1) where x=[0,...,0,1,0,...,0]
% with x_j=1 and j is column of B with largest column sum.
%
% The problem is that we don't know j (and trying all j=1...n actually means computing inv(M)).
%
% We can obtain the column sums of absolute values of B by computing z=B'*xi where each component 
% of xi is either 1 or -1, and the signs are chosen according to signs in column j of B. 
% Then j is the index where |z(j)| is maximal.
%
% The problem is that we don't know the correct signs in column j.
% Assume first that elements in column j have huge absolute values, and the other columns don't.
% Then we can find the signs by looking at the average of each row:
%
% Start with x = [1,...,1]/n and find y=B*x which gives the average of each row of B.
% Then pick the vector xi according to the signs of the vector xi.
% Then compute z=B'*xi which hopefully gives the column sums of absolute values of B.
% Let j be the component of z with the largest absolute value, 
% let x be the unit vector with x(j)=1, x(k)=0 for all other components.
% Then compute y=B*x and let ainorm:=norm(y,1).
%
% If j is really the column with largest sum of absolute values:
% y=B*x is this column, and xi:=sign(y) are signs of this column.
% Then z:=B'*xi gives column sums of B with certain signs, but z'*x=z(j)=norm(B,1) should be
% component with largest absolute value. Hence we should have norm(z,Inf)=z'*x.
% If norm(z,Inf) > z(j) this means that we can actually get a larger column sum in some other 
% column k of the matrix B. Therefore we should then let x be the vector with x(k)=1 and 0 otherwise.
% Here k is the index of component of z with largest absolute value.
%
% This algorithm is called "1-norm power method". For more information see the explanation of 
% the "p-norm power method" (Algorithm PM) in 
%    N. Higham, "Estimating the matrix p-norm", Numerische Mathematik 1991
%    citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.34.9868

