# HG changeset patch
# User Xie Rui <calvocuno@gmail.com>
# Date 1488977767 -28800
#      Wed Mar 08 20:56:07 2017 +0800
# Node ID a0135f3ba6d0792ec2931195707ff9265be151a8
# Parent  4890b1c4a6bd53f409f994d327054ddd61a5e548
add minres for sparse linear systems

diff --git a/scripts/sparse/minres.m b/scripts/sparse/minres.m
new file mode 100644
--- /dev/null
+++ b/scripts/sparse/minres.m
@@ -0,0 +1,404 @@
+## -*- texinfo -*-
+## @deftypefn  {} {@var{x} =} minres (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0})
+## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} minres (@dots{})
+##
+## Solve the linear system of equations @w{@code{@var{A} * @var{x} = @var{b}}}
+## by means of the Minimum Residual Method.
+##
+## The input arguments are
+##
+## @itemize
+## @item
+## @var{A} should be a square, real and symmetric (preferably sparse) matrix 
+## which may be both indefinite and singular or a function
+## handle, inline function or string containing the name of a function which
+## computes @w{@code{@var{A} * @var{x}}}.
+##
+## @item
+## @var{b} is the right-hand side vector.
+##
+## @item
+## @var{tol} is the required relative tolerance for the residual error,
+## @w{@code{@var{b} - @var{A} * @var{x}}}.  
+## If @var{tol} is omitted or empty then a tolerance of 1e-6 is used.
+##
+## @item
+## @var{maxit} is the maximum allowable number of iterations; if @var{maxit}
+## is omitted or empty then a value of 100 is used.
+##
+## @item
+## @var{m} = @var{m1} * @var{m2} is the preconditioning matrix, so that
+## the iteration is (theoretically) equivalent to solving by @code{minres}
+## @w{@code{@var{P} * @var{x} = @var{m} \ @var{b}}}, with
+## @w{@code{@var{P} = @var{m} \ @var{A}}}.  
+## If @var{m1} is omitted or empty @code{[]} then no preconditioning is applied.
+## If @var{m2} is omitted, @var{m} = @var{m1} will be used as a preconditioner.
+##
+## @item
+## @var{x0} is the initial guess.  If @var{x0} is omitted or empty then the
+## function sets @var{x0} to a zero vector by default.
+## @end itemize
+##
+## The output arguments are
+##
+## @itemize
+## @item
+## @var{x} is the computed approximation to the solution of
+## @w{@code{@var{A} * @var{x} = @var{b}}}.
+##
+## @item
+## @var{flag} reports on the convergence.  A value of 0 means the solution
+## converged and the tolerance criterion given by @var{tol} is satisfied.
+## A value of 1 means that the @var{maxit} limit for the iteration count was
+## reached.
+##
+## @item
+## @var{relres} is the final relative residual,
+## measured in the Euclidean norm.
+##
+## @item
+## @var{iter} is the actual number of iterations performed.
+##
+## @item
+## @var{resvec(i)} is the Euclidean norm of the residual after the
+## (@var{i}-1)-th iteration, @code{@var{i} = 1, 2, @dots{}, @var{iter}+1}.
+##
+## @end itemize
+##
+##
+## Let us consider a trivial problem with a diagonal matrix (we exploit the
+## sparsity of A)
+##
+## @example
+## @group
+## n = 10;
+## A = diag (sparse (1:n));
+## b = rand (n, 1);
+## [l, u, p] = ilu (A, struct ("droptol", 1.e-3));
+## @end group
+## @end example
+##
+## @sc{Example 1:} Simplest use of @code{minres}
+##
+## @example
+## x = minres (A, b)
+## @end example
+##
+## @sc{Example 2:} @code{minres} with a function which computes
+## @code{@var{A} * @var{x}}
+##
+## @example
+## @group
+## function y = apply_a (x)
+##   y = [1: 10]' .* x;
+## endfunction
+##
+## x = minres ("apply_a", b)
+## @end group
+## @end example
+##
+## @sc{Example 3:} @code{minres} with a preconditioner: @var{l} * @var{u}
+##
+## @example
+## x = minres (A, b, 1.e-6, 500, l*u)
+## @end example
+##
+## @sc{Example 4:} @code{minres} with a preconditioner: @var{l} * @var{u}.
+## Faster than @sc{Example 3} since lower and upper triangular matrices are
+## easier to invert
+##
+## @example
+## x = minres (A, b, 1.e-6, 500, l, u)
+## @end example
+##
+## @sc{Example 5:} @code{minres} when @var{A} is indefinite. It fails with @code{pcg}.
+##
+## @example
+## A = diag([20:-1:1, -1:-1:-20]);
+## b = sum(A,2);
+## x = minres(A, b)
+## @end example
+##
+## References:
+##
+## @enumerate
+## @item
+## C. C. PAIGE and M. A. SAUNDERS, @cite{Solution of Sparse Indefinite Systems of Linear Equations},
+## SIAM J. Numer. Anal., 1975. (the minimum residual method)
+##
+## @end enumerate
+## @end deftypefn
+
+function [x, flag, relres, iter, resvec]  = minres(A, b, tol, maxit, m1, m2, x0, varargin)
+  ## Check the inputs
+  if (nargin < 2)
+      error('Not enough inputs!');
+  endif
+
+  [mb, nb] = size(b);
+  if (nb ~= 1)
+        error('b is not a vector!');
+  endif
+  
+  Aisnum = isnumeric(A);
+  if Aisnum
+    ## If A is a matrix
+    [ma, na] = size(A);
+    
+    if (ma ~= na)
+        error('A is not a square matrix!');
+    endif
+    if (na ~= mb)
+        error('A and b are not matched!');
+    endif
+    if ~isequal(A, A.')
+        error('A is not symmetric!');
+    endif
+  endif
+
+  if (nargin < 3) || isempty(tol)
+    tol = 1e-6;
+  endif
+    
+  if (nargin < 4) || isempty(maxit)
+    maxit = min(100, mb + 5);
+  endif
+
+  if (nargin >= 5) && ~isempty(m1)
+    ## Preconditioner exists
+    m1exist = true;
+    misnum = isnumeric(m1);
+      if (nargin >= 6) && ~isempty(m2)
+        m2exist = true;
+      else
+        m2exist = false;
+      endif
+  else
+    ## Preconditioner does not exist
+      m1exist = false;
+      m2exist = false;
+  endif
+
+  if (nargin >= 7) && ~isempty(x0)
+      [mx0, nx0] = size(x0);
+      if (mx0 ~= mb) || (nx0 ~= 1)
+          error('x0 has the wrong size!');
+      endif
+  else
+      x0 = zeros(mb, 1);
+  endif
+
+        
+  ## Preallocation
+  N = maxit;
+  flag = 1;
+
+  n = length(b);
+  beta = zeros(N, 1);
+  alpha = zeros(N, 1);
+  gamma = zeros(N, 1);
+  delta = zeros(N, 1);
+  epsilon = zeros(N, 1);
+  c = zeros(N, 1);
+  s = zeros(N, 1);
+  resvec = zeros(N, 1);
+  if Aisnum
+      resvec(1) = norm(A * x0 - b);
+  else
+      resvec(1) = norm(feval(A, x0, varargin{:}) - b);
+  endif
+
+  if (isequal(b, zeros(n, 1)))
+    ## If b is a zero vector
+    x = zeros(mb, 1);
+    flag = 0;
+    relres = NaN;
+    iter = 0;
+    resvec = [resvec(1); 0];
+    return
+  endif
+  
+  if m1exist
+    if m2exist
+      if misnum
+        by = m2 \ (m1 \ b);
+      else
+        by = feval(m2, feval(m1, b, varargin{:}), varargin{:});
+      endif
+    else
+      if misnum
+        by = m1 \ b;
+      else
+        by = feval(m1, b, varargin{:});
+      endif
+    endif
+  else
+    by = b;
+  endif
+
+  ## Initiation
+  v_0 = zeros(n, 1);
+  if Aisnum
+    if m1exist
+      if m2exist
+        if misnum
+          b0 = by - m2 \ (m1 \ (A * x0));
+        else
+          b0 = by - feval(m2, feval(m1, A * x0, varargin{:}), varargin{:});
+        endif
+      else
+        if misnum
+          b0 = by - m1 \ (A * x0);
+        else
+          b0 = by - feval(m1, A * x0, varargin{:});
+        endif
+      endif
+    else
+      b0 = by - A * x0;
+    endif
+  else
+    if m1exist
+      if m2exist
+        if misnum
+          b0 = by - m2 \ (m1 \ feval(A, x0, varargin{:}));
+        else
+          b0 = by - feval(m2, feval(m1, feval(A, x0, varargin{:})), varargin{:});
+        endif
+      else
+        if misnum
+          b0 = by - m1 \ feval(A, x0, varargin{:});
+        else
+          b0 = by - feval(m1, feval(A, x0, varargin{:}), varargin{:});
+        endif
+      endif
+    else
+      b0 = by - feval(A, x0, varargin{:});
+    endif
+  endif
+  beta(1) = norm(b0);
+  relres = beta(1) / norm(by);
+  if (relres <= tol) || (beta(1) <= eps)
+    ## If x0 is already good enouph
+    x = x0;
+    flag = 0;
+    iter = 0;
+    resvec = [resvec(1)];
+    return
+  endif
+  v_1 = b0 / beta(1);
+  v(:, 1) = v_1; 
+        
+  v_p = v_0;
+  v_n = v_1;
+  temp2 = beta(1);
+  m_p = zeros(n, 1);
+  m_pp = zeros(n, 1);
+  y = zeros(n, 1);
+
+  ## Iteration
+  for k = 1: (N - 1)
+    if Aisnum
+      if m1exist
+        if m2exist
+          if misnum
+            temp0 = m2 \ (m1 \ (A * (m2 \ (m1 \ v_n))));
+          else
+            temp0 = feval(m2, feval(m1, A * feval(m2, feval(m1, v_n, varargin{:}), varargin{:}), varargin{:}), varargin{:});
+          endif
+        else
+          if misnum
+            temp0 = m1 \ (A * (m1 \ v_n));
+          else
+            temp0 = feval(m1, A * feval(m1, v_n, varargin{:}), varargin{:});
+          endif
+        endif
+      else
+        temp0 = A * v_n;
+      endif
+    else
+      if m1exist
+        if m2exist
+          if misnum
+            temp0 = m2 \ (m1 \ feval(A, m2 \ (m1 \ v_n), varargin{:}));
+          else
+            temp0 = feval(m2, feval(m1, feval(A, feval(m2, feval(m1, v_n, varargin{:}), varargin{:}), varargin{:}), varargin{:}), varargin{:});
+          endif
+        else
+          if misnum
+            temp0 = m1 \ feval(A, m1 \ v_n, varargin{:});
+          else
+            temp0 = feval(m1, feval(A, feval(m1, v_n, varargin{:}), varargin{:}), varargin{:});
+          endif
+        endif
+      else
+        temp0 = feval(A, v_n, varargin{:});
+      endif
+    endif
+    alpha(k) = v_n' * temp0;
+    temp1 = temp0 - alpha(k) * v_n - beta(k) * v_p;
+    beta(k + 1) = norm(temp1);
+    v_f = temp1 / beta(k + 1);
+    v(:, k + 1) = v_f;
+
+    if k > 2
+       epsilon(k) = s(k - 2) * beta(k);
+       gamma_h = - c(k - 2) * beta(k) * s(k - 1) - alpha(k) * c(k - 1);
+       delta(k) = - c(k - 2) * beta(k) * c(k - 1) + alpha(k) * s(k - 1);
+    elseif k == 2
+       epsilon(k) = 0;
+       gamma_h = beta(2) * s(1) - alpha(2) * c(1);
+       delta(k) = beta(2) * c(1) + alpha(2) * s(1);
+    else
+       epsilon(k) = 0;
+       gamma_h = alpha(1);
+       delta(k) = 0;
+    endif
+
+    gamma(k) = sqrt(gamma_h ^ 2 + beta(k + 1) ^ 2);
+    c(k) = gamma_h / gamma(k);
+    s(k) = beta(k + 1) / gamma(k);
+
+    m = 1 / gamma(k) * (v_n - epsilon(k) * m_pp - delta(k) * m_p);
+
+    y = y + m * temp2 * c(k);
+    temp2 = temp2 * s(k);
+
+    if m1exist
+      if m2exist
+        if misnum
+          x = m2 \ (m1 \ y) + x0;
+        else
+          x = feval(m2, feval(m1, y, varargin{:}), varargin{:}) + x0;
+        endif
+      else
+        if misnum
+          x = m1 \ y + x0;
+        else
+          x = feval(m1, y, varargin{:}) + x0;
+        endif
+      endif
+    else
+      x = y + x0;
+    endif
+    if Aisnum
+       r = norm(A * x - b);
+    else
+       r = norm(feval(A, x, varargin{:}) - b);
+    endif
+    relres = r / norm(b);
+    resvec(k + 1) = r;
+    iter = k;
+    ## Check convergence
+    if (relres <= tol) || (beta(k + 1) <= eps)
+       flag = 0;
+       resvec = resvec(1: (k + 1));
+       break
+    endif
+
+    m_pp = m_p;
+    m_p = m;
+    v_p = v_n;
+    v_n = v_f;
+
+  endfor 
+endfunction
\ No newline at end of file
