--- a/doc/ChangeLog	Wed Apr 02 13:29:19 2008 -0400
+++ b/doc/ChangeLog	Wed Apr 02 14:08:28 2008 -0400
@@ -1,3 +1,7 @@
+2008-04-02  David Bateman  <dbateman@free.fr>
+
+	* interpreter/sparse.txi: Document spaugment.
+
 2008-03-26  Rafael Laboissiere  <rafael@debian.org>
 
 	* interpreter/mkoctfile.1: Remove spurious whitespace before macros

--- a/doc/interpreter/sparse.txi	Wed Apr 02 13:29:19 2008 -0400
+++ b/doc/interpreter/sparse.txi	Wed Apr 02 14:08:28 2008 -0400
@@ -474,7 +474,7 @@
 
 @item Linear algebra:
   @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank}
-@c @dfn{spaugment}
+  @dfn{spaugment}
 @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now
 
 @item Iterative techniques:
@@ -825,6 +825,11 @@
 
 @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.
+
+@DOCSTRING(spaugment)
+
 @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse Matrices
 @section Iterative Techniques applied to sparse matrices
 

--- a/scripts/ChangeLog	Wed Apr 02 13:29:19 2008 -0400
+++ b/scripts/ChangeLog	Wed Apr 02 14:08:28 2008 -0400
@@ -1,5 +1,8 @@
 2008-04-02  David Bateman  <dbateman@free.fr>
 
+	* sparse/spaugment.m: New function
+	* sparse/Makefile.in (SOURCES): Add it here.
+	
 	* plot/__gnuplot_ginput__.m: Use the gnuplot stream itself for
 	communication rather than a chat file if mkfifo is not available.
 	* plot/gnuplot_drawnow.m: Open stream with popen2 to allow two way

--- a/scripts/sparse/Makefile.in	Wed Apr 02 13:29:19 2008 -0400
+++ b/scripts/sparse/Makefile.in	Wed Apr 02 14:08:28 2008 -0400
@@ -33,8 +33,8 @@
 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))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/spaugment.m	Wed Apr 02 14:08:28 2008 -0400
@@ -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} {@var{s} =} spaugment (@var{a}, @var{c})
+## Creates the augmented matrix of @var{a}. This is given by
+##
+## @example
+## [@var{c} * 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 
+## @code{@var{a} \\ @var{b}}, by
+## 
+## @example
+## @var{s} * [ @var{r} / @var{c}; x] = [@var{b}, 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)