# HG changeset patch
# User Cristiano Dorigo <cristiano.dorigo@hotmail.it>
# Date 1461936174 -7200
#      Fri Apr 29 15:22:54 2016 +0200
# Node ID 8208b851edb6fcbcc780bc88bfbe0aa2da955dbd
# Parent  7bdfa98a98a30c5d71a3f74f057f091586f9b4b8
Add external function to control default input for gmres.m and pcg.m

* scripts/sparse/private/__default__input__.m: new file
* scripts/sparse/gmres.m: now uses __default__input.m to control the settings of default inputs
* scripts/sparse/pcg.m:  now uses __default__input.m to control the settings of default inputs
* scripts/sparse/cgs.m:  now uses __default__input.m to control the settings of default inputs
* scripts/sparse/bicg.m:  now uses __default__input.m to control the settings of default inputs
* scripts/sparse/bicgstab.m:  now uses __default__input.m to control the settings of default inputs
* scripts/sparse/pcr.m:  now uses __default__input.m to control the settings of default inputs
* scripts/sparse/qmr.m:  now uses __default__input.m to control the settings of default inputs

diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/bicg.m
--- a/scripts/sparse/bicg.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/bicg.m	Fri Apr 29 15:22:54 2016 +0200
@@ -73,7 +73,7 @@
 ## Author: Sylvain Pelissier <sylvain.pelissier@gmail.com>
 ## Author: Carlo de Falco
 
-function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0)
+function [x, flag, res1, k, resvec] = bicg (A, b, varargin) 
 
   if (nargin >= 2 && isvector (full (b)))
 
@@ -90,18 +90,17 @@ function [x, flag, res1, k, resvec] = bi
     else
       error ("bicg: A must be a square matrix or function");
     endif
-
-    if (nargin < 3 || isempty (tol))
-      tol = 1e-6;
+    
+    if (length (varargin) < 5)
+      varargin {5} = [];
     endif
 
-    if (nargin < 4 || isempty (maxit))
-      maxit = min (rows (b), 20);
-    else
-      maxit = fix (maxit);
-    endif
+    [tol, maxit, x0] = __default__input__ ({1e-6, min(rows (b), 20), zeros( size(b))}, varargin {1}, fix (varargin {2}), varargin {5});
 
-    if (nargin < 5 || isempty (M1))
+    M1 = varargin {3};
+    M2 = varargin {4};
+
+    if (isempty (M1))
       M1m1x = @(x, ignore) x;
       M1tm1x = M1m1x;
     elseif (ischar (M1))
@@ -118,7 +117,7 @@ function [x, flag, res1, k, resvec] = bi
       error ("bicg: preconditioner M1 must be a function or matrix");
     endif
 
-    if (nargin < 6 || isempty (M2))
+    if (isempty (M2))
       M2m1x = @(x, ignore) x;
       M2tm1x = M2m1x;
     elseif (ischar (M2))
@@ -138,10 +137,6 @@ function [x, flag, res1, k, resvec] = bi
     Pm1x  = @(x) M2m1x  (M1m1x (x));
     Ptm1x = @(x) M1tm1x (M2tm1x (x));
 
-    if (nargin < 7 || isempty (x0))
-      x0 = zeros (size (b));
-    endif
-
     y = x = x0;
     c = b;
 
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/bicgstab.m
--- a/scripts/sparse/bicgstab.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/bicgstab.m	Fri Apr 29 15:22:54 2016 +0200
@@ -69,8 +69,7 @@
 ##
 ## @end deftypefn
 
-function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit,
-                                                     M1, M2, x0)
+function [x, flag, relres, iter, resvec] = bicgstab (A, b, varargin)
 
   if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
 
@@ -83,16 +82,17 @@ function [x, flag, relres, iter, resvec]
     else
       error ("bicgstab: A must be a square matrix or function");
     endif
+    
+    if (length (varargin) < 5)
+      varargin {5} = [];
+    endif
+    
+    [tol, maxit, x] = __default__input__ ({1e-6, min(rows (b),20), zeros(rows (b), 1)}, varargin {1}, varargin {2}, varargin {5});
 
-    if (nargin < 3 || isempty (tol))
-      tol = 1e-6;
-    endif
-
-    if (nargin < 4 || isempty (maxit))
-      maxit = min (rows (b), 20);
-    endif
-
-    if (nargin < 5 || isempty (M1))
+    M1 = varargin {3};
+    M2 = varargin {4};
+   
+    if (isempty (M1))
       M1m1x = @(x) x;
     elseif (ischar (M1))
       M1m1x = str2func (M1);
@@ -104,7 +104,7 @@ function [x, flag, relres, iter, resvec]
       error ("bicgstab: preconditioner M1 must be a function or matrix");
     endif
 
-    if (nargin < 6 || isempty (M2))
+    if (isempty (M2))
       M2m1x = @(x) x;
     elseif (ischar (M2))
       M2m1x = str2func (M2);
@@ -118,17 +118,6 @@ function [x, flag, relres, iter, resvec]
 
     precon = @(x) M2m1x (M1m1x (x));
 
-    if (nargin < 7 || isempty (x0))
-      x0 = zeros (size (b));
-    endif
-
-    ## specifies initial estimate x0
-    if (nargin < 7)
-      x = zeros (rows (b), 1);
-    else
-      x = x0;
-    endif
-
     norm_b = norm (b);
 
     res = b - Ax (x);
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/cgs.m
--- a/scripts/sparse/cgs.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/cgs.m	Fri Apr 29 15:22:54 2016 +0200
@@ -68,7 +68,7 @@
 ## @seealso{pcg, bicgstab, bicg, gmres, qmr}
 ## @end deftypefn
 
-function [x, flag, relres, iter, resvec] = cgs (A, b, tol, maxit, M1, M2, x0)
+function [x, flag, relres, iter, resvec] = cgs (A, b, varargin) 
 
   if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
 
@@ -81,16 +81,17 @@ function [x, flag, relres, iter, resvec]
     else
       error ("cgs: A must be a square matrix or function");
     endif
+    
+    if (length (varargin) < 5)
+      varargin {5} = [];
+    endif
+    
+    [tol, maxit, x0] = __default__input__ ({1e-6, min(rows (b), 20), zeros(size (b))}, varargin {1}, varargin {2}, varargin {5});
+    
+    M1 = varargin {3};
+    M2 = varargin {4};
 
-    if (nargin < 3 || isempty (tol))
-      tol = 1e-6;
-    endif
-
-    if (nargin < 4 || isempty (maxit))
-      maxit = min (rows (b), 20);
-    endif
-
-    if (nargin < 5 || isempty (M1))
+    if (isempty (M1))
       M1m1x = @(x) x;
     elseif (ischar (M1))
       M1m1x = str2func (M1);
@@ -102,7 +103,7 @@ function [x, flag, relres, iter, resvec]
       error ("cgs: preconditioner M1 must be a function or matrix");
     endif
 
-    if (nargin < 6 || isempty (M2))
+    if (isempty (M2))
       M2m1x = @(x) x;
     elseif (ischar (M2))
       M2m1x = str2func (M2);
@@ -116,11 +117,6 @@ function [x, flag, relres, iter, resvec]
 
     precon = @(x) M2m1x (M1m1x (x));
 
-    if (nargin < 7 || isempty (x0))
-      x0 = zeros (size (b));
-    endif
-
-
     x = x0;
 
     res = b - Ax (x);
@@ -220,4 +216,3 @@ endfunction
 %! b = sum (A, 2);
 %! [x, flag, relres, iter, resvec] = cgs (A, b, tol, [], diag (diag (A)));
 %! assert (x, ones (size (b)), 1e-7);
-
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/gmres.m
--- a/scripts/sparse/gmres.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/gmres.m	Fri Apr 29 15:22:54 2016 +0200
@@ -71,11 +71,15 @@
 ## @seealso{bicg, bicgstab, cgs, pcg, pcr, qmr}
 ## @end deftypefn
 
-function [x, flag, relres, it, resvec] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
+function [x, flag, relres, it, resvec] = gmres (A, b, varargin) 
 
-  if (nargin < 2 || nargin > 8)
+  if (nargin < 2)
     print_usage ();
   endif
+ 
+  if (length (varargin) < 6)
+    varargin{6} = [];
+  endif
 
   if (ischar (A))
     Ax = str2func (A);
@@ -87,19 +91,16 @@ function [x, flag, relres, it, resvec] =
     error ("gmres: A must be a function or square matrix");
   endif
 
-  if (nargin < 3 || isempty (restart))
-    restart = rows (b);
-  endif
+# First check restart since maxit depends from it
+  restart = __default__input__ ({rows(b)}, varargin{1}); 
 
-  if (nargin < 4 || isempty (rtol))
-    rtol = 1e-6;
-  endif
+# Check the input values of rtol, maxit, x0
+  [rtol, maxit, x0] = __default__input__ ({1e-06, min(rows (b)/restart, 10), zeros(size(b))}, varargin{2}, varargin{3}, varargin{6});
 
-  if (nargin < 5 || isempty (maxit))
-    maxit = min (rows (b)/restart, 10);
-  endif
+  M1 = varargin {4};
+  M2 = varargin {5};
 
-  if (nargin < 6 || isempty (M1))
+  if (isempty (M1))
     M1m1x = @(x) x;
   elseif (ischar (M1))
     M1m1x = str2func (M1);
@@ -111,7 +112,7 @@ function [x, flag, relres, it, resvec] =
     error ("gmres: preconditioner M1 must be a function or matrix");
   endif
 
-  if (nargin < 7 || isempty (M2))
+  if (isempty (M2))
     M2m1x = @(x) x;
   elseif (ischar (M2))
     M2m1x = str2func (M2);
@@ -125,10 +126,6 @@ function [x, flag, relres, it, resvec] =
 
   Pm1x = @(x) M2m1x (M1m1x (x));
 
-  if (nargin < 8 || isempty (x0))
-    x0 = zeros (size (b));
-  endif
-
   x_old = x0;
   x = x_old;
   prec_res = Pm1x (b - Ax (x_old));
@@ -228,9 +225,7 @@ endfunction
 
 
 %!error gmres (1)
-%!error gmres (1,2,3,4,5,6,7,8,9)
 %!error <A must be> gmres ({1},2)
 %!error <A must be a function or square matrix> gmres ({1},2)
 %!error <M1 must be a function or matrix> gmres (1,2,3,4,5,{6})
 %!error <M2 must be a function or matrix> gmres (1,2,3,4,5,6,{7})
-
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/pcg.m
--- a/scripts/sparse/pcg.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/pcg.m	Fri Apr 29 15:22:54 2016 +0200
@@ -216,38 +216,32 @@
 ## Modified by: Vittoria Rezzonico <vittoria.rezzonico@epfl.ch>
 ##  - Add the ability to provide the pre-conditioner as two separate matrices
 
-function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, tol, maxit, m1, m2, x0, varargin)
+function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, varargin) 
 
   ## M = M1*M2
-
-  if (nargin < 7 || isempty (x0))
-    x = zeros (size (b));
-  else
-    x = x0;
+  if (length (varargin) < 6)
+    varargin {6} = [];
   endif
 
-  if (nargin < 5 || isempty (m1))
+  m1 = varargin {3};
+  m2 = varargin {4};
+
+  [tol, maxit, x] = __default__input__ ({1e-06, min( rows(b), 20) , zeros( size (b))}, varargin{1}, varargin{2}, varargin{5});
+
+  if (isempty (m1))
      exist_m1 = 0;
   else
      exist_m1 = 1;
   endif
 
-  if (nargin < 6 || isempty (m2))
+  if (isempty (m2))
      exist_m2 = 0;
   else
      exist_m2 = 1;
   endif
-
-  if (nargin < 4 || isempty (maxit))
-    maxit = min (rows (b), 20);
-  endif
-
+  
   maxit += 2;
-
-  if (nargin < 3 || isempty (tol))
-    tol = 1e-6;
-  endif
-
+  
   preconditioned_residual_out = false;
   if (nargout > 5)
     T = zeros (maxit, maxit);
@@ -264,7 +258,7 @@ function [x, flag, relres, iter, resvec,
     r = b - A*x;
   else
     ## A should be a function.
-    r = b - feval (A, x, varargin{:});
+    r = b - feval (A, x, varargin{6:end});
   endif
 
   b_norm = norm (b);
@@ -277,7 +271,7 @@ function [x, flag, relres, iter, resvec,
       if (isnumeric (m1))
         y = m1 \ r;
       else
-        y = feval (m1, r, varargin{:});
+        y = feval (m1, r, varargin{6:end});
       endif
     else
       y = r;
@@ -286,7 +280,7 @@ function [x, flag, relres, iter, resvec,
       if (isnumeric (m2))
         z = m2 \ y;
       else
-        z = feval (m2, y, varargin{:});
+        z = feval (m2, y, varargin{6:end});
       endif
     else
       z = y;
@@ -301,7 +295,7 @@ function [x, flag, relres, iter, resvec,
       w = A * p;
     else
       ## A should be a function.
-      w = feval (A, p, varargin{:});
+      w = feval (A, p, varargin{6:end});
     endif
     ## Needed only for eigest.
     oldalpha = alpha;
@@ -343,7 +337,7 @@ function [x, flag, relres, iter, resvec,
       if (isnumeric (m1))
         y = m1 \ r;
       else
-        y = feval (m1, r, varargin{:});
+        y = feval (m1, r, varargin{6:end});
       endif
     else
       y = r;
@@ -352,7 +346,7 @@ function [x, flag, relres, iter, resvec,
       if (isnumeric (m2))
         z = m2 \ y;
       else
-        z = feval (m2, y, varargin{:});
+        z = feval (m2, y, varargin{6:end});
       endif
     else
       z = y;
@@ -533,6 +527,6 @@ endfunction
 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
 %! assert (norm (x - X) / norm (X), 0, 1e-6);
 %! assert (flag, 0);
-%! assert (iter, 1);  # should converge in one iteration
+%! assert (iter, 1); # should converge in one iteration
 %! assert (isnan (eigest), isnan ([NaN, NaN]));
 
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/pcr.m
--- a/scripts/sparse/pcr.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/pcr.m	Fri Apr 29 15:22:54 2016 +0200
@@ -172,30 +172,20 @@
 
 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
 
-function [x, flag, relres, iter, resvec] = pcr (A, b, tol, maxit, m, x0, varargin)
+function [x, flag, relres, iter, resvec] = pcr (A, b, varargin)
 
   breakdown = false;
 
-  if (nargin < 6 || isempty (x0))
-    x = zeros (size (b));
-  else
-    x = x0;
+  if (length (varargin) < 5)
+    varargin {5} = [];
   endif
 
-  if (nargin < 5)
-    m = [];
-  endif
+  [tol, maxit, x] = __default__input__ ({1e-6, 20, zeros(size (b))}, varargin {1}, varargin {2}, varargin {4});
 
-  if (nargin < 4 || isempty (maxit))
-    maxit = 20;
-  endif
+  m = varargin {3};
 
   maxit += 2;
 
-  if (nargin < 3 || isempty (tol))
-    tol = 1e-6;
-  endif
-
   if (nargin < 2)
     print_usage ();
   endif
@@ -204,7 +194,7 @@ function [x, flag, relres, iter, resvec]
   if (isnumeric (A))            # is A a matrix?
     r = b - A*x;
   else                          # then A should be a function!
-    r = b - feval (A, x, varargin{:});
+    r = b - feval (A, x, varargin{5:end});
   endif
 
   if (isnumeric (m))            # is M a matrix?
@@ -214,7 +204,7 @@ function [x, flag, relres, iter, resvec]
       p = m \ r;
     endif
   else                          # then M should be a function!
-    p = feval (m, r, varargin{:});
+    p = feval (m, r, varargin{5:end});
   endif
 
   iter = 2;
@@ -225,7 +215,7 @@ function [x, flag, relres, iter, resvec]
   if (isnumeric (A))            # is A a matrix?
     q = A * p;
   else                          # then A should be a function!
-    q = feval (A, p, varargin{:});
+    q = feval (A, p, varargin{5:end});
   endif
 
   resvec(1) = abs (norm (r));
@@ -240,7 +230,7 @@ function [x, flag, relres, iter, resvec]
         s = m \ q;
       endif
     else                        # then M should be a function!
-      s = feval (m, q, varargin{:});
+      s = feval (m, q, varargin{5:end});
     endif
     b_top = r' * s;
     b_bot = q' * s;
@@ -257,7 +247,7 @@ function [x, flag, relres, iter, resvec]
     if (isnumeric (A))          # is A a matrix?
       t = A*s;
     else                        # then A should be a function!
-      t = feval (A, s, varargin{:});
+      t = feval (A, s, varargin{5:end});
     endif
 
     alpha0 = (t'*s) / b_bot;
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/private/__default__input__.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/private/__default__input__.m	Fri Apr 29 15:22:54 2016 +0200
@@ -0,0 +1,50 @@
+## Copyright (C) 2016 Cristiano Dorigo 
+##
+## 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  {} {@var{[varargout]} =} __default__input__ (@var{def_val}, @var{varargin})
+## Check if the arguments in input of a function are empty or missing and in such cases sets up them in default values
+## 
+## The input argoments are:
+## @itemize @minus
+## @item @var{def_val} is a cell array that contains the values to use as default
+## @item @var{varargin} are the input argument
+## @end itemize
+##
+## The output argoments:
+## @itemize @minus
+## @item @var{varargout} all the input argument with filled the empty or missing paramenters
+## 
+## @end itemize
+##
+## @end deftypefn
+
+
+function [varargout] = __default__input__ (def_val, varargin)
+
+  m = length (def_val);
+  n = length (varargin);
+
+  for i = 1:m
+    if ( n < i || isempty (varargin {i}) )
+      varargout {i} = def_val {i};
+    else
+      varargout {i} = varargin {i};
+    endif
+  endfor
+
diff -r 7bdfa98a98a3 -r 8208b851edb6 scripts/sparse/qmr.m
--- a/scripts/sparse/qmr.m	Mon Apr 25 08:07:36 2016 -0700
+++ b/scripts/sparse/qmr.m	Fri Apr 29 15:22:54 2016 +0200
@@ -86,7 +86,7 @@
 
 ## Author: Nathan Podlich <nathan.podlich@gmail.com>
 
-function [x, flag, relres, iter, resvec] = qmr (A, b, tol, maxit, M1, M2, x0)
+function [x, flag, relres, iter, resvec] = qmr (A, b, varargin)
 
   if (nargin >= 2 && isvector (full (b)))
 
@@ -103,16 +103,15 @@ function [x, flag, relres, iter, resvec]
     else
       error ("qmr: A must be a square matrix or function");
     endif
-
-    if (nargin < 3 || isempty (tol))
-      tol = 1e-6;
+    
+    if (length (varargin) < 5)
+      varargin {5} = [];
     endif
 
-    if (nargin < 4 || isempty (maxit))
-      maxit = min (rows (b), 20);
-    else
-      maxit = fix (maxit);
-    endif
+    [tol, maxit, x] = __default__input__ ({1e-6, min(rows (b), 20), zeros(size (b))}, varargin {1}, varargin {2}, varargin {5});
+    
+    M1 = varargin {3};
+    M2 = varargin {4};
 
     if (nargin < 5 || isempty (M1))
       M1m1x = @(x, ignore) x;
@@ -148,12 +147,6 @@ function [x, flag, relres, iter, resvec]
       error ("qmr: preconditioner M2 must be a function or matrix");
     endif
 
-    if (nargin < 7 || isempty (x0))
-      x = zeros (size (b));
-    else
-      x = x0;
-    endif
-
     r = b - Ax (x);
 
     bnorm = norm (b);
