GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
amd.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2008-2017 David Bateman
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 // This is the octave interface to amd, which bore the copyright given
24 // in the help of the functions.
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <cstdlib>
31 
32 #include <string>
33 #include <vector>
34 
35 #include "ov.h"
36 #include "defun-dld.h"
37 #include "errwarn.h"
38 #include "pager.h"
39 #include "ov-re-mat.h"
40 
41 #include "ov-re-sparse.h"
42 #include "ov-cx-sparse.h"
43 #include "oct-map.h"
44 
45 #include "oct-sparse.h"
46 #include "oct-locbuf.h"
47 
48 #if defined (OCTAVE_ENABLE_64)
49 # define AMD_NAME(name) amd_l ## name
50 #else
51 # define AMD_NAME(name) amd ## name
52 #endif
53 
54 DEFUN_DLD (amd, args, nargout,
55  doc: /* -*- texinfo -*-
56 @deftypefn {} {@var{p} =} amd (@var{S})
57 @deftypefnx {} {@var{p} =} amd (@var{S}, @var{opts})
58 
59 Return the approximate minimum degree permutation of a matrix.
60 
61 This is a permutation such that the Cholesky@tie{}factorization of
62 @code{@var{S} (@var{p}, @var{p})} tends to be sparser than the
63 Cholesky@tie{}factorization of @var{S} itself. @code{amd} is typically
64 faster than @code{symamd} but serves a similar purpose.
65 
66 The optional parameter @var{opts} is a structure that controls the behavior
67 of @code{amd}. The fields of the structure are
68 
69 @table @asis
70 @item @var{opts}.dense
71 Determines what @code{amd} considers to be a dense row or column of the
72 input matrix. Rows or columns with more than @code{max (16, (dense *
73 sqrt (@var{n})))} entries, where @var{n} is the order of the matrix @var{S},
74 are ignored by @code{amd} during the calculation of the permutation.
75 The value of dense must be a positive scalar and the default value is 10.0
76 
77 @item @var{opts}.aggressive
78 If this value is a nonzero scalar, then @code{amd} performs aggressive
79 absorption. The default is not to perform aggressive absorption.
80 @end table
81 
82 The author of the code itself is Timothy A. Davis
83 @email{davis@@cise.ufl.edu}, University of Florida
84 (see @url{http://www.cise.ufl.edu/research/sparse/amd}).
85 @seealso{symamd, colamd}
86 @end deftypefn */)
87 {
88 #if defined (HAVE_AMD)
89 
90  int nargin = args.length ();
91 
92  if (nargin < 1 || nargin > 2)
93  print_usage ();
94 
95  octave_idx_type n_row, n_col;
96  const octave_idx_type *ridx, *cidx;
97  SparseMatrix sm;
99 
100  if (args(0).is_sparse_type ())
101  {
102  if (args(0).is_complex_type ())
103  {
104  scm = args(0).sparse_complex_matrix_value ();
105  n_row = scm.rows ();
106  n_col = scm.cols ();
107  ridx = scm.xridx ();
108  cidx = scm.xcidx ();
109  }
110  else
111  {
112  sm = args(0).sparse_matrix_value ();
113  n_row = sm.rows ();
114  n_col = sm.cols ();
115  ridx = sm.xridx ();
116  cidx = sm.xcidx ();
117  }
118  }
119  else
120  {
121  if (args(0).is_complex_type ())
122  sm = SparseMatrix (real (args(0).complex_matrix_value ()));
123  else
124  sm = SparseMatrix (args(0).matrix_value ());
125 
126  n_row = sm.rows ();
127  n_col = sm.cols ();
128  ridx = sm.xridx ();
129  cidx = sm.xcidx ();
130  }
131 
132  if (n_row != n_col)
133  err_square_matrix_required ("amd", "S");
134 
135  OCTAVE_LOCAL_BUFFER (double, Control, AMD_CONTROL);
136  AMD_NAME (_defaults) (Control);
137  if (nargin > 1)
138  {
139  octave_scalar_map arg1 = args(1).xscalar_map_value ("amd: OPTS argument must be a scalar structure");
140 
142 
143  tmp = arg1.getfield ("dense");
144  if (tmp.is_defined ())
145  Control[AMD_DENSE] = tmp.double_value ();
146 
147  tmp = arg1.getfield ("aggressive");
148  if (tmp.is_defined ())
149  Control[AMD_AGGRESSIVE] = tmp.double_value ();
150  }
151 
153  Matrix xinfo (AMD_INFO, 1);
154  double *Info = xinfo.fortran_vec ();
155 
156  // FIXME: how can we manage the memory allocation of amd
157  // in a cleaner manner?
158  SUITESPARSE_ASSIGN_FPTR (malloc_func, amd_malloc, malloc);
159  SUITESPARSE_ASSIGN_FPTR (free_func, amd_free, free);
160  SUITESPARSE_ASSIGN_FPTR (calloc_func, amd_calloc, calloc);
161  SUITESPARSE_ASSIGN_FPTR (realloc_func, amd_realloc, realloc);
162  SUITESPARSE_ASSIGN_FPTR (printf_func, amd_printf, printf);
163 
164  octave_idx_type result = AMD_NAME (_order) (n_col, cidx, ridx, P,
165  Control, Info);
166 
167  if (result == AMD_OUT_OF_MEMORY)
168  error ("amd: out of memory");
169  else if (result == AMD_INVALID)
170  error ("amd: matrix S is corrupted");
171 
172  Matrix Pout (1, n_col);
173  for (octave_idx_type i = 0; i < n_col; i++)
174  Pout.xelem (i) = P[i] + 1;
175 
176  if (nargout > 1)
177  return ovl (Pout, xinfo);
178  else
179  return ovl (Pout);
180 
181 #else
182 
183  octave_unused_parameter (args);
184  octave_unused_parameter (nargout);
185 
186  err_disabled_feature ("amd", "AMD");
187 
188 #endif
189 }
190 
191 /*
192 %!shared A, A2, opts
193 %! A = ones (20, 30);
194 %! A2 = ones (30, 30);
195 %!
196 %!testif HAVE_AMD
197 %! assert(amd (A2), [1:30]);
198 %! opts.dense = 25;
199 %! assert(amd (A2, opts), [1:30]);
200 %! opts.aggressive = 1;
201 %! assert(amd (A2, opts), [1:30]);
202 
203 %!error <S must be a square matrix|was unavailable or disabled> amd (A)
204 %!error amd (A2, 2)
205 %!error <matrix S is corrupted|was unavailable or disabled> amd ([])
206 */
octave_idx_type * xridx(void)
Definition: Sparse.h:536
octave_idx_type cols(void) const
Definition: Sparse.h:272
octave_idx_type rows(void) const
Definition: Sparse.h:271
ar P
Definition: __luinc__.cc:49
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
nd group nd example oindent or xample printf("%s\n", nthargout(2,"system","cmd"))
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
bool is_defined(void) const
Definition: ov.h:536
octave_idx_type * xcidx(void)
Definition: Sparse.h:549
void error(const char *fmt,...)
Definition: error.cc:570
#define AMD_NAME(name)
Definition: amd.cc:51
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:112
JNIEnv void * args
Definition: ov-java.cc:67
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:935
int nargin
Definition: graphics.cc:10115
double tmp
Definition: data.cc:6300
Definition: dMatrix.h:37
With real return the complex result
Definition: data.cc:3375
T & xelem(octave_idx_type n)
Definition: Array.h:455
void free(void *)
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
octave_value getfield(const std::string &key) const
Definition: oct-map.cc:171
#define DEFUN_DLD(name, args_name, nargout_name, doc)
Definition: defun-dld.h:45
void * malloc(size_t)
const T * fortran_vec(void) const
Definition: Array.h:584
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:136
double double_value(bool frc_str_conv=false) const
Definition: ov.h:775
void err_disabled_feature(const std::string &fcn, const std::string &feature, const std::string &pkg)
Definition: errwarn.cc:50