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