Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/colamd.cc @ 8920:eb63fbe60fab
update copyright notices
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Sat, 07 Mar 2009 10:41:27 -0500 |
parents | 25bc2d31e1bf |
children | 7c02ec148a3c |
rev | line source |
---|---|
5164 | 1 /* |
2 | |
8920 | 3 Copyright (C) 2004, 2005, 2006, 2007, 2008 David Bateman |
7016 | 4 Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 Andy Adler |
5 | |
6 This file is part of Octave. | |
5164 | 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 | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
5164 | 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 | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
5164 | 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 | |
5451 | 44 #include "oct-sparse.h" |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7924
diff
changeset
|
45 #include "oct-locbuf.h" |
5164 | 46 |
5440 | 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 | |
5164 | 55 // The symmetric column elimination tree code take from the Davis LDL code. |
56 // Copyright given elsewhere in this file. | |
7924 | 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) | |
5164 | 60 { |
5440 | 61 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n); |
62 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0)); | |
5164 | 63 if (P) |
64 // If P is present then compute Pinv, the inverse of P | |
5440 | 65 for (octave_idx_type k = 0 ; k < n ; k++) |
5164 | 66 Pinv [P [k]] = k ; |
67 | |
5440 | 68 for (octave_idx_type k = 0 ; k < n ; k++) |
5164 | 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 | |
5440 | 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++) | |
5164 | 76 { |
77 // A (i,k) is nonzero (original or permuted A) | |
5440 | 78 octave_idx_type i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ; |
5164 | 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 | |
7924 | 95 static inline octave_idx_type |
96 make_set (octave_idx_type i, octave_idx_type *pp) | |
5164 | 97 { |
98 pp[i] = i; | |
99 return i; | |
100 } | |
101 | |
7924 | 102 static inline octave_idx_type |
103 link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp) | |
5164 | 104 { |
105 pp[s] = t; | |
106 return t; | |
107 } | |
108 | |
7924 | 109 static inline octave_idx_type |
110 find (octave_idx_type i, octave_idx_type *pp) | |
5164 | 111 { |
5440 | 112 register octave_idx_type p, gp; |
5164 | 113 |
114 p = pp[i]; | |
115 gp = pp[p]; | |
7924 | 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; | |
5164 | 126 } |
127 | |
7924 | 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) | |
5164 | 132 { |
7924 | 133 for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w]) |
5164 | 134 postnum = etdfs (w, first_kid, next_kid, post, postnum); |
7924 | 135 |
5164 | 136 post[postnum++] = v; |
137 | |
138 return postnum; | |
139 } | |
140 | |
7924 | 141 static void |
142 tree_postorder (octave_idx_type n, octave_idx_type *parent, | |
143 octave_idx_type *post) | |
5164 | 144 { |
145 // Allocate storage for working arrays and results | |
5440 | 146 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1); |
147 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1); | |
5164 | 148 |
149 // Set up structure describing children | |
7924 | 150 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1) |
151 /* do nothing */; | |
152 | |
5440 | 153 for (octave_idx_type v = n-1; v >= 0; v--) |
5164 | 154 { |
5440 | 155 octave_idx_type dad = parent[v]; |
5164 | 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 | |
7924 | 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) | |
5164 | 168 { |
5440 | 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); | |
5164 | 172 |
173 // Compute firstcol[row] = first nonzero column in row | |
7924 | 174 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc) |
175 /* do nothing */; | |
176 | |
5440 | 177 for (octave_idx_type col = 0; col < nc; col++) |
178 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) | |
5164 | 179 { |
5440 | 180 octave_idx_type row = ridx[p]; |
5164 | 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. | |
5440 | 189 for (octave_idx_type col = 0; col < nc; col++) |
5164 | 190 { |
5440 | 191 octave_idx_type cset = make_set (col, pp); |
5164 | 192 root[cset] = col; |
193 parent[col] = nc; | |
5440 | 194 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) |
5164 | 195 { |
5440 | 196 octave_idx_type row = firstcol[ridx[p]]; |
5164 | 197 if (row >= col) |
198 continue; | |
5440 | 199 octave_idx_type rset = find (row, pp); |
200 octave_idx_type rroot = root[rset]; | |
5164 | 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. @code{@var{p} = colamd\n\ | |
219 (@var{s})} returns the column approximate minimum degree permutation\n\ | |
220 vector for the sparse matrix @var{s}. For a non-symmetric matrix @var{s},\n\ | |
221 @code{@var{s} (:,@var{p})} tends to have sparser LU factors than @var{s}.\n\ | |
222 The Cholesky factorization of @code{@var{s} (:,@var{p})' * @var{s}\n\ | |
223 (:,@var{p})} also tends to be sparser than that of @code{@var{s}' *\n\ | |
224 @var{s}}.\n\ | |
225 \n\ | |
5440 | 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))} entries\n\ | |
228 are ignored. Columns with more than @code{max(16,knobs(2)*sqrt(min(m,n)))}\n\ | |
229 entries are removed prior to ordering, and ordered last in the output\n\ | |
230 permutation @var{p}. Only completely dense rows or columns are removed\n\ | |
231 if @code{@var{knobs} (1)} and @code{@var{knobs} (2)} are < 0, respectively.\n\ | |
232 If @code{@var{knobs} (3)} is nonzero, @var{stats} and @var{knobs} are\n\ | |
233 printed. The default is @code{@var{knobs} = [10 10 0]}. Note that\n\ | |
234 @var{knobs} differs from earlier versions of colamd\n\ | |
5164 | 235 \n\ |
236 @var{stats} is an optional 20-element output vector that provides data\n\ | |
237 about the ordering and the validity of the input matrix @var{s}. Ordering\n\ | |
238 statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1)} and\n\ | |
239 @code{@var{stats} (2)} are the number of dense or empty rows and columns\n\ | |
240 ignored by COLAMD and @code{@var{stats} (3)} is the number of garbage\n\ | |
241 collections performed on the internal data structure used by COLAMD\n\ | |
242 (roughly of size @code{2.2 * nnz(@var{s}) + 4 * @var{m} + 7 * @var{n}}\n\ | |
243 integers).\n\ | |
244 \n\ | |
245 Octave built-in functions are intended to generate valid sparse matrices,\n\ | |
246 with no duplicate entries, with ascending row indices of the nonzeros\n\ | |
247 in each column, with a non-negative number of entries in each column (!)\n\ | |
248 and so on. If a matrix is invalid, then COLAMD may or may not be able\n\ | |
249 to continue. If there are duplicate entries (a row index appears two or\n\ | |
250 more times in the same column) or if the row indices in a column are out\n\ | |
251 of order, then COLAMD can correct these errors by ignoring the duplicate\n\ | |
252 entries and sorting each column of its internal copy of the matrix\n\ | |
253 @var{s} (the input matrix @var{s} is not repaired, however). If a matrix\n\ | |
254 is invalid in other ways then COLAMD cannot continue, an error message is\n\ | |
255 printed, and no output arguments (@var{p} or @var{stats}) are returned.\n\ | |
256 COLAMD is thus a simple way to check a sparse matrix to see if it's\n\ | |
257 valid.\n\ | |
258 \n\ | |
259 @code{@var{stats} (4:7)} provide information if COLAMD was able to\n\ | |
260 continue. The matrix is OK if @code{@var{stats} (4)} is zero, or 1 if\n\ | |
261 invalid. @code{@var{stats} (5)} is the rightmost column index that is\n\ | |
262 unsorted or contains duplicate entries, or zero if no such column exists.\n\ | |
263 @code{@var{stats} (6)} is the last seen duplicate or out-of-order row\n\ | |
264 index in the column index given by @code{@var{stats} (5)}, or zero if no\n\ | |
265 such row index exists. @code{@var{stats} (7)} is the number of duplicate\n\ | |
266 or out-of-order row indices. @code{@var{stats} (8:20)} is always zero in\n\ | |
267 the current version of COLAMD (reserved for future use).\n\ | |
268 \n\ | |
269 The ordering is followed by a column elimination tree post-ordering.\n\ | |
270 \n\ | |
271 The authors of the code itself are Stefan I. Larimore and Timothy A.\n\ | |
272 Davis (davis@@cise.ufl.edu), University of Florida. The algorithm was\n\ | |
273 developed in collaboration with John Gilbert, Xerox PARC, and Esmond\n\ | |
274 Ng, Oak Ridge National Laboratory. (see\n\ | |
275 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\ | |
5642 | 276 @seealso{colperm, symamd}\n\ |
277 @end deftypefn") | |
5164 | 278 { |
279 octave_value_list retval; | |
5451 | 280 |
281 #ifdef HAVE_COLAMD | |
282 | |
5164 | 283 int nargin = args.length (); |
284 int spumoni = 0; | |
285 | |
286 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2) | |
5823 | 287 print_usage (); |
5164 | 288 else |
289 { | |
290 // Get knobs | |
5440 | 291 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); |
292 COLAMD_NAME (_set_defaults) (knobs); | |
5164 | 293 |
294 // Check for user-passed knobs | |
295 if (nargin == 2) | |
296 { | |
297 NDArray User_knobs = args(1).array_value (); | |
298 int nel_User_knobs = User_knobs.length (); | |
299 | |
300 if (nel_User_knobs > 0) | |
5440 | 301 knobs [COLAMD_DENSE_ROW] = User_knobs (0); |
5164 | 302 if (nel_User_knobs > 1) |
5440 | 303 knobs [COLAMD_DENSE_COL] = User_knobs (1) ; |
5164 | 304 if (nel_User_knobs > 2) |
5760 | 305 spumoni = static_cast<int> (User_knobs (2)); |
5440 | 306 |
307 // print knob settings if spumoni is set | |
308 if (spumoni) | |
309 { | |
310 | |
311 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "." | |
312 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n"; | |
5164 | 313 |
5440 | 314 if (knobs [COLAMD_DENSE_ROW] >= 0) |
315 octave_stdout << "knobs(1): " << User_knobs (0) | |
316 << ", rows with > max(16," | |
317 << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))" | |
318 << " entries removed\n"; | |
319 else | |
320 octave_stdout << "knobs(1): " << User_knobs (0) | |
321 << ", only completely dense rows removed\n"; | |
322 | |
323 if (knobs [COLAMD_DENSE_COL] >= 0) | |
324 octave_stdout << "knobs(2): " << User_knobs (1) | |
325 << ", cols with > max(16," | |
326 << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))" | |
327 << " entries removed\n"; | |
328 else | |
329 octave_stdout << "knobs(2): " << User_knobs (1) | |
330 << ", only completely dense columns removed\n"; | |
331 | |
332 octave_stdout << "knobs(3): " << User_knobs (2) | |
333 << ", statistics and knobs printed\n"; | |
334 | |
335 } | |
5164 | 336 } |
337 | |
5440 | 338 octave_idx_type n_row, n_col, nnz; |
339 octave_idx_type *ridx, *cidx; | |
5164 | 340 SparseComplexMatrix scm; |
341 SparseMatrix sm; | |
342 | |
5631 | 343 if (args(0).is_sparse_type ()) |
5164 | 344 { |
345 if (args(0).is_complex_type ()) | |
346 { | |
347 scm = args(0). sparse_complex_matrix_value (); | |
348 n_row = scm.rows (); | |
349 n_col = scm.cols (); | |
5604 | 350 nnz = scm.nzmax (); |
5164 | 351 ridx = scm.xridx (); |
352 cidx = scm.xcidx (); | |
353 } | |
354 else | |
355 { | |
356 sm = args(0).sparse_matrix_value (); | |
357 | |
358 n_row = sm.rows (); | |
359 n_col = sm.cols (); | |
5604 | 360 nnz = sm.nzmax (); |
5164 | 361 ridx = sm.xridx (); |
362 cidx = sm.xcidx (); | |
363 } | |
364 } | |
365 else | |
366 { | |
367 if (args(0).is_complex_type ()) | |
368 sm = SparseMatrix (real (args(0).complex_matrix_value ())); | |
369 else | |
370 sm = SparseMatrix (args(0).matrix_value ()); | |
371 | |
372 n_row = sm.rows (); | |
373 n_col = sm.cols (); | |
5604 | 374 nnz = sm.nzmax (); |
5164 | 375 ridx = sm.xridx (); |
376 cidx = sm.xcidx (); | |
377 } | |
378 | |
379 // Allocate workspace for colamd | |
5440 | 380 OCTAVE_LOCAL_BUFFER (octave_idx_type, p, n_col+1); |
381 for (octave_idx_type i = 0; i < n_col+1; i++) | |
5164 | 382 p[i] = cidx [i]; |
383 | |
5440 | 384 octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col); |
385 OCTAVE_LOCAL_BUFFER (octave_idx_type, A, Alen); | |
386 for (octave_idx_type i = 0; i < nnz; i++) | |
5164 | 387 A[i] = ridx [i]; |
388 | |
389 // Order the columns (destroys A) | |
5440 | 390 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS); |
391 if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats)) | |
5164 | 392 { |
5440 | 393 COLAMD_NAME (_report) (stats) ; |
5164 | 394 error ("colamd: internal error!"); |
395 return retval; | |
396 } | |
397 | |
398 // column elimination tree post-ordering (reuse variables) | |
5440 | 399 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1); |
400 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col + 1); | |
401 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); | |
5164 | 402 |
5440 | 403 for (octave_idx_type i = 0; i < n_col; i++) |
5164 | 404 { |
405 colbeg[i] = cidx[p[i]]; | |
406 colend[i] = cidx[p[i]+1]; | |
407 } | |
408 | |
409 coletree (ridx, colbeg, colend, etree, n_row, n_col); | |
410 | |
411 // Calculate the tree post-ordering | |
7924 | 412 tree_postorder (n_col, etree, colbeg); |
5164 | 413 |
414 // return the permutation vector | |
415 NDArray out_perm (dim_vector (1, n_col)); | |
5440 | 416 for (octave_idx_type i = 0; i < n_col; i++) |
5164 | 417 out_perm(i) = p [colbeg [i]] + 1; |
418 | |
7924 | 419 retval(0) = out_perm; |
5164 | 420 |
421 // print stats if spumoni > 0 | |
422 if (spumoni > 0) | |
5440 | 423 COLAMD_NAME (_report) (stats) ; |
5164 | 424 |
425 // Return the stats vector | |
426 if (nargout == 2) | |
427 { | |
428 NDArray out_stats (dim_vector (1, COLAMD_STATS)); | |
5440 | 429 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++) |
5164 | 430 out_stats (i) = stats [i] ; |
431 retval(1) = out_stats; | |
432 | |
433 // fix stats (5) and (6), for 1-based information on | |
434 // jumbled matrix. note that this correction doesn't | |
435 // occur if symamd returns FALSE | |
436 out_stats (COLAMD_INFO1) ++ ; | |
437 out_stats (COLAMD_INFO2) ++ ; | |
438 } | |
439 } | |
440 | |
5451 | 441 #else |
442 | |
443 error ("colamd: not available in this version of Octave"); | |
444 | |
445 #endif | |
446 | |
5164 | 447 return retval; |
448 } | |
449 | |
450 DEFUN_DLD (symamd, args, nargout, | |
451 "-*- texinfo -*-\n\ | |
452 @deftypefn {Loadable Function} {@var{p} =} symamd (@var{s})\n\ | |
453 @deftypefnx {Loadable Function} {@var{p} =} symamd (@var{s}, @var{knobs})\n\ | |
454 @deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} symamd (@var{s})\n\ | |
455 @deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} symamd (@var{s}, @var{knobs})\n\ | |
456 \n\ | |
457 For a symmetric positive definite matrix @var{s}, returns the permutation\n\ | |
458 vector p such that @code{@var{s} (@var{p}, @var{p})} tends to have a\n\ | |
459 sparser Cholesky factor than @var{s}. Sometimes SYMAMD works well for\n\ | |
460 symmetric indefinite matrices too. The matrix @var{s} is assumed to be\n\ | |
461 symmetric; only the strictly lower triangular part is referenced. @var{s}\n\ | |
462 must be square.\n\ | |
463 \n\ | |
5440 | 464 @var{knobs} is an optional one- to two-element input vector. If @var{s} is\n\ |
465 n-by-n, then rows and columns with more than\n\ | |
466 @code{max(16,@var{knobs}(1)*sqrt(n))} entries are removed prior to ordering,\n\ | |
467 and ordered last in the output permutation @var{p}. No rows/columns are\n\ | |
468 removed if @code{@var{knobs}(1) < 0}. If @code{@var{knobs} (2)} is nonzero,\n\ | |
469 @code{stats} and @var{knobs} are printed. The default is @code{@var{knobs} \n\ | |
470 = [10 0]}. Note that @var{knobs} differs from earlier versions of symamd.\n\ | |
5164 | 471 \n\ |
472 @var{stats} is an optional 20-element output vector that provides data\n\ | |
473 about the ordering and the validity of the input matrix @var{s}. Ordering\n\ | |
474 statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1) =\n\ | |
475 @var{stats} (2)} is the number of dense or empty rows and columns\n\ | |
476 ignored by SYMAMD and @code{@var{stats} (3)} is the number of garbage\n\ | |
477 collections performed on the internal data structure used by SYMAMD\n\ | |
478 (roughly of size @code{8.4 * nnz (tril (@var{s}, -1)) + 9 * @var{n}}\n\ | |
479 integers).\n\ | |
480 \n\ | |
481 Octave built-in functions are intended to generate valid sparse matrices,\n\ | |
482 with no duplicate entries, with ascending row indices of the nonzeros\n\ | |
483 in each column, with a non-negative number of entries in each column (!)\n\ | |
484 and so on. If a matrix is invalid, then SYMAMD may or may not be able\n\ | |
485 to continue. If there are duplicate entries (a row index appears two or\n\ | |
486 more times in the same column) or if the row indices in a column are out\n\ | |
487 of order, then SYMAMD can correct these errors by ignoring the duplicate\n\ | |
488 entries and sorting each column of its internal copy of the matrix S (the\n\ | |
489 input matrix S is not repaired, however). If a matrix is invalid in\n\ | |
490 other ways then SYMAMD cannot continue, an error message is printed, and\n\ | |
491 no output arguments (@var{p} or @var{stats}) are returned. SYMAMD is\n\ | |
492 thus a simple way to check a sparse matrix to see if it's valid.\n\ | |
493 \n\ | |
494 @code{@var{stats} (4:7)} provide information if SYMAMD was able to\n\ | |
495 continue. The matrix is OK if @code{@var{stats} (4)} is zero, or 1\n\ | |
496 if invalid. @code{@var{stats} (5)} is the rightmost column index that\n\ | |
497 is unsorted or contains duplicate entries, or zero if no such column\n\ | |
498 exists. @code{@var{stats} (6)} is the last seen duplicate or out-of-order\n\ | |
499 row index in the column index given by @code{@var{stats} (5)}, or zero\n\ | |
500 if no such row index exists. @code{@var{stats} (7)} is the number of\n\ | |
501 duplicate or out-of-order row indices. @code{@var{stats} (8:20)} is\n\ | |
502 always zero in the current version of SYMAMD (reserved for future use).\n\ | |
503 \n\ | |
504 The ordering is followed by a column elimination tree post-ordering.\n\ | |
505 \n\ | |
506 \n\ | |
507 The authors of the code itself are Stefan I. Larimore and Timothy A.\n\ | |
508 Davis (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\ | |
5642 | 512 @seealso{colperm, colamd}\n\ |
513 @end deftypefn") | |
5164 | 514 { |
515 octave_value_list retval; | |
5451 | 516 |
517 #ifdef HAVE_COLAMD | |
518 | |
5164 | 519 int nargin = args.length (); |
520 int spumoni = 0; | |
521 | |
522 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2) | |
5823 | 523 print_usage (); |
5164 | 524 else |
525 { | |
526 // Get knobs | |
527 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); | |
5440 | 528 COLAMD_NAME (_set_defaults) (knobs); |
5164 | 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) | |
5760 | 539 spumoni = static_cast<int> (User_knobs (1)); |
5164 | 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 | |
5440 | 547 octave_idx_type n_row, n_col, nnz; |
548 octave_idx_type *ridx, *cidx; | |
5164 | 549 SparseMatrix sm; |
550 SparseComplexMatrix scm; | |
551 | |
5631 | 552 if (args(0).is_sparse_type ()) |
5164 | 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 (); | |
5604 | 559 nnz = scm.nzmax (); |
5164 | 560 ridx = scm.xridx (); |
561 cidx = scm.xcidx (); | |
562 } | |
563 else | |
564 { | |
565 sm = args(0).sparse_matrix_value (); | |
566 n_row = sm.rows (); | |
567 n_col = sm.cols (); | |
5604 | 568 nnz = sm.nzmax (); |
5164 | 569 ridx = sm.xridx (); |
570 cidx = sm.xcidx (); | |
571 } | |
572 } | |
573 else | |
574 { | |
575 if (args(0).is_complex_type ()) | |
576 sm = SparseMatrix (real (args(0).complex_matrix_value ())); | |
577 else | |
578 sm = SparseMatrix (args(0).matrix_value ()); | |
579 | |
580 n_row = sm.rows (); | |
581 n_col = sm.cols (); | |
5604 | 582 nnz = sm.nzmax (); |
5164 | 583 ridx = sm.xridx (); |
584 cidx = sm.xcidx (); | |
585 } | |
586 | |
587 if (n_row != n_col) | |
588 { | |
589 error ("symamd: matrix must be square"); | |
590 return retval; | |
591 } | |
592 | |
593 // Allocate workspace for symamd | |
5440 | 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, knobs, stats, &calloc, &free)) | |
5164 | 597 { |
5440 | 598 SYMAMD_NAME (_report) (stats) ; |
5164 | 599 error ("symamd: internal error!") ; |
600 return retval; | |
601 } | |
602 | |
603 // column elimination tree post-ordering | |
5440 | 604 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); |
5164 | 605 symetree (ridx, cidx, etree, perm, n_col); |
606 | |
607 // Calculate the tree post-ordering | |
5440 | 608 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1); |
7924 | 609 tree_postorder (n_col, etree, post); |
5164 | 610 |
611 // return the permutation vector | |
612 NDArray out_perm (dim_vector (1, n_col)); | |
5440 | 613 for (octave_idx_type i = 0; i < n_col; i++) |
5164 | 614 out_perm(i) = perm [post [i]] + 1; |
615 | |
7924 | 616 retval(0) = out_perm; |
5164 | 617 |
618 // print stats if spumoni > 0 | |
619 if (spumoni > 0) | |
5440 | 620 SYMAMD_NAME (_report) (stats) ; |
5164 | 621 |
622 // Return the stats vector | |
623 if (nargout == 2) | |
624 { | |
625 NDArray out_stats (dim_vector (1, COLAMD_STATS)); | |
5440 | 626 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++) |
5164 | 627 out_stats (i) = stats [i] ; |
628 retval(1) = out_stats; | |
629 | |
630 // fix stats (5) and (6), for 1-based information on | |
631 // jumbled matrix. note that this correction doesn't | |
632 // occur if symamd returns FALSE | |
633 out_stats (COLAMD_INFO1) ++ ; | |
634 out_stats (COLAMD_INFO2) ++ ; | |
635 } | |
636 } | |
637 | |
5451 | 638 #else |
639 | |
640 error ("symamd: not available in this version of Octave"); | |
641 | |
642 #endif | |
643 | |
5164 | 644 return retval; |
645 } | |
646 | |
647 DEFUN_DLD (etree, args, nargout, | |
648 "-*- texinfo -*-\n\ | |
649 @deftypefn {Loadable Function} {@var{p} =} etree (@var{s})\n\ | |
650 @deftypefnx {Loadable Function} {@var{p} =} etree (@var{s}, @var{typ})\n\ | |
651 @deftypefnx {Loadable Function} {[@var{p}, @var{q}] =} etree (@var{s}, @var{typ})\n\ | |
652 \n\ | |
653 Returns the elimination tree for the matrix @var{s}. By default @var{s}\n\ | |
654 is assumed to be symmetric and the symmetric elimination tree is\n\ | |
655 returned. The argument @var{typ} controls whether a symmetric or\n\ | |
656 column elimination tree is returned. Valid values of @var{typ} are\n\ | |
657 'sym' or 'col', for symmetric or column elimination tree respectively\n\ | |
658 \n\ | |
659 Called with a second argument, @dfn{etree} also returns the postorder\n\ | |
660 permutations on the tree.\n\ | |
661 @end deftypefn") | |
662 { | |
663 octave_value_list retval; | |
5297 | 664 |
5164 | 665 int nargin = args.length (); |
666 | |
667 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2) | |
5823 | 668 print_usage (); |
5164 | 669 else |
670 { | |
5440 | 671 octave_idx_type n_row, n_col, nnz; |
672 octave_idx_type *ridx, *cidx; | |
5164 | 673 bool is_sym = true; |
674 SparseMatrix sm; | |
675 SparseComplexMatrix scm; | |
676 | |
5631 | 677 if (args(0).is_sparse_type ()) |
5164 | 678 { |
679 if (args(0).is_complex_type ()) | |
680 { | |
681 scm = args(0).sparse_complex_matrix_value (); | |
682 n_row = scm.rows (); | |
683 n_col = scm.cols (); | |
5604 | 684 nnz = scm.nzmax (); |
5164 | 685 ridx = scm.xridx (); |
686 cidx = scm.xcidx (); | |
687 } | |
688 else | |
689 { | |
690 sm = args(0).sparse_matrix_value (); | |
691 n_row = sm.rows (); | |
692 n_col = sm.cols (); | |
5604 | 693 nnz = sm.nzmax (); |
5164 | 694 ridx = sm.xridx (); |
695 cidx = sm.xcidx (); | |
696 } | |
697 | |
698 } | |
699 else | |
700 { | |
701 error ("etree: must be called with a sparse matrix"); | |
702 return retval; | |
703 } | |
704 | |
705 if (nargin == 2) | |
6484 | 706 { |
707 if (args(1).is_string ()) | |
708 { | |
709 std::string str = args(1).string_value (); | |
710 if (str.find ("C") == 0 || str.find ("c") == 0) | |
711 is_sym = false; | |
712 } | |
713 else | |
714 { | |
715 error ("etree: second argument must be a string"); | |
716 return retval; | |
717 } | |
718 } | |
7924 | 719 |
5164 | 720 // column elimination tree post-ordering (reuse variables) |
5440 | 721 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); |
5164 | 722 |
723 if (is_sym) | |
724 { | |
725 if (n_row != n_col) | |
726 { | |
727 error ("etree: matrix is marked as symmetric, but not square"); | |
728 return retval; | |
729 } | |
7924 | 730 |
7520 | 731 symetree (ridx, cidx, etree, 0, n_col); |
5164 | 732 } |
733 else | |
734 { | |
5440 | 735 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col); |
736 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col); | |
5164 | 737 |
5440 | 738 for (octave_idx_type i = 0; i < n_col; i++) |
5164 | 739 { |
740 colbeg[i] = cidx[i]; | |
741 colend[i] = cidx[i+1]; | |
742 } | |
743 | |
744 coletree (ridx, colbeg, colend, etree, n_row, n_col); | |
745 } | |
746 | |
747 NDArray tree (dim_vector (1, n_col)); | |
5440 | 748 for (octave_idx_type i = 0; i < n_col; i++) |
5164 | 749 // We flag a root with n_col while Matlab does it with zero |
750 // Convert for matlab compatiable output | |
751 if (etree[i] == n_col) | |
7924 | 752 tree(i) = 0; |
5164 | 753 else |
7924 | 754 tree(i) = etree[i] + 1; |
5164 | 755 |
7924 | 756 retval(0) = tree; |
5164 | 757 |
758 if (nargout == 2) | |
759 { | |
760 // Calculate the tree post-ordering | |
5440 | 761 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1); |
7924 | 762 tree_postorder (n_col, etree, post); |
5164 | 763 |
764 NDArray postorder (dim_vector (1, n_col)); | |
5440 | 765 for (octave_idx_type i = 0; i < n_col; i++) |
7924 | 766 postorder(i) = post[i] + 1; |
5164 | 767 |
7924 | 768 retval(1) = postorder; |
5164 | 769 } |
770 } | |
771 | |
772 return retval; | |
773 } | |
774 | |
775 /* | |
776 ;;; Local Variables: *** | |
777 ;;; mode: C++ *** | |
778 ;;; End: *** | |
779 */ |