[Chageset]: Add the spaugment function

[David Bateman]

Another easy function to implement. The help string and the example took
longer to write :-)

D.

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

# HG changeset patch
# User David Bateman <address@hidden>
# Date 1207154950 -7200
# Node ID 05a11ba8bfe58e9b8c8dcad9a44785acbde9a130
# Parent  08b6dd8d9187e122924617bc401d915cd17c89c1
Add the spaugment function

diff --git a/doc/ChangeLog b/doc/ChangeLog
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,7 @@ 2008-03-26  Rafael Laboissiere  <rafael@
+2008-04-02  David Bateman  <address@hidden>
+
+       * interpreter/sparse.txi: Document spaugment.
+
 2008-03-26  Rafael Laboissiere  <address@hidden>
 
        * interpreter/mkoctfile.1: Remove spurious whitespace before macros
diff --git a/doc/interpreter/sparse.txi b/doc/interpreter/sparse.txi
--- a/doc/interpreter/sparse.txi
+++ b/doc/interpreter/sparse.txi
@@ -474,7 +474,7 @@ used.
 
 @item Linear algebra:
   @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank}
address@hidden @dfn{spaugment}
+  @dfn{spaugment}
 @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now
 
 @item Iterative techniques:
@@ -824,6 +824,11 @@ used with care.
 @DOCSTRING(sprank)
 
 @DOCSTRING(symbfact)
+
+For non square matrices, the user can also utilize the @code{spaugment}
+function to find a least squares solution to a linear equation.
+
address@hidden(spaugment)
 
 @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse 
Matrices
 @section Iterative Techniques applied to sparse matrices
diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,8 @@ 2008-04-01  Julian Schnidder  <j.schnidd
+2008-04-02  David Bateman  <address@hidden>
+
+       * sparse/spaugment.m: New function
+       * sparse/Makefile.in (SOURCES): Add it here.
+       
 2008-04-01  Julian Schnidder  <address@hidden>
 
        * miscellaneous/perl.m: New function.
diff --git a/scripts/sparse/Makefile.in b/scripts/sparse/Makefile.in
--- a/scripts/sparse/Makefile.in
+++ b/scripts/sparse/Makefile.in
@@ -33,8 +33,8 @@ INSTALL_DATA = @INSTALL_DATA@
 INSTALL_DATA = @INSTALL_DATA@
 
 SOURCES = colperm.m etreeplot.m gplot.m nonzeros.m normest.m \
-  pcg.m pcr.m spalloc.m spconvert.m spdiags.m speye.m spfun.m \
-  sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \
+  pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \
+  spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \
   spvcat.m spy.m treeplot.m
 
 DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
diff --git a/scripts/sparse/spaugment.m b/scripts/sparse/spaugment.m
new file mode 100644
--- /dev/null
+++ b/scripts/sparse/spaugment.m
@@ -0,0 +1,97 @@
+## Copyright (C) 2008  David Bateman
+##
+## 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} address@hidden =} spaugment (@var{a}, @var{c})
+## Creates the augmented matrix of @var{a}. This is given by
+##
+## @example
+## address@hidden * eye(@var{m}, @var{m}),@var{a}; @var{a}', zeros(@var{n},
+## @var{n})]
+## @end example
+##
+## @noindent
+## This is related to the leasted squared solution of 
+## @address@hidden \\ @var{b}}, by
+## 
+## @example
+## @var{s} * [ @var{r} / @var{c}; x] = address@hidden, zeros(@var{n},
+## columns(@var{b})]
+## @end example
+##
+## @noindent
+## where @var{r} is the residual error
+##
+## @example
+## @var{r} = @var{b} - @var{a} * @var{x}
+## @end example
+##
+## As the matrix @var{s} is symmetric indefinite it can be factorized
+## with @code{lu}, and the minimum norm solution can therefore be found
+## without the need for a @code{qr} factorization. As the residual
+## error will be @code{zeros (@var{m}, @var{m})} for under determined
+## problems, and example can be 
+##
+## @example
+## @group
+## m = 11; n = 10; mn = max(m ,n);
+## a = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],[-1,0,1], m, n);
+## x0 = a \ ones (m,1);
+## s = spaugment (a);
+## [L, U, P, Q] = lu (s);
+## x1 = Q * (U \ (L \ (P  * [ones(m,1); zeros(n,1)])));
+## x1 = x1(end - n + 1 : end);
+## @end group
+## @end example
+##
+## To find the solution of an overdetermined problem needs an estimate
+## of the residual error @var{r} and so it is more complex to formulate
+## a minimum norm solution using the @code{spaugment} function.
+##
+## In general the left division operator is more stable and faster than
+## using the @code{spaugment} function.
+## @end deftypefn
+
+function s = spaugment (a, c)
+  if (nargin < 2)
+    if (issparse (a))
+      c = max (max (abs (a))) / 1000;
+    else
+      if (ndims (a) != 2)
+       error ("spaugment: expecting 2-dimenisional matrix")
+      else
+       c = max (abs (a(:))) / 1000;
+      endif
+    endif
+  elseif (!isscalar (c))
+    error ("spaugment: c must be a scalar");
+  endif
+
+  [m, n] = size (a);
+  s = [ c * speye(m, m), a; a', sparse(n, n)];
+endfunction
+
+%!test
+%! m = 11; n = 10; mn = max(m ,n);
+%! a = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],[-1,0,1], m, n);
+%! x0 = a \ ones (m,1);
+%! s = spaugment (a);
+%! [L, U, P, Q] = lu (s);
+%! x1 = Q * (U \ (L \ (P  * [ones(m,1); zeros(n,1)])));
+%! x1 = x1(end - n + 1 : end);
+%! assert (x1, x0, 1e-10)