Mercurial > octave-nkf

comparison liboctave/UMFPACK/AMD/OCTAVE/amd.cc @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 /*
2
3 Copyright (C) 2004 David Bateman
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful, but WITHOUT
11 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; see the file COPYING. If not, write to the Free
17 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18
19 In addition to the terms of the GPL, you are permitted to link
20 this program with any Open Source program, as defined by the
21 Open Source Initiative (www.opensource.org)
22
23 */
24
25 /*
26
27 This is the Octave interface to the UMFPACK AMD code, which bore the following
28 copyright
29
30 AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis,
31 Patrick R. Amestoy, and Iain S. Duff. See ../README for License.
32 email: davis@cise.ufl.edu CISE Department, Univ. of Florida.
33 web: http://www.cise.ufl.edu/research/sparse/amd
34 --------------------------------------------------------------------------
35
36 Acknowledgements: This work was supported by the National Science
37 Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270.
38
39 */
40
41 #include <cstdlib>
42 #include <string>
43
44 #include <octave/config.h>
45 #include <octave/ov.h>
46 #include <octave/defun-dld.h>
47 #include <octave/pager.h>
48 #include <octave/ov-re-mat.h>
49
50 #include "ov-re-sparse.h"
51 #include "ov-cx-sparse.h"
52
53 // External AMD functions in C
54 extern "C" {
55 #include "amd.h"
56 }
57
58 DEFUN_DLD (amd, args, nargout,
59 "-*- texinfo -*-\n\
60 @deftypefn {Loadable Function} {@var{p} =} amd (@var{s})\n\
61 @deftypefnx {Loadable Function} {@var{Control} =} amd ()\n\
62 @deftypefnx {Loadable Function} {[@var{p}, @var{info}] =} amd (@var{s})\n\
63 \n\
64 AMD Approximate minimum degree permutation. Returns the approximate\n\
65 minimum degree permutation vector for the sparse matrix\n\
66 @code{@var{c} = @var{S} + @var{S}'}. The Cholesky factorization of\n\
67 @code{@var{c} (@var{p}, @var{p})}, or @code{@var{s} (@var{p}, @var{p})},\n\
68 tends to be sparser than that of @var{c} or @var{s}.\n\
69 @var{s} must be square. If @var{s} is full, @code{amd (@var{S})} is\n\
70 equivalent to @code{amd (sparse (@var{s}))}.\n\
71 \n\
72 @table @asis\n\
73 @item @var{Control} (1)\n\
74 If S is n-by-n, then rows/columns with more than\n\
75 @code{@dfn{max} (16, (@var{Control} (1)) * @dfn{sqrt} (@var{n}))} entries\n\
76 in @code{@var{s} + @var{S}'} are considered @emph{dense}, and ignored during\n\
77 ordering. They are placed last in the output permutation. The default is\n\
78 10.0 if @var{Control} is not present.\n\
79 @item @var{Control} (2)\n\
80 If nonzero, then aggressive absorption is performed. This is the default if\n\
81 @var{Control} is not present.\n\
82 @item @var{Control} (3)\n\
83 If nonzero, print statistics about the ordering.\n\
84 @end table\n\
85 \n\
86 @table @asis\n\
87 @item @var{Info} (1)\n\
88 status (0: ok, -1: out of memory, -2: matrix invalid)\n\
89 @item @var{Info} (2)\n\
90 @code{@var{n} = size (@var{a}, 1)}\n\
91 @item @var{Info} (3)\n\
92 @code{nnz (A)}\n\
93 @item @var{Info} (4)\n\
94 The symmetry of the matrix @var{s} (0.0 means purely unsymmetric, 1.0 means\n\
95 purely symmetric). Computed as: @code{@var{b} = tril (@var{s}, -1) +\n\
96 triu (@var{s}, 1); @var{symmetry} = nnz (@var{b} & @var{b}')\n\
97 / nnz (@var{b});}\n\
98 @item @var{Info} (5)\n\
99 @code{nnz (diag (@var{s}))}\n\
100 @item @var{Info} (6)\n\
101 @dfn{nnz} in @code{@var{s} + @var{s}'}, excluding the diagonal\n\
102 (= @code{nnz (@var{b} + @var{b}')})\n\
103 @item @var{Info} (7)\n\
104 Number of @emph{dense} rows/columns in @code{@var{s} + @var{s}'}\n\
105 @item @var{Info} (8)\n\
106 The amount of memory used by AMD, in bytes\n\
107 @item @var{Info} (9)\n\
108 The number of memory compactions performed by AMD\n\
109 @end table\n\
110 \n\
111 The following statistics are slight upper bounds because of the\n\
112 approximate degree in AMD. The bounds are looser if @emph{dense}\n\
113 rows/columns are ignored during ordering @code{(@var{Info} (7) > 0)}.\n\
114 The statistics are for a subsequent factorization of the matrix\n\
115 @code{@var{c} (@var{p},@var{p})}. The LU factorization statistics assume\n \
116 no pivoting.\n\
117 \n\
118 @table @asis\n\
119 @item @var{Info} (10)\n\
120 The number of nonzeros in L, excluding the diagonal\n\
121 @item @var{Info} (11)\n\
122 The number of divide operations for LL', LDL', or LU\n\
123 @item @var{Info (12)}\n\
124 The number of multiply-subtract pairs for LL' or LDL'\n\
125 @item @var{Info} (13)\n\
126 The number of multiply-subtract pairs for LU\n\
127 @item @var{Info} (14)\n\
128 The max number of nonzeros in any column of L (incl. diagonal)\n\
129 @item @var{Info} (15:20)\n\
130 unused, reserved for future use\n\
131 @end table\n\
132 \n\
133 An assembly tree post-ordering is performed, which is typically the same\n\
134 as an elimination tree post-ordering. It is not always identical because\n\
135 of the approximate degree update used, and because @emph{dense} rows/columns\n\
136 do not take part in the post-order. It well-suited for a subsequent\n\
137 @dfn{chol}, however. If you require a precise elimination tree\n\
138 post-ordering, then do:\n\
139 \n\
140 @group\n\
141 @var{p} = @dfn{amd} (@var{s});\n\
142 % skip this if S already symmetric\n\
143 @var{c} = spones (@var{s}) + spones (@var{s}');\n\
144 [@var{ignore}, @var{q}] = sparsfun ('symetree', @var{c} (@var{p}, @var{p}));\n\
145 @var{p} = @var{p} (@var{q});\n\
146 @end group\n\
147 \n\
148 AMD Version 1.1 (Jan. 21, 2004), Copyright @copyright{} 2004 by\n\
149 Timothy A. Davis, Patrick R. Amestoy, and Iain S. Duff.\n\
150 \n\
151 email: davis@@cise.ufl.edu (CISE Department, Univ. of Florida).\n\
152 \n\
153 web: http://www.cise.ufl.edu/research/sparse/amd\n\
154 \n\
155 Acknowledgements: This work was supported by the National Science\n\
156 Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270.\n\
157 @end deftypefn")
158 {
159 int nargin = args.length ();
160 octave_value_list retval;
161 int spumoni = 0;
162
163 if (nargin > 2 || nargout > 2)
164 usage ("p = amd (A) or [p, Info] = amd (A, Control)");
165 else if (nargin == 0)
166 {
167 // Get the default control parameter, and return
168 NDArray control (dim_vector (AMD_CONTROL, 1));
169 double *control_ptr = control.fortran_vec ();
170 amd_defaults (control_ptr);
171 if (nargout == 0)
172 amd_control (control_ptr);
173 else
174 retval(0) = control;
175 }
176 else
177 {
178 NDArray control (dim_vector (AMD_CONTROL, 1));
179 double *control_ptr = control.fortran_vec ();
180 amd_defaults (control_ptr);
181
182 if (nargin > 1)
183 {
184 NDArray control_in = args(1).array_value();
185
186 if (error_state)
187 {
188 error ("amd: could not read control vector");
189 return retval;
190 }
191
192 dim_vector dv = control_in.dims ();
193 if (dv.length() > 2 || (dv(0) != 1 && dv(1) != 1))
194 {
195 error ("amd: control vector isn't a vector");
196 return retval;
197 }
198
199 int nc = dv.numel ();
200 control (AMD_DENSE) = (nc > 0 ? control_in (AMD_DENSE) :
201 AMD_DEFAULT_DENSE);
202 control (AMD_AGGRESSIVE) = (nc > 1 ? control_in (AMD_AGGRESSIVE) :
203 AMD_DEFAULT_AGGRESSIVE);
204 spumoni = (nc > 2 ? (control_in (2) != 0) : 0);
205 }
206
207 if (spumoni > 0)
208 amd_control (control_ptr);
209
210 int *Ap, *Ai;
211 int n, m, nz;
212
213 // These are here only so that the C++ destructors don't prematurally
214 // remove the underlying data we are interested in
215 SparseMatrix sm;
216 SparseComplexMatrix scm;
217 Matrix mm;
218 ComplexMatrix cm;
219
220 if (args(0).class_name () != "sparse" && spumoni > 0)
221 octave_stdout << " input matrix A is full (sparse copy"
222 << " of A will be created)" << std::endl;
223
224 if (args(0).is_complex_type ())
225 {
226 scm = args(0).sparse_complex_matrix_value ();
227 Ai = scm.ridx ();
228 Ap = scm.cidx ();
229 m = scm.rows ();
230 n = scm.cols ();
231 nz = scm.nnz ();
232 }
233 else
234 {
235 sm = args(0).sparse_matrix_value ();
236 Ai = sm.ridx ();
237 Ap = sm.cidx ();
238 m = sm.rows ();
239 n = sm.cols ();
240 nz = sm.nnz ();
241 }
242
243 if (spumoni > 0)
244 octave_stdout << " input matrix A is " << m << "-by-" << n
245 << std::endl;
246
247 if (m != n)
248 {
249 error ("amd: A must be square");
250 return retval;
251 }
252
253 if (spumoni > 0)
254 octave_stdout << " input matrix A has " << nz <<
255 " nonzero entries" << std::endl;
256
257 // allocate workspace for output permutation
258 Array<int> P(n+1);
259 int *P_ptr = P.fortran_vec ();
260 NDArray info (dim_vector (AMD_INFO, 1));
261 double *info_ptr = info.fortran_vec ();
262 int result;
263
264 // order the matrix
265 result = amd_order (n, Ap, Ai, P_ptr, control_ptr, info_ptr);
266
267 // print results (including return value)
268 if (spumoni > 0)
269 amd_info (info_ptr);
270
271 // check error conditions
272 if (result == AMD_OUT_OF_MEMORY)
273 error ("amd: out of memory");
274 else if (result == AMD_INVALID)
275 error ("amd: input matrix A is corrupted");
276 else
277 {
278 // copy the outputs to Octave
279
280 // output permutation, P
281 NDArray perm (dim_vector (1, n));
282 for (int i = 0; i < n; i++)
283 perm (i) = double (P(i) + 1); // 1-based indexing for Octave
284
285 retval (0) = perm;
286
287 // Info
288 if (nargout > 1)
289 retval (1) = info;
290 }
291 }
292 return retval;
293 }
294
295 /*
296 ;;; Local Variables: ***
297 ;;; mode: C++ ***
298 ;;; End: ***
299 */