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
colamd.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2017 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 #if defined (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 "errwarn.h"
39 #include "pager.h"
40 #include "ov-re-mat.h"
41 
42 #include "ov-re-sparse.h"
43 #include "ov-cx-sparse.h"
44 
45 #include "oct-sparse.h"
46 #include "oct-locbuf.h"
47 
48 #if defined (OCTAVE_ENABLE_64)
49 # define COLAMD_NAME(name) colamd_l ## name
50 # define SYMAMD_NAME(name) symamd_l ## name
51 #else
52 # define COLAMD_NAME(name) colamd ## name
53 # define SYMAMD_NAME(name) symamd ## name
54 #endif
55 
56 // The symmetric column elimination tree code take from the Davis LDL code.
57 // Copyright given elsewhere in this file.
58 static void
59 symetree (const octave_idx_type *ridx, const octave_idx_type *cidx,
61 {
63  OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
64  if (P)
65  // If P is present then compute Pinv, the inverse of P
66  for (octave_idx_type k = 0 ; k < n ; k++)
67  Pinv[P[k]] = k;
68 
69  for (octave_idx_type k = 0 ; k < n ; k++)
70  {
71  // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k)
72  Parent[k] = n ; // parent of k is not yet known
73  Flag[k] = k ; // mark node k as visited
74  octave_idx_type kk = (P) ? P[k] // kth original, or permuted, column
75  : (k);
76  octave_idx_type p2 = cidx[kk+1];
77  for (octave_idx_type p = cidx[kk] ; p < p2 ; p++)
78  {
79  // A (i,k) is nonzero (original or permuted A)
80  octave_idx_type i = (Pinv) ? (Pinv[ridx[p]]) : (ridx[p]);
81  if (i < k)
82  {
83  // follow path from i to root of etree, stop at flagged node
84  for ( ; Flag[i] != k ; i = Parent[i])
85  {
86  // find parent of i if not yet determined
87  if (Parent[i] == n)
88  Parent[i] = k;
89  Flag[i] = k ; // mark i as visited
90  }
91  }
92  }
93  }
94 }
95 
96 // The elimination tree post-ordering code below is taken from SuperLU
97 static inline octave_idx_type
99 {
100  pp[i] = i;
101  return i;
102 }
103 
104 static inline octave_idx_type
106 {
107  pp[s] = t;
108  return t;
109 }
110 
111 static inline octave_idx_type
113 {
114  octave_idx_type p = pp[i];
115  octave_idx_type 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
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
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,
168 {
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  doc: /* -*- texinfo -*-
213 @deftypefn {} {@var{p} =} colamd (@var{S})
214 @deftypefnx {} {@var{p} =} colamd (@var{S}, @var{knobs})
215 @deftypefnx {} {[@var{p}, @var{stats}] =} colamd (@var{S})
216 @deftypefnx {} {[@var{p}, @var{stats}] =} colamd (@var{S}, @var{knobs})
217 
218 Compute the column approximate minimum degree permutation.
219 
220 @code{@var{p} = colamd (@var{S})} returns the column approximate minimum
221 degree permutation vector for the sparse matrix @var{S}. For a
222 non-symmetric matrix @var{S}, @code{@var{S}(:,@var{p})} tends to have
223 sparser LU@tie{}factors than @var{S}. The Cholesky@tie{}factorization of
224 @code{@var{S}(:,@var{p})' * @var{S}(:,@var{p})} also tends to be sparser
225 than that of @code{@var{S}' * @var{S}}.
226 
227 @var{knobs} is an optional one- to three-element input vector. If @var{S}
228 is m-by-n, then rows with more than @code{max(16,@var{knobs}(1)*sqrt(n))}
229 entries are ignored. Columns with more than
230 @code{max (16,@var{knobs}(2)*sqrt(min(m,n)))} entries are removed prior to
231 ordering, and ordered last in the output permutation @var{p}. Only
232 completely dense rows or columns are removed if @code{@var{knobs}(1)} and
233 @code{@var{knobs}(2)} are < 0, respectively. If @code{@var{knobs}(3)} is
234 nonzero, @var{stats} and @var{knobs} are printed. The default is
235 @code{@var{knobs} = [10 10 0]}. Note that @var{knobs} differs from earlier
236 versions of colamd.
237 
238 @var{stats} is an optional 20-element output vector that provides data
239 about the ordering and the validity of the input matrix @var{S}. Ordering
240 statistics are in @code{@var{stats}(1:3)}. @code{@var{stats}(1)} and
241 @code{@var{stats}(2)} are the number of dense or empty rows and columns
242 ignored by @sc{colamd} and @code{@var{stats}(3)} is the number of garbage
243 collections performed on the internal data structure used by @sc{colamd}
244 (roughly of size @code{2.2 * nnz(@var{S}) + 4 * @var{m} + 7 * @var{n}}
245 integers).
246 
247 Octave built-in functions are intended to generate valid sparse matrices,
248 with no duplicate entries, with ascending row indices of the nonzeros
249 in each column, with a non-negative number of entries in each column (!)
250 and so on. If a matrix is invalid, then @sc{colamd} may or may not be able
251 to continue. If there are duplicate entries (a row index appears two or
252 more times in the same column) or if the row indices in a column are out
253 of order, then @sc{colamd} can correct these errors by ignoring the
254 duplicate entries and sorting each column of its internal copy of the
255 matrix @var{S} (the input matrix @var{S} is not repaired, however). If a
256 matrix is invalid in other ways then @sc{colamd} cannot continue, an error
257 message is printed, and no output arguments (@var{p} or @var{stats}) are
258 returned.
259 @sc{colamd} is thus a simple way to check a sparse matrix to see if it's
260 valid.
261 
262 @code{@var{stats}(4:7)} provide information if @sc{colamd} was able to
263 continue. The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if
264 invalid. @code{@var{stats}(5)} is the rightmost column index that is
265 unsorted or contains duplicate entries, or zero if no such column exists.
266 @code{@var{stats}(6)} is the last seen duplicate or out-of-order row
267 index in the column index given by @code{@var{stats}(5)}, or zero if no
268 such row index exists. @code{@var{stats}(7)} is the number of duplicate
269 or out-of-order row indices. @code{@var{stats}(8:20)} is always zero in
270 the current version of @sc{colamd} (reserved for future use).
271 
272 The ordering is followed by a column elimination tree post-ordering.
273 
274 The authors of the code itself are @nospell{Stefan I. Larimore} and
275 @nospell{Timothy A. Davis @email{davis@@cise.ufl.edu}}, University of Florida. The algorithm was developed in collaboration with @nospell{John Gilbert},
276 Xerox PARC, and @nospell{Esmond Ng}, Oak Ridge National Laboratory. (see
277 @url{http://www.cise.ufl.edu/research/sparse/colamd})
278 @seealso{colperm, symamd, ccolamd}
279 @end deftypefn */)
280 {
281 #if defined (HAVE_COLAMD)
282 
283  int nargin = args.length ();
284 
285  if (nargin < 1 || nargin > 2)
286  print_usage ();
287 
288  octave_value_list retval (nargout == 2 ? 2 : 1);
289  int spumoni = 0;
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.numel ();
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
314  << ", " << COLAMD_DATE << ":\n";
315 
316  if (knobs[COLAMD_DENSE_ROW] >= 0)
317  octave_stdout << "knobs(1): " << User_knobs (0)
318  << ", rows with > max (16,"
319  << knobs[COLAMD_DENSE_ROW] << "*sqrt (size(A,2)))"
320  << " entries removed\n";
321  else
322  octave_stdout << "knobs(1): " << User_knobs (0)
323  << ", only completely dense rows removed\n";
324 
325  if (knobs[COLAMD_DENSE_COL] >= 0)
326  octave_stdout << "knobs(2): " << User_knobs (1)
327  << ", cols with > max (16,"
328  << knobs[COLAMD_DENSE_COL] << "*sqrt (size(A)))"
329  << " entries removed\n";
330  else
331  octave_stdout << "knobs(2): " << User_knobs (1)
332  << ", only completely dense columns removed\n";
333 
334  octave_stdout << "knobs(3): " << User_knobs (2)
335  << ", statistics and knobs printed\n";
336 
337  }
338  }
339 
340  octave_idx_type n_row, n_col, nnz;
341  octave_idx_type *ridx, *cidx;
343  SparseMatrix sm;
344 
345  if (args(0).is_sparse_type ())
346  {
347  if (args(0).is_complex_type ())
348  {
349  scm = args(0).sparse_complex_matrix_value ();
350  n_row = scm.rows ();
351  n_col = scm.cols ();
352  nnz = scm.nnz ();
353  ridx = scm.xridx ();
354  cidx = scm.xcidx ();
355  }
356  else
357  {
358  sm = args(0).sparse_matrix_value ();
359 
360  n_row = sm.rows ();
361  n_col = sm.cols ();
362  nnz = sm.nnz ();
363  ridx = sm.xridx ();
364  cidx = sm.xcidx ();
365  }
366  }
367  else
368  {
369  if (args(0).is_complex_type ())
370  sm = SparseMatrix (real (args(0).complex_matrix_value ()));
371  else
372  sm = SparseMatrix (args(0).matrix_value ());
373 
374  n_row = sm.rows ();
375  n_col = sm.cols ();
376  nnz = sm.nnz ();
377  ridx = sm.xridx ();
378  cidx = sm.xcidx ();
379  }
380 
381  // Allocate workspace for colamd
383  for (octave_idx_type i = 0; i < n_col+1; i++)
384  p[i] = cidx[i];
385 
386  octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col);
388  for (octave_idx_type i = 0; i < nnz; i++)
389  A[i] = ridx[i];
390 
391  // Order the columns (destroys A)
392  OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
393  if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats))
394  {
395  COLAMD_NAME (_report)(stats);
396 
397  error ("colamd: internal error!");
398  }
399 
400  // column elimination tree post-ordering (reuse variables)
401  OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1);
402  OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col + 1);
403  OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
404 
405  for (octave_idx_type i = 0; i < n_col; i++)
406  {
407  colbeg[i] = cidx[p[i]];
408  colend[i] = cidx[p[i]+1];
409  }
410 
411  coletree (ridx, colbeg, colend, etree, n_row, n_col);
412 
413  // Calculate the tree post-ordering
414  tree_postorder (n_col, etree, colbeg);
415 
416  // return the permutation vector
417  NDArray out_perm (dim_vector (1, n_col));
418  for (octave_idx_type i = 0; i < n_col; i++)
419  out_perm(i) = p[colbeg[i]] + 1;
420 
421  retval(0) = out_perm;
422 
423  // print stats if spumoni > 0
424  if (spumoni > 0)
425  COLAMD_NAME (_report)(stats);
426 
427  // Return the stats vector
428  if (nargout == 2)
429  {
430  NDArray out_stats (dim_vector (1, COLAMD_STATS));
431  for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
432  out_stats(i) = stats[i];
433  retval(1) = out_stats;
434 
435  // fix stats (5) and (6), for 1-based information on
436  // jumbled matrix. note that this correction doesn't
437  // occur if symamd returns FALSE
438  out_stats(COLAMD_INFO1)++;
439  out_stats(COLAMD_INFO2)++;
440  }
441 
442  return retval;
443 
444 #else
445 
446  octave_unused_parameter (args);
447  octave_unused_parameter (nargout);
448 
449  err_disabled_feature ("colamd", "COLAMD");
450 
451 #endif
452 }
453 
454 DEFUN_DLD (symamd, args, nargout,
455  doc: /* -*- texinfo -*-
456 @deftypefn {} {@var{p} =} symamd (@var{S})
457 @deftypefnx {} {@var{p} =} symamd (@var{S}, @var{knobs})
458 @deftypefnx {} {[@var{p}, @var{stats}] =} symamd (@var{S})
459 @deftypefnx {} {[@var{p}, @var{stats}] =} symamd (@var{S}, @var{knobs})
460 
461 For a symmetric positive definite matrix @var{S}, returns the permutation
462 vector p such that @code{@var{S}(@var{p}, @var{p})} tends to have a
463 sparser Cholesky@tie{}factor than @var{S}.
464 
465 Sometimes @code{symamd} works well for symmetric indefinite matrices too.
466 The matrix @var{S} is assumed to be symmetric; only the strictly lower
467 triangular part is referenced. @var{S} must be square.
468 
469 @var{knobs} is an optional one- to two-element input vector. If @var{S} is
470 n-by-n, then rows and columns with more than
471 @code{max (16,@var{knobs}(1)*sqrt(n))} entries are removed prior to
472 ordering, and ordered last in the output permutation @var{p}. No
473 rows/columns are removed if @code{@var{knobs}(1) < 0}. If
474 @code{@var{knobs} (2)} is nonzero, @code{stats} and @var{knobs} are
475 printed. The default is @code{@var{knobs} = [10 0]}. Note that
476 @var{knobs} differs from earlier versions of @code{symamd}.
477 
478 @var{stats} is an optional 20-element output vector that provides data
479 about the ordering and the validity of the input matrix @var{S}. Ordering
480 statistics are in @code{@var{stats}(1:3)}.
481 @code{@var{stats}(1) = @var{stats}(2)} is the number of dense or empty rows
482 and columns ignored by SYMAMD and @code{@var{stats}(3)} is the number of
483 garbage collections performed on the internal data structure used by SYMAMD
484 (roughly of size @code{8.4 * nnz (tril (@var{S}, -1)) + 9 * @var{n}}
485 integers).
486 
487 Octave built-in functions are intended to generate valid sparse matrices,
488 with no duplicate entries, with ascending row indices of the nonzeros
489 in each column, with a non-negative number of entries in each column (!)
490 and so on. If a matrix is invalid, then SYMAMD may or may not be able
491 to continue. If there are duplicate entries (a row index appears two or
492 more times in the same column) or if the row indices in a column are out
493 of order, then SYMAMD can correct these errors by ignoring the duplicate
494 entries and sorting each column of its internal copy of the matrix S (the
495 input matrix S is not repaired, however). If a matrix is invalid in
496 other ways then SYMAMD cannot continue, an error message is printed, and
497 no output arguments (@var{p} or @var{stats}) are returned. SYMAMD is
498 thus a simple way to check a sparse matrix to see if it's valid.
499 
500 @code{@var{stats}(4:7)} provide information if SYMAMD was able to
501 continue. The matrix is OK if @code{@var{stats} (4)} is zero, or 1
502 if invalid. @code{@var{stats}(5)} is the rightmost column index that
503 is unsorted or contains duplicate entries, or zero if no such column
504 exists. @code{@var{stats}(6)} is the last seen duplicate or out-of-order
505 row index in the column index given by @code{@var{stats}(5)}, or zero
506 if no such row index exists. @code{@var{stats}(7)} is the number of
507 duplicate or out-of-order row indices. @code{@var{stats}(8:20)} is
508 always zero in the current version of SYMAMD (reserved for future use).
509 
510 The ordering is followed by a column elimination tree post-ordering.
511 
512 The authors of the code itself are @nospell{Stefan I. Larimore} and
513 @nospell{Timothy A. Davis @email{davis@@cise.ufl.edu}}, University of Florida. The algorithm was developed in collaboration with @nospell{John Gilbert},
514 Xerox PARC, and @nospell{Esmond Ng}, Oak Ridge National Laboratory. (see
515 @url{http://www.cise.ufl.edu/research/sparse/colamd})
516 @seealso{colperm, colamd}
517 @end deftypefn */)
518 {
519 #if defined (HAVE_COLAMD)
520 
521  int nargin = args.length ();
522 
523  if (nargin < 1 || nargin > 2)
524  print_usage ();
525 
526  octave_value_list retval (nargin == 2 ? 2 : 1);
527  int spumoni = 0;
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.numel ();
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  err_square_matrix_required ("symamd", "S");
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,
594  knobs, stats, &calloc, &free))
595  {
596  SYMAMD_NAME (_report)(stats);
597 
598  error ("symamd: internal error!");
599  }
600 
601  // column elimination tree post-ordering
602  OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
603  symetree (ridx, cidx, etree, perm, n_col);
604 
605  // Calculate the tree post-ordering
606  OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
607  tree_postorder (n_col, etree, post);
608 
609  // return the permutation vector
610  NDArray out_perm (dim_vector (1, n_col));
611  for (octave_idx_type i = 0; i < n_col; i++)
612  out_perm(i) = perm[post[i]] + 1;
613 
614  retval(0) = out_perm;
615 
616  // print stats if spumoni > 0
617  if (spumoni > 0)
618  SYMAMD_NAME (_report)(stats);
619 
620  // Return the stats vector
621  if (nargout == 2)
622  {
623  NDArray out_stats (dim_vector (1, COLAMD_STATS));
624  for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
625  out_stats(i) = stats[i];
626  retval(1) = out_stats;
627 
628  // fix stats (5) and (6), for 1-based information on
629  // jumbled matrix. note that this correction doesn't
630  // occur if symamd returns FALSE
631  out_stats(COLAMD_INFO1)++;
632  out_stats(COLAMD_INFO2)++;
633  }
634 
635  return retval;
636 
637 #else
638 
639  octave_unused_parameter (args);
640  octave_unused_parameter (nargout);
641 
642  err_disabled_feature ("symamd", "COLAMD");
643 
644 #endif
645 }
646 
647 DEFUN_DLD (etree, args, nargout,
648  doc: /* -*- texinfo -*-
649 @deftypefn {} {@var{p} =} etree (@var{S})
650 @deftypefnx {} {@var{p} =} etree (@var{S}, @var{typ})
651 @deftypefnx {} {[@var{p}, @var{q}] =} etree (@var{S}, @var{typ})
652 
653 Return the elimination tree for the matrix @var{S}.
654 
655 By default @var{S} is assumed to be symmetric and the symmetric elimination
656 tree is returned. The argument @var{typ} controls whether a symmetric or
657 column elimination tree is returned. Valid values of @var{typ} are
658 @qcode{"sym"} or @qcode{"col"}, for symmetric or column elimination tree
659 respectively.
660 
661 Called with a second argument, @code{etree} also returns the postorder
662 permutations on the tree.
663 @end deftypefn */)
664 {
665  int nargin = args.length ();
666 
667  if (nargin < 1 || nargin > 2)
668  print_usage ();
669 
670  octave_value_list retval (nargout == 2 ? 2 : 1);
671 
672  octave_idx_type n_row = 0;
673  octave_idx_type n_col = 0;
674  octave_idx_type *ridx = 0;
675  octave_idx_type *cidx = 0;
676 
677  if (args(0).is_sparse_type ())
678  error ("etree: S must be a sparse matrix");
679 
680  if (args(0).is_complex_type ())
681  {
682  SparseComplexMatrix scm = args(0).sparse_complex_matrix_value ();
683 
684  n_row = scm.rows ();
685  n_col = scm.cols ();
686  ridx = scm.xridx ();
687  cidx = scm.xcidx ();
688  }
689  else
690  {
691  SparseMatrix sm = args(0).sparse_matrix_value ();
692 
693  n_row = sm.rows ();
694  n_col = sm.cols ();
695  ridx = sm.xridx ();
696  cidx = sm.xcidx ();
697  }
698 
699  bool is_sym = true;
700 
701  if (nargin == 2)
702  {
703  std::string str = args(1).xstring_value ("etree: TYP must be a string");
704  if (str.find ("C") == 0 || str.find ("c") == 0)
705  is_sym = false;
706  }
707 
708  // column elimination tree post-ordering (reuse variables)
709  OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
710 
711  if (is_sym)
712  {
713  if (n_row != n_col)
714  error ("etree: S is marked as symmetric, but is not square");
715 
716  symetree (ridx, cidx, etree, 0, n_col);
717  }
718  else
719  {
720  OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
721  OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col);
722 
723  for (octave_idx_type i = 0; i < n_col; i++)
724  {
725  colbeg[i] = cidx[i];
726  colend[i] = cidx[i+1];
727  }
728 
729  coletree (ridx, colbeg, colend, etree, n_row, n_col);
730  }
731 
732  NDArray tree (dim_vector (1, n_col));
733  for (octave_idx_type i = 0; i < n_col; i++)
734  // We flag a root with n_col while Matlab does it with zero
735  // Convert for matlab compatiable output
736  if (etree[i] == n_col)
737  tree(i) = 0;
738  else
739  tree(i) = etree[i] + 1;
740 
741  retval(0) = tree;
742 
743  if (nargout == 2)
744  {
745  // Calculate the tree post-ordering
746  OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
747  tree_postorder (n_col, etree, post);
748 
749  NDArray postorder (dim_vector (1, n_col));
750  for (octave_idx_type i = 0; i < n_col; i++)
751  postorder(i) = post[i] + 1;
752 
753  retval(1) = postorder;
754  }
755 
756  return retval;
757 }
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
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
for large enough k
Definition: lu.cc:606
octave_idx_type * xcidx(void)
Definition: Sparse.h:549
void error(const char *fmt,...)
Definition: error.cc:570
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function t
Definition: ov-usr-fcn.cc:935
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:112
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:165
s
Definition: file-io.cc:2682
static void tree_postorder(octave_idx_type n, octave_idx_type *parent, octave_idx_type *post)
Definition: colamd.cc:142
octave_idx_type nnz(void) const
Actual number of nonzero terms.
Definition: Sparse.h:253
static octave_idx_type make_set(octave_idx_type i, octave_idx_type *pp)
Definition: colamd.cc:98
static octave_idx_type link(octave_idx_type s, octave_idx_type t, octave_idx_type *pp)
Definition: colamd.cc:105
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
std::complex< double > w(std::complex< double > z, double relerr=0)
int nargin
Definition: graphics.cc:10115
std::string str
Definition: hash.cc:118
octave_value retval
Definition: data.cc:6294
static octave_idx_type find(octave_idx_type i, octave_idx_type *pp)
Definition: colamd.cc:112
#define octave_stdout
Definition: pager.h:146
void free(void *)
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
p
Definition: lu.cc:138
Definition: pt.h:39
#define SYMAMD_NAME(name)
Definition: colamd.cc:53
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
#define DEFUN_DLD(name, args_name, nargout_name, doc)
Definition: defun-dld.h:45
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:59
#define COLAMD_NAME(name)
Definition: colamd.cc:52
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:136
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:129
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
void err_disabled_feature(const std::string &fcn, const std::string &feature, const std::string &pkg)
Definition: errwarn.cc:50
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:854
where the brackets indicate optional arguments and and character or cell array For character arrays the conversion is repeated for every row
Definition: str2double.cc:342
F77_RET_T const F77_INT F77_CMPLX * A