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