Mercurial > octave-nkf

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