Mercurial > octave-nkf

comparison libinterp/dldfcn/colamd.cc @ 15195:2fc554ffbc28

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
split libinterp from src * libinterp: New directory. Move all files from src directory here except Makefile.am, main.cc, main-cli.cc, mkoctfile.in.cc, mkoctfilr.in.sh, octave-config.in.cc, octave-config.in.sh. * libinterp/Makefile.am: New file, extracted from src/Makefile.am. * src/Makefile.am: Delete everything except targets and definitions needed to build and link main and utility programs. * Makefile.am (SUBDIRS): Include libinterp in the list. * autogen.sh: Run config-module.sh in libinterp/dldfcn directory, not src/dldfcn directory. * configure.ac (AC_CONFIG_SRCDIR): Use libinterp/octave.cc, not src/octave.cc. (DL_LDFLAGS, LIBOCTINTERP): Use libinterp, not src. (AC_CONFIG_FILES): Include libinterp/Makefile in the list. * find-docstring-files.sh: Look in libinterp, not src. * gui/src/Makefile.am (liboctgui_la_CPPFLAGS): Find header files in libinterp, not src.
author John W. Eaton <jwe@octave.org>
date Sat, 18 Aug 2012 16:23:39 -0400
parents src/dldfcn/colamd.cc@000587f92082
children 6aafe87a3144
comparison
equal deleted inserted replaced
15194:0f0b795044c3 15195:2fc554ffbc28
1 /*
2
3 Copyright (C) 2004-2012 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21
22 */
23
24 // This is the octave interface to colamd, which bore the copyright given
25 // in the help of the functions.
26
27 #ifdef HAVE_CONFIG_H
28 #include <config.h>
29 #endif
30
31 #include <cstdlib>
32
33 #include <string>
34 #include <vector>
35
36 #include "ov.h"
37 #include "defun-dld.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
44 #include "oct-sparse.h"
45 #include "oct-locbuf.h"
46
47 #ifdef IDX_TYPE_LONG
48 #define COLAMD_NAME(name) colamd_l ## name
49 #define SYMAMD_NAME(name) symamd_l ## name
50 #else
51 #define COLAMD_NAME(name) colamd ## name
52 #define SYMAMD_NAME(name) symamd ## name
53 #endif
54
55 // The symmetric column elimination tree code take from the Davis LDL code.
56 // Copyright given elsewhere in this file.
57 static void
58 symetree (const octave_idx_type *ridx, const octave_idx_type *cidx,
59 octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n)
60 {
61 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n);
62 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
63 if (P)
64 // If P is present then compute Pinv, the inverse of P
65 for (octave_idx_type k = 0 ; k < n ; k++)
66 Pinv[P[k]] = k ;
67
68 for (octave_idx_type k = 0 ; k < n ; k++)
69 {
70 // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k)
71 Parent[k] = n ; // parent of k is not yet known
72 Flag[k] = k ; // mark node k as visited
73 octave_idx_type kk = (P) ? (P[k]) : (k) ; // kth original, or permuted, column
74 octave_idx_type p2 = cidx[kk+1] ;
75 for (octave_idx_type p = cidx[kk] ; p < p2 ; p++)
76 {
77 // A (i,k) is nonzero (original or permuted A)
78 octave_idx_type i = (Pinv) ? (Pinv[ridx[p]]) : (ridx[p]) ;
79 if (i < k)
80 {
81 // follow path from i to root of etree, stop at flagged node
82 for ( ; Flag[i] != k ; i = Parent[i])
83 {
84 // find parent of i if not yet determined
85 if (Parent[i] == n)
86 Parent[i] = k ;
87 Flag[i] = k ; // mark i as visited
88 }
89 }
90 }
91 }
92 }
93
94 // The elimination tree post-ordering code below is taken from SuperLU
95 static inline octave_idx_type
96 make_set (octave_idx_type i, octave_idx_type *pp)
97 {
98 pp[i] = i;
99 return i;
100 }
101
102 static inline octave_idx_type
103 link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp)
104 {
105 pp[s] = t;
106 return t;
107 }
108
109 static inline octave_idx_type
110 find (octave_idx_type i, octave_idx_type *pp)
111 {
112 register octave_idx_type p, gp;
113
114 p = pp[i];
115 gp = pp[p];
116
117 while (gp != p)
118 {
119 pp[i] = gp;
120 i = gp;
121 p = pp[i];
122 gp = pp[p];
123 }
124
125 return p;
126 }
127
128 static octave_idx_type
129 etdfs (octave_idx_type v, octave_idx_type *first_kid,
130 octave_idx_type *next_kid, octave_idx_type *post,
131 octave_idx_type postnum)
132 {
133 for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w])
134 postnum = etdfs (w, first_kid, next_kid, post, postnum);
135
136 post[postnum++] = v;
137
138 return postnum;
139 }
140
141 static void
142 tree_postorder (octave_idx_type n, octave_idx_type *parent,
143 octave_idx_type *post)
144 {
145 // Allocate storage for working arrays and results
146 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1);
147 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1);
148
149 // Set up structure describing children
150 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1)
151 /* do nothing */;
152
153 for (octave_idx_type v = n-1; v >= 0; v--)
154 {
155 octave_idx_type dad = parent[v];
156 next_kid[v] = first_kid[dad];
157 first_kid[dad] = v;
158 }
159
160 // Depth-first search from dummy root vertex #n
161 etdfs (n, first_kid, next_kid, post, 0);
162 }
163
164 static void
165 coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg,
166 octave_idx_type *colend, octave_idx_type *parent,
167 octave_idx_type nr, octave_idx_type nc)
168 {
169 OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc);
170 OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc);
171 OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr);
172
173 // Compute firstcol[row] = first nonzero column in row
174 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc)
175 /* do nothing */;
176
177 for (octave_idx_type col = 0; col < nc; col++)
178 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
179 {
180 octave_idx_type row = ridx[p];
181 if (firstcol[row] > col)
182 firstcol[row] = col;
183 }
184
185 // Compute etree by Liu's algorithm for symmetric matrices,
186 // except use (firstcol[r],c) in place of an edge (r,c) of A.
187 // Thus each row clique in A'*A is replaced by a star
188 // centered at its first vertex, which has the same fill.
189 for (octave_idx_type col = 0; col < nc; col++)
190 {
191 octave_idx_type cset = make_set (col, pp);
192 root[cset] = col;
193 parent[col] = nc;
194 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
195 {
196 octave_idx_type row = firstcol[ridx[p]];
197 if (row >= col)
198 continue;
199 octave_idx_type rset = find (row, pp);
200 octave_idx_type rroot = root[rset];
201 if (rroot != col)
202 {
203 parent[rroot] = col;
204 cset = link (cset, rset, pp);
205 root[cset] = col;
206 }
207 }
208 }
209 }
210
211 DEFUN_DLD (colamd, args, nargout,
212 "-*- texinfo -*-\n\
213 @deftypefn {Loadable Function} {@var{p} =} colamd (@var{S})\n\
214 @deftypefnx {Loadable Function} {@var{p} =} colamd (@var{S}, @var{knobs})\n\
215 @deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} colamd (@var{S})\n\
216 @deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} colamd (@var{S}, @var{knobs})\n\
217 \n\
218 Column approximate minimum degree permutation.\n\
219 @code{@var{p} = colamd (@var{S})} returns the column approximate minimum\n\
220 degree permutation vector for the sparse matrix @var{S}. For a\n\
221 non-symmetric matrix @var{S}, @code{@var{S}(:,@var{p})} tends to have\n\
222 sparser LU@tie{}factors than @var{S}. The Cholesky@tie{}factorization of\n\
223 @code{@var{S}(:,@var{p})' * @var{S}(:,@var{p})} also tends to be sparser\n\
224 than that of @code{@var{S}' * @var{S}}.\n\
225 \n\
226 @var{knobs} is an optional one- to three-element input vector. If @var{S} is\n\
227 m-by-n, then rows with more than @code{max(16,@var{knobs}(1)*sqrt(n))}\n\
228 entries are ignored. Columns with more than\n\
229 @code{max (16,@var{knobs}(2)*sqrt(min(m,n)))} entries are removed prior to\n\
230 ordering, and ordered last in the output permutation @var{p}. Only\n\
231 completely dense rows or columns are removed if @code{@var{knobs}(1)} and\n\
232 @code{@var{knobs}(2)} are < 0, respectively. If @code{@var{knobs}(3)} is\n\
233 nonzero, @var{stats} and @var{knobs} are printed. The default is\n\
234 @code{@var{knobs} = [10 10 0]}. Note that @var{knobs} differs from earlier\n\
235 versions of colamd.\n\
236 \n\
237 @var{stats} is an optional 20-element output vector that provides data\n\
238 about the ordering and the validity of the input matrix @var{S}. Ordering\n\
239 statistics are in @code{@var{stats}(1:3)}. @code{@var{stats}(1)} and\n\
240 @code{@var{stats}(2)} are the number of dense or empty rows and columns\n\
241 ignored by @sc{colamd} and @code{@var{stats}(3)} is the number of garbage\n\
242 collections performed on the internal data structure used by @sc{colamd}\n\
243 (roughly of size @code{2.2 * nnz(@var{S}) + 4 * @var{m} + 7 * @var{n}}\n\
244 integers).\n\
245 \n\
246 Octave built-in functions are intended to generate valid sparse matrices,\n\
247 with no duplicate entries, with ascending row indices of the nonzeros\n\
248 in each column, with a non-negative number of entries in each column (!)\n\
249 and so on. If a matrix is invalid, then @sc{colamd} may or may not be able\n\
250 to continue. If there are duplicate entries (a row index appears two or\n\
251 more times in the same column) or if the row indices in a column are out\n\
252 of order, then @sc{colamd} can correct these errors by ignoring the duplicate\n\
253 entries and sorting each column of its internal copy of the matrix\n\
254 @var{S} (the input matrix @var{S} is not repaired, however). If a matrix\n\
255 is invalid in other ways then @sc{colamd} cannot continue, an error message\n\
256 is printed, and no output arguments (@var{p} or @var{stats}) are returned.\n\
257 @sc{colamd} is thus a simple way to check a sparse matrix to see if it's\n\
258 valid.\n\
259 \n\
260 @code{@var{stats}(4:7)} provide information if COLAMD was able to\n\
261 continue. The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if\n\
262 invalid. @code{@var{stats}(5)} is the rightmost column index that is\n\
263 unsorted or contains duplicate entries, or zero if no such column exists.\n\
264 @code{@var{stats}(6)} is the last seen duplicate or out-of-order row\n\
265 index in the column index given by @code{@var{stats}(5)}, or zero if no\n\
266 such row index exists. @code{@var{stats}(7)} is the number of duplicate\n\
267 or out-of-order row indices. @code{@var{stats}(8:20)} is always zero in\n\
268 the current version of @sc{colamd} (reserved for future use).\n\
269 \n\
270 The ordering is followed by a column elimination tree post-ordering.\n\
271 \n\
272 The authors of the code itself are Stefan I. Larimore and Timothy A.\n\
273 Davis @email{davis@@cise.ufl.edu}, University of Florida. The algorithm was\n\
274 developed in collaboration with John Gilbert, Xerox PARC, and Esmond\n\
275 Ng, Oak Ridge National Laboratory. (see\n\
276 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\
277 @seealso{colperm, symamd, ccolamd}\n\
278 @end deftypefn")
279 {
280 octave_value_list retval;
281
282 #ifdef HAVE_COLAMD
283
284 int nargin = args.length ();
285 int spumoni = 0;
286
287 if (nargout > 2 || nargin < 1 || nargin > 2)
288 print_usage ();
289 else
290 {
291 // Get knobs
292 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
293 COLAMD_NAME (_set_defaults) (knobs);
294
295 // Check for user-passed knobs
296 if (nargin == 2)
297 {
298 NDArray User_knobs = args(1).array_value ();
299 int nel_User_knobs = User_knobs.length ();
300
301 if (nel_User_knobs > 0)
302 knobs[COLAMD_DENSE_ROW] = User_knobs(0);
303 if (nel_User_knobs > 1)
304 knobs[COLAMD_DENSE_COL] = User_knobs(1) ;
305 if (nel_User_knobs > 2)
306 spumoni = static_cast<int> (User_knobs(2));
307
308 // print knob settings if spumoni is set
309 if (spumoni)
310 {
311
312 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "."
313 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n";
314
315 if (knobs[COLAMD_DENSE_ROW] >= 0)
316 octave_stdout << "knobs(1): " << User_knobs (0)
317 << ", rows with > max (16,"
318 << knobs[COLAMD_DENSE_ROW] << "*sqrt (size(A,2)))"
319 << " entries removed\n";
320 else
321 octave_stdout << "knobs(1): " << User_knobs (0)
322 << ", only completely dense rows removed\n";
323
324 if (knobs[COLAMD_DENSE_COL] >= 0)
325 octave_stdout << "knobs(2): " << User_knobs (1)
326 << ", cols with > max (16,"
327 << knobs[COLAMD_DENSE_COL] << "*sqrt (size(A)))"
328 << " entries removed\n";
329 else
330 octave_stdout << "knobs(2): " << User_knobs (1)
331 << ", only completely dense columns removed\n";
332
333 octave_stdout << "knobs(3): " << User_knobs (2)
334 << ", statistics and knobs printed\n";
335
336 }
337 }
338
339 octave_idx_type n_row, n_col, nnz;
340 octave_idx_type *ridx, *cidx;
341 SparseComplexMatrix scm;
342 SparseMatrix sm;
343
344 if (args(0).is_sparse_type ())
345 {
346 if (args(0).is_complex_type ())
347 {
348 scm = args(0). sparse_complex_matrix_value ();
349 n_row = scm.rows ();
350 n_col = scm.cols ();
351 nnz = scm.nnz ();
352 ridx = scm.xridx ();
353 cidx = scm.xcidx ();
354 }
355 else
356 {
357 sm = args(0).sparse_matrix_value ();
358
359 n_row = sm.rows ();
360 n_col = sm.cols ();
361 nnz = sm.nnz ();
362 ridx = sm.xridx ();
363 cidx = sm.xcidx ();
364 }
365 }
366 else
367 {
368 if (args(0).is_complex_type ())
369 sm = SparseMatrix (real (args(0).complex_matrix_value ()));
370 else
371 sm = SparseMatrix (args(0).matrix_value ());
372
373 n_row = sm.rows ();
374 n_col = sm.cols ();
375 nnz = sm.nnz ();
376 ridx = sm.xridx ();
377 cidx = sm.xcidx ();
378 }
379
380 // Allocate workspace for colamd
381 OCTAVE_LOCAL_BUFFER (octave_idx_type, p, n_col+1);
382 for (octave_idx_type i = 0; i < n_col+1; i++)
383 p[i] = cidx[i];
384
385 octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col);
386 OCTAVE_LOCAL_BUFFER (octave_idx_type, A, Alen);
387 for (octave_idx_type i = 0; i < nnz; i++)
388 A[i] = ridx[i];
389
390 // Order the columns (destroys A)
391 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
392 if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats))
393 {
394 COLAMD_NAME (_report) (stats) ;
395 error ("colamd: internal error!");
396 return retval;
397 }
398
399 // column elimination tree post-ordering (reuse variables)
400 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1);
401 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col + 1);
402 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
403
404 for (octave_idx_type i = 0; i < n_col; i++)
405 {
406 colbeg[i] = cidx[p[i]];
407 colend[i] = cidx[p[i]+1];
408 }
409
410 coletree (ridx, colbeg, colend, etree, n_row, n_col);
411
412 // Calculate the tree post-ordering
413 tree_postorder (n_col, etree, colbeg);
414
415 // return the permutation vector
416 NDArray out_perm (dim_vector (1, n_col));
417 for (octave_idx_type i = 0; i < n_col; i++)
418 out_perm(i) = p[colbeg[i]] + 1;
419
420 retval(0) = out_perm;
421
422 // print stats if spumoni > 0
423 if (spumoni > 0)
424 COLAMD_NAME (_report) (stats) ;
425
426 // Return the stats vector
427 if (nargout == 2)
428 {
429 NDArray out_stats (dim_vector (1, COLAMD_STATS));
430 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
431 out_stats(i) = stats[i] ;
432 retval(1) = out_stats;
433
434 // fix stats (5) and (6), for 1-based information on
435 // jumbled matrix. note that this correction doesn't
436 // occur if symamd returns FALSE
437 out_stats (COLAMD_INFO1) ++ ;
438 out_stats (COLAMD_INFO2) ++ ;
439 }
440 }
441
442 #else
443
444 error ("colamd: not available in this version of Octave");
445
446 #endif
447
448 return retval;
449 }
450
451 DEFUN_DLD (symamd, args, nargout,
452 "-*- texinfo -*-\n\
453 @deftypefn {Loadable Function} {@var{p} =} symamd (@var{S})\n\
454 @deftypefnx {Loadable Function} {@var{p} =} symamd (@var{S}, @var{knobs})\n\
455 @deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} symamd (@var{S})\n\
456 @deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} symamd (@var{S}, @var{knobs})\n\
457 \n\
458 For a symmetric positive definite matrix @var{S}, returns the permutation\n\
459 vector p such that @code{@var{S}(@var{p}, @var{p})} tends to have a\n\
460 sparser Cholesky@tie{}factor than @var{S}. Sometimes @code{symamd} works\n\
461 well for symmetric indefinite matrices too. The matrix @var{S} is assumed\n\
462 to be symmetric; only the strictly lower triangular part is referenced.\n\
463 @var{S} must be square.\n\
464 \n\
465 @var{knobs} is an optional one- to two-element input vector. If @var{S} is\n\
466 n-by-n, then rows and columns with more than\n\
467 @code{max (16,@var{knobs}(1)*sqrt(n))} entries are removed prior to ordering,\n\
468 and ordered last in the output permutation @var{p}. No rows/columns are\n\
469 removed if @code{@var{knobs}(1) < 0}. If @code{@var{knobs} (2)} is nonzero,\n\
470 @code{stats} and @var{knobs} are printed. The default is @code{@var{knobs}\n\
471 = [10 0]}. Note that @var{knobs} differs from earlier versions of symamd.\n\
472 \n\
473 @var{stats} is an optional 20-element output vector that provides data\n\
474 about the ordering and the validity of the input matrix @var{S}. Ordering\n\
475 statistics are in @code{@var{stats}(1:3)}. @code{@var{stats}(1) =\n\
476 @var{stats}(2)} is the number of dense or empty rows and columns\n\
477 ignored by SYMAMD and @code{@var{stats}(3)} is the number of garbage\n\
478 collections performed on the internal data structure used by SYMAMD\n\
479 (roughly of size @code{8.4 * nnz (tril (@var{S}, -1)) + 9 * @var{n}}\n\
480 integers).\n\
481 \n\
482 Octave built-in functions are intended to generate valid sparse matrices,\n\
483 with no duplicate entries, with ascending row indices of the nonzeros\n\
484 in each column, with a non-negative number of entries in each column (!)\n\
485 and so on. If a matrix is invalid, then SYMAMD may or may not be able\n\
486 to continue. If there are duplicate entries (a row index appears two or\n\
487 more times in the same column) or if the row indices in a column are out\n\
488 of order, then SYMAMD can correct these errors by ignoring the duplicate\n\
489 entries and sorting each column of its internal copy of the matrix S (the\n\
490 input matrix S is not repaired, however). If a matrix is invalid in\n\
491 other ways then SYMAMD cannot continue, an error message is printed, and\n\
492 no output arguments (@var{p} or @var{stats}) are returned. SYMAMD is\n\
493 thus a simple way to check a sparse matrix to see if it's valid.\n\
494 \n\
495 @code{@var{stats}(4:7)} provide information if SYMAMD was able to\n\
496 continue. The matrix is OK if @code{@var{stats} (4)} is zero, or 1\n\
497 if invalid. @code{@var{stats}(5)} is the rightmost column index that\n\
498 is unsorted or contains duplicate entries, or zero if no such column\n\
499 exists. @code{@var{stats}(6)} is the last seen duplicate or out-of-order\n\
500 row index in the column index given by @code{@var{stats}(5)}, or zero\n\
501 if no such row index exists. @code{@var{stats}(7)} is the number of\n\
502 duplicate or out-of-order row indices. @code{@var{stats}(8:20)} is\n\
503 always zero in the current version of SYMAMD (reserved for future use).\n\
504 \n\
505 The ordering is followed by a column elimination tree post-ordering.\n\
506 \n\
507 The authors of the code itself are Stefan I. Larimore and Timothy A.\n\
508 Davis @email{davis@@cise.ufl.edu}, University of Florida. The algorithm was\n\
509 developed in collaboration with John Gilbert, Xerox PARC, and Esmond\n\
510 Ng, Oak Ridge National Laboratory. (see\n\
511 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\
512 @seealso{colperm, colamd}\n\
513 @end deftypefn")
514 {
515 octave_value_list retval;
516
517 #ifdef HAVE_COLAMD
518
519 int nargin = args.length ();
520 int spumoni = 0;
521
522 if (nargout > 2 || nargin < 1 || nargin > 2)
523 print_usage ();
524 else
525 {
526 // Get knobs
527 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
528 COLAMD_NAME (_set_defaults) (knobs);
529
530 // Check for user-passed knobs
531 if (nargin == 2)
532 {
533 NDArray User_knobs = args(1).array_value ();
534 int nel_User_knobs = User_knobs.length ();
535
536 if (nel_User_knobs > 0)
537 knobs[COLAMD_DENSE_ROW] = User_knobs(COLAMD_DENSE_ROW);
538 if (nel_User_knobs > 1)
539 spumoni = static_cast<int> (User_knobs (1));
540 }
541
542 // print knob settings if spumoni is set
543 if (spumoni > 0)
544 octave_stdout << "symamd: dense row/col fraction: "
545 << knobs[COLAMD_DENSE_ROW] << std::endl;
546
547 octave_idx_type n_row, n_col;
548 octave_idx_type *ridx, *cidx;
549 SparseMatrix sm;
550 SparseComplexMatrix scm;
551
552 if (args(0).is_sparse_type ())
553 {
554 if (args(0).is_complex_type ())
555 {
556 scm = args(0).sparse_complex_matrix_value ();
557 n_row = scm.rows ();
558 n_col = scm.cols ();
559 ridx = scm.xridx ();
560 cidx = scm.xcidx ();
561 }
562 else
563 {
564 sm = args(0).sparse_matrix_value ();
565 n_row = sm.rows ();
566 n_col = sm.cols ();
567 ridx = sm.xridx ();
568 cidx = sm.xcidx ();
569 }
570 }
571 else
572 {
573 if (args(0).is_complex_type ())
574 sm = SparseMatrix (real (args(0).complex_matrix_value ()));
575 else
576 sm = SparseMatrix (args(0).matrix_value ());
577
578 n_row = sm.rows ();
579 n_col = sm.cols ();
580 ridx = sm.xridx ();
581 cidx = sm.xcidx ();
582 }
583
584 if (n_row != n_col)
585 {
586 error ("symamd: matrix S must be square");
587 return retval;
588 }
589
590 // Allocate workspace for symamd
591 OCTAVE_LOCAL_BUFFER (octave_idx_type, perm, n_col+1);
592 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
593 if (!SYMAMD_NAME () (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free))
594 {
595 SYMAMD_NAME (_report) (stats) ;
596 error ("symamd: internal error!") ;
597 return retval;
598 }
599
600 // column elimination tree post-ordering
601 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
602 symetree (ridx, cidx, etree, perm, n_col);
603
604 // Calculate the tree post-ordering
605 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
606 tree_postorder (n_col, etree, post);
607
608 // return the permutation vector
609 NDArray out_perm (dim_vector (1, n_col));
610 for (octave_idx_type i = 0; i < n_col; i++)
611 out_perm(i) = perm[post[i]] + 1;
612
613 retval(0) = out_perm;
614
615 // print stats if spumoni > 0
616 if (spumoni > 0)
617 SYMAMD_NAME (_report) (stats) ;
618
619 // Return the stats vector
620 if (nargout == 2)
621 {
622 NDArray out_stats (dim_vector (1, COLAMD_STATS));
623 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
624 out_stats(i) = stats[i] ;
625 retval(1) = out_stats;
626
627 // fix stats (5) and (6), for 1-based information on
628 // jumbled matrix. note that this correction doesn't
629 // occur if symamd returns FALSE
630 out_stats (COLAMD_INFO1) ++ ;
631 out_stats (COLAMD_INFO2) ++ ;
632 }
633 }
634
635 #else
636
637 error ("symamd: not available in this version of Octave");
638
639 #endif
640
641 return retval;
642 }
643
644 DEFUN_DLD (etree, args, nargout,
645 "-*- texinfo -*-\n\
646 @deftypefn {Loadable Function} {@var{p} =} etree (@var{S})\n\
647 @deftypefnx {Loadable Function} {@var{p} =} etree (@var{S}, @var{typ})\n\
648 @deftypefnx {Loadable Function} {[@var{p}, @var{q}] =} etree (@var{S}, @var{typ})\n\
649 \n\
650 Return the elimination tree for the matrix @var{S}. By default @var{S}\n\
651 is assumed to be symmetric and the symmetric elimination tree is\n\
652 returned. The argument @var{typ} controls whether a symmetric or\n\
653 column elimination tree is returned. Valid values of @var{typ} are\n\
654 \"sym\" or \"col\", for symmetric or column elimination tree respectively\n\
655 \n\
656 Called with a second argument, @code{etree} also returns the postorder\n\
657 permutations on the tree.\n\
658 @end deftypefn")
659 {
660 octave_value_list retval;
661
662 int nargin = args.length ();
663
664 if (nargout > 2 || nargin < 1 || nargin > 2)
665 print_usage ();
666 else
667 {
668 octave_idx_type n_row, n_col;
669 octave_idx_type *ridx, *cidx;
670 bool is_sym = true;
671 SparseMatrix sm;
672 SparseComplexMatrix scm;
673
674 if (args(0).is_sparse_type ())
675 {
676 if (args(0).is_complex_type ())
677 {
678 scm = args(0).sparse_complex_matrix_value ();
679 n_row = scm.rows ();
680 n_col = scm.cols ();
681 ridx = scm.xridx ();
682 cidx = scm.xcidx ();
683 }
684 else
685 {
686 sm = args(0).sparse_matrix_value ();
687 n_row = sm.rows ();
688 n_col = sm.cols ();
689 ridx = sm.xridx ();
690 cidx = sm.xcidx ();
691 }
692
693 }
694 else
695 {
696 error ("etree: S must be a sparse matrix");
697 return retval;
698 }
699
700 if (nargin == 2)
701 {
702 if (args(1).is_string ())
703 {
704 std::string str = args(1).string_value ();
705 if (str.find ("C") == 0 || str.find ("c") == 0)
706 is_sym = false;
707 }
708 else
709 {
710 error ("etree: TYP must be a string");
711 return retval;
712 }
713 }
714
715 // column elimination tree post-ordering (reuse variables)
716 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
717
718 if (is_sym)
719 {
720 if (n_row != n_col)
721 {
722 error ("etree: S is marked as symmetric, but is not square");
723 return retval;
724 }
725
726 symetree (ridx, cidx, etree, 0, n_col);
727 }
728 else
729 {
730 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
731 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col);
732
733 for (octave_idx_type i = 0; i < n_col; i++)
734 {
735 colbeg[i] = cidx[i];
736 colend[i] = cidx[i+1];
737 }
738
739 coletree (ridx, colbeg, colend, etree, n_row, n_col);
740 }
741
742 NDArray tree (dim_vector (1, n_col));
743 for (octave_idx_type i = 0; i < n_col; i++)
744 // We flag a root with n_col while Matlab does it with zero
745 // Convert for matlab compatiable output
746 if (etree[i] == n_col)
747 tree(i) = 0;
748 else
749 tree(i) = etree[i] + 1;
750
751 retval(0) = tree;
752
753 if (nargout == 2)
754 {
755 // Calculate the tree post-ordering
756 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
757 tree_postorder (n_col, etree, post);
758
759 NDArray postorder (dim_vector (1, n_col));
760 for (octave_idx_type i = 0; i < n_col; i++)
761 postorder(i) = post[i] + 1;
762
763 retval(1) = postorder;
764 }
765 }
766
767 return retval;
768 }