Mercurial > octave-nkf

changeset 7619:56012914972a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add the amd function
author David Bateman <dbateman@free.fr>
date Fri, 21 Mar 2008 13:20:11 -0400
parents 3209a584e1ac
children 36594d5bbe13
files ChangeLog configure.in doc/ChangeLog doc/interpreter/sparse.txi liboctave/ChangeLog liboctave/oct-sparse.h src/ChangeLog src/DLD-FUNCTIONS/amd.cc src/Makefile.in
diffstat 9 files changed, 270 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog	Thu Mar 20 20:54:19 2008 +0100
+++ b/ChangeLog	Fri Mar 21 13:20:11 2008 -0400
@@ -1,3 +1,7 @@
+2008-03-21  David Bateman  <dbateman@free.fr>
+
+	* configure.in (HAVE_AMD): Complete test for presence of amd.
+
 2008-03-18  David Bateman  <dbateman@free.fr>
 
 	* configure.in (AC_CHECK_FUNCS): Also check lgamma_r.
--- a/configure.in	Thu Mar 20 20:54:19 2008 +0100
+++ b/configure.in	Fri Mar 21 13:20:11 2008 -0400
@@ -734,7 +734,26 @@
 # 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 @@
   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
--- a/doc/ChangeLog	Thu Mar 20 20:54:19 2008 +0100
+++ b/doc/ChangeLog	Fri Mar 21 13:20:11 2008 -0400
@@ -1,3 +1,7 @@
+2008-03-21  David Bateman  <dbateman@free.fr>
+
+	* interpreter/sparse.txi: Document amd function.
+
 2008-03-19  Michael D. Godfrey  <godfrey@isl.stanford.edu>
 
 	* interpreter/plot.txi: Reorder symbol character table.
--- a/doc/interpreter/sparse.txi	Thu Mar 20 20:54:19 2008 +0100
+++ b/doc/interpreter/sparse.txi	Fri Mar 21 13:20:11 2008 -0400
@@ -469,7 +469,7 @@
 @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 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
-@dfn{colamd} or @dfn{ccolamd} should be used. For completeness
+@dfn{amd}, @dfn{colamd} or @dfn{ccolamd} should be used. For completeness
 the reordering functions @dfn{colperm} and @dfn{randperm} are
 also available.
 
@@ -707,6 +707,8 @@
 and ldiv (\) operators, and so no the user does not need to explicitly
 reorder the matrix to maximize performance.
 
+@DOCSTRING(amd)
+
 @DOCSTRING(ccolamd)
 
 @DOCSTRING(colamd)
--- a/liboctave/ChangeLog	Thu Mar 20 20:54:19 2008 +0100
+++ b/liboctave/ChangeLog	Fri Mar 21 13:20:11 2008 -0400
@@ -1,3 +1,7 @@
+2008-03-21  David Bateman  <dbateman@free.fr>
+
+	* oct-sparse.h: Add headers for amd.h.
+
 2008-03-19  John W. Eaton  <jwe@octave.org>
 
 	* oct-env.cc (octave_env::do_base_pathname): Also handle rooted
--- a/liboctave/oct-sparse.h	Thu Mar 20 20:54:19 2008 +0100
+++ b/liboctave/oct-sparse.h	Fri Mar 21 13:20:11 2008 -0400
@@ -27,6 +27,16 @@
 #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)
 #include <suitesparse/umfpack.h>
 #elif defined (HAVE_UFSPARSE_UMFPACK_H)
--- a/src/ChangeLog	Thu Mar 20 20:54:19 2008 +0100
+++ b/src/ChangeLog	Fri Mar 21 13:20:11 2008 -0400
@@ -1,3 +1,8 @@
+2008-03-21  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/amd.cc: New file.
+	* Makefile.in (DLD_XSRC): Add amd.cc.
+
 2008-03-20  David Bateman  <dbateman@free.fr>
 
 	* data.cc (static octave_value make_diag (const Cell&,
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/amd.cc	Fri Mar 21 13:20:11 2008 -0400
@@ -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\
+@deftypefn {Loadable Function} {@var{p} =} amd (@var{s})\n\
+@deftypefnx {Loadable Function} {@var{p} =} 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 @code{@var{s} (@var{p},\n\
+@var{p})} 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\
+@table @asis\n\
+@item 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\
+@item 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\
+@end 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\
+@seealso{symamd, colamd}\n\
+@end 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;
+}
--- a/src/Makefile.in	Thu Mar 20 20:54:19 2008 +0100
+++ b/src/Makefile.in	Fri Mar 21 13:20:11 2008 -0400
@@ -62,8 +62,8 @@
 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 \