## Copyright (C) 2016 Marco Caliari
##
## This file is part of Octave.
##
## Octave 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 3 of the License, or (at
## your option) any later version.
##
## Octave 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 Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{c} =} normest1 (@var{a})
## @deftypefnx {Function File} {@var{c} =} normest1 (@var{a}, @var{t})
## @deftypefnx {Function File} {@var{c} =} normest1 (@var{a}, @var{t}, @var{x0})
## @deftypefnx {Function File} {@var{c} =} normest1 (@var{afun}, @var{t}, @var{x0}, @var{p1}, @var{p2}, @dots{})
## @deftypefnx  {Function File} {[@var{c}, @var{v}] =} normest1 (@var{a}, @dots{})
## @deftypefnx  {Function File} {[@var{c}, @var{v}, @var{w}] =} normest1 (@var{a}, @dots{})
## @deftypefnx  {Function File} {[@var{c}, @var{v}, @var{w}, @var{it}] =} normest1 (@var{a}, @dots{})
##
## Estimate the 1-norm of the matrix @var{a} using a block algorithm.
##
## For a medium size matrix @var{a}, @code{norm (@var{a}, 1)} should be
## used instead. For a large sparse matrix, when only an estimate of the norm
## is needed, @code{normest1 (@var{a})} might be faster. Moreover, it can be
## used for the estimate of the 1-norm of a linear
## operator @var{a} when matrix-vector products @code{@var{a} * @var{x}} and
## @code{@var{a}' * @var{x}} can be cheaply computed. In this case,
## instead of the matrix @var{a}, a function
## @code{@var{afun} (@var{flag}, @var{x})} can be used. It should returns
## @itemize @bullet
## @item
## the dimension @var{n} of @var{a}, if @var{flag} is "dim"
## @item
## true if @var{a} is a real operator, if @var{flag} is "real"
## @item
## the result @code{@var{a} * @var{x}}, if @var{flag} is "notransp"
## @item
## the result @code{@var{a}' * @var{x}}, if @var{flag} is "transp"
## @end itemize
## A typical case is @var{a} defined by @code{@var{b} ^ @var{m}},
## in which the result @code{@var{a} * @var{x}} can be computed without
## even forming explicitely @code{@var{b} ^ @var{m}} by
##
## @code{@var{y} = @var{x};, for @var{i} = 1:@var{m}, @var{y} = @var{b} * @var{y}; endfor}
##
## The parameters @var{p1}, @var{p2}, {@dots} are arguments of
## @code{@var{afun} (@var{flag}, @var{x}, @var{p1}, @var{p2}, @dots)}.
##
## The default value for @var{t} is 2. The algorithm requires
## matrix-matrix products with sizes @var{n} x @var{n} and
## @var{n} x @var{t}.
##
## The initial matrix @var{x0} should have columns of unit 1-norm.
## The default initial matrix @var{x0} has the first column
## @code{ones (@var{n}, 1) / @var{n}}
## and, if @var{t} >  1, the remaining columns with random elements
## @code{-1 / @var{n}}, @code{1 / @var{n}}, divided by @var{n}.
## Therefore, if consistent results are required, the "state" of the
## random generator should be fixed before invoking @code{normest1}.
##
## On output, @var{c} is the desired estimate, @var{v} and @var{w}
## vectors such that @code{@var{w} = @var{a} * @var{v}}, with
## @code{norm (@var{w}, 1)} = @code{@var{c} * norm (@var{v}, 1)}.
## @var{it} contains in @code{@var{it}(1)} the number of iterations
## (the maximum number is hardcoded to 5) and in  @code{@var{it}(2)}
## the total number of products @code{@var{a} * @var{x}} or
## @code{@var{a}' * @var{x}} performed by the algorithm.
##
## Reference: @nospell{N. J. Higham and F. Tisseur},
## @cite{A block algorithm for matrix 1-norm estimation, with and
## application to 1-norm pseudospectra}, SIAM J. Matrix Anal. Appl.,
## pp. 1185--1201, Vol 21, No. 4, 2000.
##
## @seealso{normest, rand}
## @end deftypefn

function [est, v, w, k] = normest1(A, varargin)

  if ((nargin >= 3) && ! isnumeric (A) && ! isnumeric (varargin{1}))
     # onenormest syntax, deprecated in 4.2
     warning("Octave:deprecated-option",
             "This syntax is deprecated, call normest1 (AFUN, T, X0, P1, P2, ...) instead.");
     [est, v, w, k] = onenormest (A, varargin{:});
     return
  endif

  if (nargin < 2 || isempty (varargin{1}))
    t = 2;
  else
    t = varargin{1};
  endif

 ## FIXME: t < 0 should print trace information
  if isnumeric (A)
    ## A is a matrix
    Aisnum = true;
    n = length (A);
    if ((n <= 4) || (t == n))
    ## compute directly
      [est, idx] = max (sum (abs (A), 1), [] ,2);
      v = zeros(n, 1);
      v(idx) = 1;
      w = A(:, idx);
      k = [0, 1];
      return
    else
      realm = isreal (A);
    endif
  else
    Aisnum = false;
    n = A ("dim", [], varargin{3:end});
    realm = A ("real", [], varargin{3:end});
    Afun = @(x) A ("notransp", x, varargin{3:end});
    A1fun = @(x) A ("transp", x, varargin{3:end});
  endif

  t = min (t, n);

  if ((nargin <= 2) || isempty (varargin{2}))
    X = [ones(n, 1), sign(2 * rand(n, t - 1) - 1)];
    i = 2;
    imax = min (t, 2 ^ (n - 1));
    ##  There are at most 2^(n-1) unparallel columns, see later.
    while (i <= imax)
      if (any (abs (X(:,i)' * X(:,1:i-1)) == n))
        ## column i is parallel to a colum 1:i-1. Change it.
        X(:,i) = sign (2 * rand (n, 1) - 1);
      else
        i++;
      endif
    endwhile
    X /= n;
  else
    X = varargin{2};
    if (columns (X) < t)
       error (sprintf ("X0 should have %d columns", t));
    endif
  endif

  itmax = 5;
  ind_hist = zeros (n, 1);
  est_old = 0;
  ind = zeros (n, 1);
  S = zeros (n, t);
  k = [0; 0];
  conv = false;
  while ((! conv) && (k(1) < itmax))
    k(1)++;
    if (Aisnum)
      Y = A * X;
    else
      Y = Afun (X);
    endif
    k(2)++;
    [est, j] = max (sum (abs (Y), 1), [], 2);
    if ((est > est_old) || (k(1) == 2))
      ind_best = ind(j);
      w = Y(:, j); ## there is an error in Algorithm 2.4
    endif
    if ((est <= est_old) && (k(1) >= 2)) ## (1) of Algorithm 2.4
      est = est_old;
      break; ## while
    endif
    est_old = est;
    Sold = S;
    S = sign (Y);
    S(S==0) = 1;
    possible_break = false;
    if (realm)
      ## test parallel (only real case)
      if (all (any (abs (Sold' * S) == n))) ## (2) of Algorithm 2.4
        ## all columns of S parallel to a column of Sold, exit
        possible_break = true;
        conv = true;
      else
        if (t > 1)
          ## at least two columns of S are not parallel
          i = 1;
          imax = min (t, 2 ^ (n - 1));
          while (i <= imax)
            ## The maximum number of parallel columns of length n with entries
            ## {-1,1} is 2^(n-1). Therefore, if the number of columns of S is
            ## greater than 2^(n-1), for sure some of them are parallel to some
            ## columns of Sold. Don't even try to change them (i <= 2^(n-1)).
            ## Now, check if S(:,i) is parallel to any previous column of S
            p = (any (abs (S(:,i)' * S(:,1:i-1)) == n));
            if (p || (any (abs (S(:,i)' * Sold) == n)))
              ## i-th column of S parallel to a previous or to a
              ## column of Sold: change it.
              S(:,i) = sign (2*rand (n, 1)-1);
            else
              i++;
            endif
          endwhile
        endif
      endif
    endif
    if (! possible_break)
      if (Aisnum)
        Z = A' * S;
      else
        Z = A1fun (S); ## (3) of Algorithm 2.4
      endif
      k(2)++;
      h = max (abs (Z), [], 2);
      ind = (1:n)';
      if ((k(1) >= 2) && (max (h, [], 1) == h(ind_best))) ## (4) of Alg. 2.4
        break; ## while
      endif
      [h, ind] = sort (h, "descend");
      if (t > 1)
       if (all (ind_hist(ind(1:t)))) ## (5) of Algorithm 2.4
          break ## while
        endif
        ind = ind(! ind_hist(ind));
        ## length(ind) could be less than t, especially if t is not << n.
        ## This is not considered in point (5) of Algorithm 2.4.
        tmax = min (length (ind), t);
        ind = ind(1:tmax);
      else
        tmax = 1;
      endif
      X = zeros (n,tmax);
      X(sub2ind (size (X), ind(1:tmax), (1:tmax)')) = 1;
      ind_hist(ind(1:tmax)) = 1;
    endif
  endwhile
  v = zeros (n, 1);
  v(ind_best) = 1;
endfunction

function [est, v, w, iter] = onenormest (varargin) # to be removed
                                # after 4.2

  if (nargin < 1 || nargin > 4)
    print_usage ();
  endif

  default_t = 5;
  itmax = 10;

  if (isnumeric (varargin{1}))
    [n, nc] = size (varargin{1});
    if (n != nc)
      error ("onenormest: matrix must be square");
    endif
    apply = @(x) varargin{1} * x;
    apply_t = @(x) varargin{1}' * x;
    if (nargin > 1)
      t = varargin{2};
    else
      t = min (n, default_t);
    endif
    issing = isa (varargin{1}, "single");
  else
    if (nargin < 3)
      print_usage ();
    endif
    apply = varargin{1};
    apply_t = varargin{2};
    n = varargin{3};
    if (nargin > 3)
      t = varargin{4};
    else
      t = default_t;
    endif
    issing = isa (n, "single");
  endif

  ## Initial test vectors X.
  X = rand (n, t);
  X ./= ones (n,1) * sum (abs (X), 1);

  ## Track if a vertex has been visited.
  been_there = zeros (n, 1);

  ## To check if the estimate has increased.
  est_old = 0;

  ## Normalized vector of signs.
  S = zeros (n, t);

  if (issing)
    myeps = eps ("single");
    X = single (X);
  else
    myeps = eps;
  endif

  for iter = 1 : itmax + 1
    Y = feval (apply, X);

    ## Find the initial estimate as the largest A*x.
    [est, ind_best] = max (sum (abs (Y), 1));
    if (est > est_old || iter == 2)
      w = Y(:,ind_best);
    endif
    if (iter >= 2 && est < est_old)
      ## No improvement, so stop.
      est = est_old;
      break;
    endif

    est_old = est;
    S_old = S;
    if (iter > itmax),
      ## Gone too far.  Stop.
      break;
    endif

    S = sign (Y);

    ## Test if any of S are approximately parallel to previous S
    ## vectors or current S vectors.  If everything is parallel,
    ## stop.  Otherwise, replace any parallel vectors with
    ## rand{-1,+1}.
    partest = any (abs (S_old' * S - n) < 4*eps*n);
    if (all (partest))
      ## All the current vectors are parallel to old vectors.
      ## We've hit a cycle, so stop.
      break;
    endif
    if (any (partest))
      ## Some vectors are parallel to old ones and are cycling,
      ## but not all of them.  Replace the parallel vectors with
      ## rand{-1,+1}.
      numpar = sum (partest);
      replacements = 2*(rand (n,numpar) < 0.5) - 1;
      S(:,partest) = replacements;
    endif
    ## Now test for parallel vectors within S.
    partest = any ((S' * S - eye (t)) == n);
    if (any (partest))
      numpar = sum (partest);
      replacements = 2*(rand (n,numpar) < 0.5) - 1;
      S(:,partest) = replacements;
    endif

    Z = feval (apply_t, S);

    ## Now find the largest non-previously-visted index per vector.
    h = max (abs (Z),2);
    [mh, mhi] = max (h);
    if (iter >= 2 && mhi == ind_best)
      ## Hit a cycle, stop.
      break;
    endif
    [h, ind] = sort (h, 'descend');
    if (t > 1)
      firstind = ind(1:t);
      if (all (been_there(firstind)))
        ## Visited all these before, so stop.
        break;
      endif
      ind = ind(! been_there(ind));
      if (length (ind) < t)
        ## There aren't enough new vectors, so we're practically
        ## in a cycle.  Stop.
        break;
      endif
    endif

    ## Visit the new indices.
    X = zeros (n, t);
    for zz = 1 : t
      X(ind(zz),zz) = 1;
    endfor
    been_there(ind(1 : t)) = 1;
  endfor

  ## The estimate est and vector w are set in the loop above.
  ## The vector v selects the ind_best column of A.
  v = zeros (n, 1);
  v(ind_best) = 1;

endfunction

%!function z = afun_A (flag, x, A, n)
%! switch flag
%! case {"dim"}
%!   z = n;
%! case {"real"}
%!   z = isreal (A);
%! case {"transp"}
%!   z = A' * x;
%! case {"notransp"}
%!   z = A * x;
%! endswitch
%!endfunction
%!function z = afun_A_P (flag, x, A, m)
%! switch flag
%! case "dim"
%!   z = length (A);
%! case "real"
%!   z = isreal (A);
%! case "transp"
%!   z = x; for i = 1:m, z = A' * z;, endfor
%! case "notransp"
%!   z = x; for i = 1:m, z = A * z;, endfor
%! endswitch
%!endfunction
%!test
%! ## simple test
%! A = reshape ((1:16)-8, 4, 4);
%! assert (normest1 (A), norm (A, 1), eps)
%!test
%! ## test t=1
%! A = rand (4); % for positive matrices always work
%! assert (normest1(A, 1), norm(A,1), 2 * eps)
%!test
%! ## test t=3
%! A = [-0.21825   0.16598   0.19388   0.75297;...
%!      -1.47732   0.78443  -1.04254   0.42240;...
%!       1.39857  -0.34046   2.28617   0.68089;...
%!       0.31205   1.50529  -0.75804  -1.22476];
%! X = [1,1,-1;1,1,1;1,1,-1;1,-1,-1]/3;
%! assert (normest1 (A, 3, X), norm (A, 1), 2 * eps)
%!test
%! ## test afun
%! A = rand (10);
%! n = length (A);
%! afun = @(flag, x) afun_A (flag, x, A, n);
%! assert (normest1 (afun), norm (A, 1), 2 * eps)
%!test
%! ## test afun with parameters
%! A = rand (10);
%! assert (normest1 (@afun_A_P, [], [], A, 3), norm (A ^ 3, 1), 2 * eps)
%!test
%! ## test output
%! A = reshape(1:16,4,4);
%! [c, v, w, it] = normest1 (A);
%! assert(norm (w, 1), c * norm(v, 1), eps);
%! ## test output
%!test
%! A = rand (100);
%! A(A <= 1/3) = -1;
%! A(A > 1/3 & A <= 2/3) = 0;
%! A(A > 2/3) = 1;
%! [c, v, w, it] = normest1 (A, 10);
%! assert(w, A * v, eps);
%!test
%! A = rand (5);
%! c = normest1 (A, 6);
%! assert (c, norm (A,1), eps);
%!test
%! A = rand (5);
%! c = normest1 (A, 2, ones (5,2) / 5);
%! assert (c, norm (A,1), eps);

%!demo # to be removed after 4.2
%! N = 100;
%! A = randn (N) + eye (N);
%! [L,U,P] = lu (A);
%! nm1inv = normest1 (@(x) U\(L\(P*x)), @(x) P'*(L'\(U'\x)), N, 30)
%! norm (inv (A), 1)

%!test
%! N = 10;
%! A = ones (N);
%! [nm1, v1, w1] = normest1 (A);
%! [nminf, vinf, winf] = normest1 (A', 6);
%! assert (nm1, N, -2*eps);
%! assert (nminf, N, -2*eps);
%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);

%!test # to be removed after 4.2
%! N = 10;
%! A = ones (N);
%! [nm1, v1, w1] = normest1 (@(x) A*x, @(x) A'*x, N, 3);
%! [nminf, vinf, winf] = normest1 (@(x) A'*x, @(x) A*x, N, 3);
%! assert (nm1, N, -2*eps);
%! assert (nminf, N, -2*eps);
%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);

%!test
%! N = 5;
%! A = hilb (N);
%! [nm1, v1, w1] = normest1 (A);
%! [nminf, vinf, winf] = normest1 (A', 6);
%! assert (nm1, norm (A, 1), -2*eps);
%! assert (nminf, norm (A, inf), -2*eps);
%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);

## Only likely to be within a factor of 10.
%!test
%! old_state = rand ("state");
%! restore_state = onCleanup (@() rand ("state", old_state));
%! rand ("state", 42);  # Initialize to guarantee reproducible results
%! N = 100;
%! A = rand (N);
%! [nm1, v1, w1] = normest1 (A);
%! [nminf, vinf, winf] = normest1 (A', 6);
%! assert (nm1, norm (A, 1), -.1);
%! assert (nminf, norm (A, inf), -.1);
%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
