[CHANGESET] Add the amd function

[David Bateman]

As the amd code is already used with umfpack, which Octave uses for the
sparse LU function, and Tim Davis supplies a mex file for the equivalent
of the amd function in matlab under the GPL, it is trivial to write an
Octave version of the same code.. Attached is a matlab compatible
version of the amd function for Octave.

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 1206102654 -3600
# Node ID 5429d6a96b4e39a1db5fbaf3692e059bed23be28
# Parent  6e930aebeebffd2fcd91fe94d7c47269df38d76a
Add the amd function

diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@ 2008-03-18  David Bateman  <address@hidden
+2008-03-21  David Bateman  <address@hidden>
+
+       * configure.in (HAVE_AMD): Complete test for presence of amd.
+
 2008-03-18  David Bateman  <address@hidden>
 
        * configure.in (AC_CHECK_FUNCS): Also check lgamma_r.
diff --git a/configure.in b/configure.in
--- a/configure.in
+++ b/configure.in
@@ -734,7 +734,26 @@ AC_SUBST(LAPACK_DIR)
 # Check for AMD library
 AMD_LIBS=
 AC_SUBST(AMD_LIBS)
-AC_CHECK_LIB(amd, amd_postorder, [AMD_LIBS="-lamd"; 
with_amd=yes],[with_amd=no])
+
+AC_ARG_WITH(amd,
+  [AS_HELP_STRING([--without-amd],
+     [don't use AMD, disable some sparse functionality])],
+  with_amd=$withval, with_amd=yes)
+
+warn_amd="AMD not found. This will result in some lack of functionality for 
sparse matrices."
+if test "$with_amd" = yes; then
+  with_amd=no
+  AC_CHECK_HEADERS([suitesparse/amd.h ufsparse/amd.h amd/amd.h amd.h], [
+    AC_CHECK_LIB(amd, amd_postorder, [AMD_LIBS="-lamd"; with_amd=yes])
+    if test "$with_amd" = yes; then
+      AC_DEFINE(HAVE_AMD, 1, [Define if the AMD library is used.])
+      warn_amd=
+    fi
+    break])
+fi 
+if test -n "$warn_amd"; then
+  AC_MSG_WARN($warn_amd)
+fi
 
 # Check for CAMD library
 CAMD_LIBS=
@@ -1948,6 +1967,11 @@ if test -n "$warn_umfpack"; then
   warn_msg_printed=true
 fi
 
+if test -n "$warn_amd"; then
+  AC_MSG_WARN($warn_amd)
+  warn_msg_printed=true
+fi
+
 if test -n "$warn_colamd"; then
   AC_MSG_WARN($warn_colamd)
   warn_msg_printed=true
diff --git a/doc/ChangeLog b/doc/ChangeLog
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,7 @@ 2008-03-19  Michael D. Godfrey  <godfrey
+2008-03-21  David Bateman  <address@hidden>
+
+       * interpreter/sparse.txi: Document amd function.
+
 2008-03-19  Michael D. Godfrey  <address@hidden>
 
        * interpreter/plot.txi: Reorder symbol character table.
diff --git a/doc/interpreter/sparse.txi b/doc/interpreter/sparse.txi
--- a/doc/interpreter/sparse.txi
+++ b/doc/interpreter/sparse.txi
@@ -469,7 +469,7 @@ used.
 @c @dfn{treelayout}
 
 @item Sparse matrix reordering:
-  @dfn{ccolamd}, @dfn{colamd}, @dfn{colperm}, @dfn{csymamd},
+  @dfn{amd}, @dfn{ccolamd}, @dfn{colamd}, @dfn{colperm}, @dfn{csymamd},
   @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm}
 
 @item Linear algebra:
@@ -621,7 +621,7 @@ Several functions are available to reord
 Several functions are available to reorder depending on the type of the
 matrix to be factorized. If the matrix is symmetric positive-definite,
 then @dfn{symamd} or @dfn{csymamd} should be used. Otherwise
address@hidden or @dfn{ccolamd} should be used. For completeness
address@hidden, @dfn{colamd} or @dfn{ccolamd} should be used. For completeness
 the reordering functions @dfn{colperm} and @dfn{randperm} are
 also available.
 
@@ -706,6 +706,8 @@ Finally, Octave implicitly reorders the 
 Finally, Octave implicitly reorders the matrix when using the div (/)
 and ldiv (\) operators, and so no the user does not need to explicitly
 reorder the matrix to maximize performance.
+
address@hidden(amd)
 
 @DOCSTRING(ccolamd)
 
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,7 @@ 2008-03-20  David Bateman  <address@hidden
+2008-03-21  David Bateman  <address@hidden>
+
+       * oct-sparse.h: Add headers for amd.h.
+
 2008-03-20  David Bateman  <address@hidden>
 
        * Array.cc (Array<T> Array<T>::diag (octave_idx_type) const): New
diff --git a/liboctave/oct-sparse.h b/liboctave/oct-sparse.h
--- a/liboctave/oct-sparse.h
+++ b/liboctave/oct-sparse.h
@@ -25,6 +25,16 @@ along with Octave; see the file COPYING.
 
 #ifdef HAVE_CONFIG_H
 #include <config.h>
+#endif
+
+#if defined (HAVE_SUITESPARSE_AMD_H)
+#include <suitesparse/amd.h>
+#elif defined (HAVE_UFSPARSE_AMD_H)
+#include <ufsparse/amd.h>
+#elif defined (HAVE_AMD_AMD_H)
+#include <amd/amd.h>
+#elif defined (HAVE_AMD_H)
+#include <amd.h>
 #endif
 
 #if defined (HAVE_SUITESPARSE_UMFPACK_H)
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@ 2008-03-20  David Bateman  <address@hidden
+2008-03-21  David Bateman  <address@hidden>
+
+       * DLD-FUNCTIONS/amd.cc: New file.
+       * Makefile.in (DLD_XSRC): Add amd.cc.
+
 2008-03-20  David Bateman  <address@hidden>
 
        * Cell.cc (Cell Cell::diag (void) const): delete.
diff --git a/src/DLD-FUNCTIONS/amd.cc b/src/DLD-FUNCTIONS/amd.cc
new file mode 100644
--- /dev/null
+++ b/src/DLD-FUNCTIONS/amd.cc
@@ -0,0 +1,212 @@
+/*
+
+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/>.
+
+*/
+
+// This is the octave interface to amd, which bore the copyright given
+// in the help of the functions.
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cstdlib>
+
+#include <string>
+#include <vector>
+
+#include "ov.h"
+#include "defun-dld.h"
+#include "pager.h"
+#include "ov-re-mat.h"
+
+#include "ov-re-sparse.h"
+#include "ov-cx-sparse.h"
+#include "oct-map.h"
+
+#include "oct-sparse.h"
+
+#ifdef IDX_TYPE_LONG
+#define AMD_NAME(name) amd_l ## name
+#else
+#define AMD_NAME(name) amd ## name
+#endif
+
+DEFUN_DLD (amd, args, nargout,
+    "-*- texinfo -*-\n\
address@hidden {Loadable Function} address@hidden =} amd (@var{s})\n\
address@hidden {Loadable Function} address@hidden =} amd (@var{s}, 
@var{opts})\n\
+\n\
+Returns the approximate minimum degree permutation of a matrix. This\n\
+permutation such that the Cholesky factorization of @address@hidden 
(@var{p},\n\
address@hidden)} tends to be sparser than the Cholesky factorization of 
@var{s}\n\
+itself. @code{amd} is typically faster than @code{symamd} but serves a\n\
+similar purpose.\n\
+\n\
+The optional parameter @var{opts} is a structure that controls the\n\
+behavior of @code{amd}. The fields of these structure are\n\
+\n\
address@hidden @asis\n\
address@hidden opts.dense\n\
+Determines what @code{amd} considers to be a dense row or column of the\n\
+input matrix. Rows or columns with more that @code{max(16, (dense *\n\
+sqrt (@var{n})} entries, where @var{n} is the order of the matrix @var{s},\n\
+are igorned by @code{amd} during the calculation of the permutation\n\
+The value of dense must be a positive scalar and its default value is 10.0\n\
+\n\
address@hidden opts.aggressive\n\
+If this value is a non zero scalar, then @code{amd} performs agressive\n\
+absorption. The default is not to perform agressive absorption.\n\
address@hidden table\n\
+\n\
+The author of the code itself is Timothy A. Davis (davis@@cise.ufl.edu),\n\
+University of Florida (see 
@url{http://www.cise.ufl.edu/research/sparse/amd}).\n\
address@hidden, colamd}\n\
address@hidden deftypefn")
+{
+  octave_value_list retval;
+
+#ifdef HAVE_AMD
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+  else
+    {
+      octave_idx_type n_row, n_col, nnz;
+      const octave_idx_type *ridx, *cidx;
+      SparseMatrix sm;
+      SparseComplexMatrix scm;
+
+      if (args(0).is_sparse_type ())
+       {
+         if (args(0).is_complex_type ())
+           {
+             scm = args(0).sparse_complex_matrix_value ();
+             n_row = scm.rows ();
+             n_col = scm.cols ();
+             nnz = scm.nzmax ();
+             ridx = scm.xridx ();
+             cidx = scm.xcidx ();
+           }
+         else
+           {
+             sm = args(0).sparse_matrix_value ();
+             n_row = sm.rows ();
+             n_col = sm.cols ();
+             nnz = sm.nzmax ();
+             ridx = sm.xridx ();
+             cidx = sm.xcidx ();
+           }
+       }
+      else
+       {
+         if (args(0).is_complex_type ())
+           sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+         else
+           sm = SparseMatrix (args(0).matrix_value ());
+         
+         n_row = sm.rows ();
+         n_col = sm.cols ();
+         nnz = sm.nzmax ();
+         ridx = sm.xridx ();
+         cidx = sm.xcidx ();
+       }
+
+      if (!error_state && n_row != n_col)
+       error ("amd: input matrix must be square");
+
+      if (!error_state)
+       {
+         OCTAVE_LOCAL_BUFFER (double, Control, AMD_CONTROL);
+         AMD_NAME (_defaults) (Control) ;
+         if (nargin > 1)
+           {
+             Octave_map arg1 = args(1).map_value ();
+         
+             if (!error_state)
+               {
+                 if (arg1.contains ("dense"))
+                   {
+                     Cell c = arg1.contents ("dense");
+                     if (c.length() == 1)
+                       Control[AMD_DENSE] = c.elem(0).double_value ();
+                     else
+                       error ("amd: invalid options structure");
+                   }
+                 if (arg1.contains ("aggressive"))
+                   {
+                     Cell c = arg1.contents ("aggressive");
+                     if (c.length() == 1)
+                       Control[AMD_AGGRESSIVE] = c.elem(0).double_value ();
+                     else
+                       error ("amd: invalid options structure");
+                   }
+               }
+           }
+
+         if (!error_state)
+           {
+             OCTAVE_LOCAL_BUFFER (octave_idx_type, P, n_col);
+             Matrix xinfo (AMD_INFO, 1);
+             double *Info = xinfo.fortran_vec ();
+
+             // FIXME: How can we manage the memory allocation of amd in 
+             // a cleaner manner? 
+             amd_malloc = malloc;
+             amd_free = free;
+             amd_calloc = calloc;
+             amd_realloc = realloc;
+             amd_printf = printf;
+
+             octave_idx_type result = AMD_NAME (_order) (n_col, cidx, ridx, P,
+                                                         Control, Info);
+
+             switch (result)
+               {
+               case AMD_OUT_OF_MEMORY:
+                 error ("amd: out of memory");
+                 break;
+               case AMD_INVALID:
+                 error ("amd: input matrix is corrupted");
+                 break;
+               default:
+                 {
+                   if (nargout > 1)
+                     retval(1) = xinfo;
+
+                   Matrix Pout (1, n_col);
+                   for (octave_idx_type i = 0; i < n_col; i++)
+                     Pout.xelem (i) = P[i] + 1;
+
+                   retval (0) = Pout;
+                 }
+               }
+           }
+       }
+    }
+#else
+
+  error ("amd: not available in this version of Octave");
+
+#endif
+
+  return retval;
+}
diff --git a/src/Makefile.in b/src/Makefile.in
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -62,8 +62,8 @@ OPT_IN := $(addprefix ../liboctave/, $(a
 OPT_IN := $(addprefix ../liboctave/, $(addsuffix .in, $(OPT_BASE)))
 OPT_INC := $(addprefix ../liboctave/, $(addsuffix .h, $(OPT_BASE)))
 
-DLD_XSRC := balance.cc besselj.cc betainc.cc bsxfun.cc cellfun.cc chol.cc \
-       ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
+DLD_XSRC := amd.cc balance.cc besselj.cc betainc.cc bsxfun.cc cellfun.cc \
+       chol.cc ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
        dasrt.cc dassl.cc det.cc dispatch.cc dlmread.cc dmperm.cc eig.cc \
        expm.cc fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc fsolve.cc \
        gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \