Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/colamd.cc @ 8377:25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 29 Oct 2008 16:52:10 +0100 |
| parents | 4976f66d469b |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 5164 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 2004, 2005, 2006, 2007 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 */ |