Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/colamd.cc @ 5440:b73d469ef0c9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-09-04 12:31:45 by dbateman] ChangeLog
author dbateman
date Sun, 04 Sep 2005 12:31:45 +0000
parents 4c8a2e4e0717
children ed08548b9054
comparison
equal deleted inserted replaced
5439:ca7ebff13a16 5440:b73d469ef0c9
38 #include "ov-re-mat.h" 38 #include "ov-re-mat.h"
39 39
40 #include "ov-re-sparse.h" 40 #include "ov-re-sparse.h"
41 #include "ov-cx-sparse.h" 41 #include "ov-cx-sparse.h"
42 42
43 #if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
44
45 // External COLAMD functions in C 43 // External COLAMD functions in C
46 extern "C" { 44 extern "C" {
47 #include "COLAMD/colamd.h" 45 #include "COLAMD/colamd.h"
48 } 46 }
49 47
48 #ifdef IDX_TYPE_LONG
49 #define COLAMD_NAME(name) colamd_l ## name
50 #define SYMAMD_NAME(name) symamd_l ## name
51 #else
52 #define COLAMD_NAME(name) colamd ## name
53 #define SYMAMD_NAME(name) symamd ## name
54 #endif
55
50 // The symmetric column elimination tree code take from the Davis LDL code. 56 // The symmetric column elimination tree code take from the Davis LDL code.
51 // Copyright given elsewhere in this file. 57 // Copyright given elsewhere in this file.
52 static 58 static
53 void symetree (const int *ridx, const int *cidx, int *Parent, int *P, int n) 59 void symetree (const octave_idx_type *ridx, const octave_idx_type *cidx,
54 { 60 octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n)
55 OCTAVE_LOCAL_BUFFER (int, Flag, n); 61 {
56 OCTAVE_LOCAL_BUFFER (int, Pinv, (P ? n : 0)); 62 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n);
63 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
57 if (P) 64 if (P)
58 // If P is present then compute Pinv, the inverse of P 65 // If P is present then compute Pinv, the inverse of P
59 for (int k = 0 ; k < n ; k++) 66 for (octave_idx_type k = 0 ; k < n ; k++)
60 Pinv [P [k]] = k ; 67 Pinv [P [k]] = k ;
61 68
62 for (int k = 0 ; k < n ; k++) 69 for (octave_idx_type k = 0 ; k < n ; k++)
63 { 70 {
64 // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) 71 // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k)
65 Parent [k] = n ; // parent of k is not yet known 72 Parent [k] = n ; // parent of k is not yet known
66 Flag [k] = k ; // mark node k as visited 73 Flag [k] = k ; // mark node k as visited
67 int kk = (P) ? (P [k]) : (k) ; // kth original, or permuted, column 74 octave_idx_type kk = (P) ? (P [k]) : (k) ; // kth original, or permuted, column
68 int p2 = cidx [kk+1] ; 75 octave_idx_type p2 = cidx [kk+1] ;
69 for (int p = cidx [kk] ; p < p2 ; p++) 76 for (octave_idx_type p = cidx [kk] ; p < p2 ; p++)
70 { 77 {
71 // A (i,k) is nonzero (original or permuted A) 78 // A (i,k) is nonzero (original or permuted A)
72 int i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ; 79 octave_idx_type i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ;
73 if (i < k) 80 if (i < k)
74 { 81 {
75 // follow path from i to root of etree, stop at flagged node 82 // follow path from i to root of etree, stop at flagged node
76 for ( ; Flag [i] != k ; i = Parent [i]) 83 for ( ; Flag [i] != k ; i = Parent [i])
77 { 84 {
85 } 92 }
86 } 93 }
87 94
88 // The elimination tree post-ordering code below is taken from SuperLU 95 // The elimination tree post-ordering code below is taken from SuperLU
89 static inline 96 static inline
90 int make_set (int i, int *pp) 97 octave_idx_type make_set (octave_idx_type i, octave_idx_type *pp)
91 { 98 {
92 pp[i] = i; 99 pp[i] = i;
93 return i; 100 return i;
94 } 101 }
95 102
96 static inline 103 static inline
97 int link (int s, int t, int *pp) 104 octave_idx_type link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp)
98 { 105 {
99 pp[s] = t; 106 pp[s] = t;
100 return t; 107 return t;
101 } 108 }
102 109
103 static inline 110 static inline
104 int find (int i, int *pp) 111 octave_idx_type find (octave_idx_type i, octave_idx_type *pp)
105 { 112 {
106 register int p, gp; 113 register octave_idx_type p, gp;
107 114
108 p = pp[i]; 115 p = pp[i];
109 gp = pp[p]; 116 gp = pp[p];
110 while (gp != p) { 117 while (gp != p) {
111 pp[i] = gp; 118 pp[i] = gp;
115 } 122 }
116 return (p); 123 return (p);
117 } 124 }
118 125
119 static 126 static
120 int etdfs (int v, int *first_kid, int *next_kid, int *post, int postnum) 127 octave_idx_type etdfs (octave_idx_type v, octave_idx_type *first_kid,
121 { 128 octave_idx_type *next_kid, octave_idx_type *post,
122 for (int w = first_kid[v]; w != -1; w = next_kid[w]) { 129 octave_idx_type postnum)
130 {
131 for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w]) {
123 postnum = etdfs (w, first_kid, next_kid, post, postnum); 132 postnum = etdfs (w, first_kid, next_kid, post, postnum);
124 } 133 }
125 post[postnum++] = v; 134 post[postnum++] = v;
126 135
127 return postnum; 136 return postnum;
128 } 137 }
129 138
130 static 139 static
131 void TreePostorder(int n, int *parent, int *post) 140 void TreePostorder(octave_idx_type n, octave_idx_type *parent, octave_idx_type *post)
132 { 141 {
133 // Allocate storage for working arrays and results 142 // Allocate storage for working arrays and results
134 OCTAVE_LOCAL_BUFFER (int, first_kid, n+1); 143 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1);
135 OCTAVE_LOCAL_BUFFER (int, next_kid, n+1); 144 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1);
136 145
137 // Set up structure describing children 146 // Set up structure describing children
138 for (int v = 0; v <= n; first_kid[v++] = -1); 147 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1);
139 for (int v = n-1; v >= 0; v--) 148 for (octave_idx_type v = n-1; v >= 0; v--)
140 { 149 {
141 int dad = parent[v]; 150 octave_idx_type dad = parent[v];
142 next_kid[v] = first_kid[dad]; 151 next_kid[v] = first_kid[dad];
143 first_kid[dad] = v; 152 first_kid[dad] = v;
144 } 153 }
145 154
146 // Depth-first search from dummy root vertex #n 155 // Depth-first search from dummy root vertex #n
147 etdfs (n, first_kid, next_kid, post, 0); 156 etdfs (n, first_kid, next_kid, post, 0);
148 } 157 }
149 158
150 static 159 static
151 void coletree (const int *ridx, const int *colbeg, int *colend, 160 void coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg,
152 int *parent, int nr, int nc) 161 octave_idx_type *colend, octave_idx_type *parent,
153 { 162 octave_idx_type nr, octave_idx_type nc)
154 OCTAVE_LOCAL_BUFFER (int, root, nc); 163 {
155 OCTAVE_LOCAL_BUFFER (int, pp, nc); 164 OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc);
156 OCTAVE_LOCAL_BUFFER (int, firstcol, nr); 165 OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc);
166 OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr);
157 167
158 // Compute firstcol[row] = first nonzero column in row 168 // Compute firstcol[row] = first nonzero column in row
159 for (int row = 0; row < nr; firstcol[row++] = nc); 169 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc);
160 for (int col = 0; col < nc; col++) 170 for (octave_idx_type col = 0; col < nc; col++)
161 for (int p = colbeg[col]; p < colend[col]; p++) 171 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
162 { 172 {
163 int row = ridx[p]; 173 octave_idx_type row = ridx[p];
164 if (firstcol[row] > col) 174 if (firstcol[row] > col)
165 firstcol[row] = col; 175 firstcol[row] = col;
166 } 176 }
167 177
168 // Compute etree by Liu's algorithm for symmetric matrices, 178 // Compute etree by Liu's algorithm for symmetric matrices,
169 // except use (firstcol[r],c) in place of an edge (r,c) of A. 179 // except use (firstcol[r],c) in place of an edge (r,c) of A.
170 // Thus each row clique in A'*A is replaced by a star 180 // Thus each row clique in A'*A is replaced by a star
171 // centered at its first vertex, which has the same fill. 181 // centered at its first vertex, which has the same fill.
172 for (int col = 0; col < nc; col++) 182 for (octave_idx_type col = 0; col < nc; col++)
173 { 183 {
174 int cset = make_set (col, pp); 184 octave_idx_type cset = make_set (col, pp);
175 root[cset] = col; 185 root[cset] = col;
176 parent[col] = nc; 186 parent[col] = nc;
177 for (int p = colbeg[col]; p < colend[col]; p++) 187 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
178 { 188 {
179 int row = firstcol[ridx[p]]; 189 octave_idx_type row = firstcol[ridx[p]];
180 if (row >= col) 190 if (row >= col)
181 continue; 191 continue;
182 int rset = find (row, pp); 192 octave_idx_type rset = find (row, pp);
183 int rroot = root[rset]; 193 octave_idx_type rroot = root[rset];
184 if (rroot != col) 194 if (rroot != col)
185 { 195 {
186 parent[rroot] = col; 196 parent[rroot] = col;
187 cset = link (cset, rset, pp); 197 cset = link (cset, rset, pp);
188 root[cset] = col; 198 root[cset] = col;
189 } 199 }
190 } 200 }
191 } 201 }
192 } 202 }
193
194 #endif
195 203
196 DEFUN_DLD (colamd, args, nargout, 204 DEFUN_DLD (colamd, args, nargout,
197 "-*- texinfo -*-\n\ 205 "-*- texinfo -*-\n\
198 @deftypefn {Loadable Function} {@var{p} =} colamd (@var{s})\n\ 206 @deftypefn {Loadable Function} {@var{p} =} colamd (@var{s})\n\
199 @deftypefnx {Loadable Function} {@var{p} =} colamd (@var{s}, @var{knobs})\n\ 207 @deftypefnx {Loadable Function} {@var{p} =} colamd (@var{s}, @var{knobs})\n\
206 @code{@var{s} (:,@var{p})} tends to have sparser LU factors than @var{s}.\n\ 214 @code{@var{s} (:,@var{p})} tends to have sparser LU factors than @var{s}.\n\
207 The Cholesky factorization of @code{@var{s} (:,@var{p})' * @var{s}\n\ 215 The Cholesky factorization of @code{@var{s} (:,@var{p})' * @var{s}\n\
208 (:,@var{p})} also tends to be sparser than that of @code{@var{s}' *\n\ 216 (:,@var{p})} also tends to be sparser than that of @code{@var{s}' *\n\
209 @var{s}}.\n\ 217 @var{s}}.\n\
210 \n\ 218 \n\
211 @var{knobs} is an optional two-element input vector. If @var{s} is\n\ 219 @var{knobs} is an optional one- to three-element input vector. If @var{s} is\n\
212 m-by-n, then rows with more than @code{(@var{knobs} (1)) * @var{n}}\n\ 220 m-by-n, then rows with more than @code{max(16,@var{knobs}(1)*sqrt(n))} entries\n\
213 entries are ignored. Columns with more than @code{(@var{knobs} (2)) *\n\ 221 are ignored. Columns with more than @code{max(16,knobs(2)*sqrt(min(m,n)))}\n\
214 @var{m}} entries are removed prior to ordering, and ordered last in the\n\ 222 entries are removed prior to ordering, and ordered last in the output\n\
215 output permutation @var{p}. If the knobs parameter is not present, then\n\ 223 permutation @var{p}. Only completely dense rows or columns are removed\n\
216 0.5 is used instead, for both @code{@var{knobs} (1)} and\n\ 224 if @code{@var{knobs} (1)} and @code{@var{knobs} (2)} are < 0, respectively.\n\
217 @code{@var{knobs} (2)}. @code{@var{knobs} (3)} controls the printing of\n\ 225 If @code{@var{knobs} (3)} is nonzero, @var{stats} and @var{knobs} are\n\
218 statistics and error messages.\n\ 226 printed. The default is @code{@var{knobs} = [10 10 0]}. Note that\n\
227 @var{knobs} differs from earlier versions of colamd\n\
219 \n\ 228 \n\
220 @var{stats} is an optional 20-element output vector that provides data\n\ 229 @var{stats} is an optional 20-element output vector that provides data\n\
221 about the ordering and the validity of the input matrix @var{s}. Ordering\n\ 230 about the ordering and the validity of the input matrix @var{s}. Ordering\n\
222 statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1)} and\n\ 231 statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1)} and\n\
223 @code{@var{stats} (2)} are the number of dense or empty rows and columns\n\ 232 @code{@var{stats} (2)} are the number of dense or empty rows and columns\n\
259 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\ 268 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\
260 @end deftypefn\n\ 269 @end deftypefn\n\
261 @seealso{colperm, symamd}") 270 @seealso{colperm, symamd}")
262 { 271 {
263 octave_value_list retval; 272 octave_value_list retval;
264
265 #if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
266
267 int nargin = args.length (); 273 int nargin = args.length ();
268 int spumoni = 0; 274 int spumoni = 0;
269 275
270 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2) 276 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2)
271 usage ("colamd: incorrect number of input and/or output arguments"); 277 usage ("colamd: incorrect number of input and/or output arguments");
272 else 278 else
273 { 279 {
274 // Get knobs 280 // Get knobs
275 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); 281 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
276 colamd_set_defaults (knobs); 282 COLAMD_NAME (_set_defaults) (knobs);
277 283
278 // Check for user-passed knobs 284 // Check for user-passed knobs
279 if (nargin == 2) 285 if (nargin == 2)
280 { 286 {
281 NDArray User_knobs = args(1).array_value (); 287 NDArray User_knobs = args(1).array_value ();
282 int nel_User_knobs = User_knobs.length (); 288 int nel_User_knobs = User_knobs.length ();
283 289
284 if (nel_User_knobs > 0) 290 if (nel_User_knobs > 0)
285 knobs [COLAMD_DENSE_ROW] = User_knobs (COLAMD_DENSE_ROW); 291 knobs [COLAMD_DENSE_ROW] = User_knobs (0);
286 if (nel_User_knobs > 1) 292 if (nel_User_knobs > 1)
287 knobs [COLAMD_DENSE_COL] = User_knobs (COLAMD_DENSE_COL) ; 293 knobs [COLAMD_DENSE_COL] = User_knobs (1) ;
288 if (nel_User_knobs > 2) 294 if (nel_User_knobs > 2)
289 spumoni = (int) User_knobs (2); 295 spumoni = (int) User_knobs (2);
290 } 296
291 297 // print knob settings if spumoni is set
292 // print knob settings if spumoni is set 298 if (spumoni)
293 if (spumoni > 0) 299 {
294 { 300
295 octave_stdout << "colamd: dense row fraction: " 301 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "."
296 << knobs [COLAMD_DENSE_ROW] << std::endl; 302 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n";
297 octave_stdout << "colamd: dense col fraction: " 303
298 << knobs [COLAMD_DENSE_COL] << std::endl; 304 if (knobs [COLAMD_DENSE_ROW] >= 0)
305 octave_stdout << "knobs(1): " << User_knobs (0)
306 << ", rows with > max(16,"
307 << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))"
308 << " entries removed\n";
309 else
310 octave_stdout << "knobs(1): " << User_knobs (0)
311 << ", only completely dense rows removed\n";
312
313 if (knobs [COLAMD_DENSE_COL] >= 0)
314 octave_stdout << "knobs(2): " << User_knobs (1)
315 << ", cols with > max(16,"
316 << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))"
317 << " entries removed\n";
318 else
319 octave_stdout << "knobs(2): " << User_knobs (1)
320 << ", only completely dense columns removed\n";
321
322 octave_stdout << "knobs(3): " << User_knobs (2)
323 << ", statistics and knobs printed\n";
324
325 }
299 } 326 }
300 327
301 int n_row, n_col, nnz; 328 octave_idx_type n_row, n_col, nnz;
302 int *ridx, *cidx; 329 octave_idx_type *ridx, *cidx;
303 SparseComplexMatrix scm; 330 SparseComplexMatrix scm;
304 SparseMatrix sm; 331 SparseMatrix sm;
305 332
306 if (args(0).class_name () == "sparse") 333 if (args(0).class_name () == "sparse")
307 { 334 {
338 ridx = sm.xridx (); 365 ridx = sm.xridx ();
339 cidx = sm.xcidx (); 366 cidx = sm.xcidx ();
340 } 367 }
341 368
342 // Allocate workspace for colamd 369 // Allocate workspace for colamd
343 OCTAVE_LOCAL_BUFFER (int, p, n_col+1); 370 OCTAVE_LOCAL_BUFFER (octave_idx_type, p, n_col+1);
344 for (int i = 0; i < n_col+1; i++) 371 for (octave_idx_type i = 0; i < n_col+1; i++)
345 p[i] = cidx [i]; 372 p[i] = cidx [i];
346 373
347 int Alen = colamd_recommended (nnz, n_row, n_col); 374 octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col);
348 OCTAVE_LOCAL_BUFFER (int, A, Alen); 375 OCTAVE_LOCAL_BUFFER (octave_idx_type, A, Alen);
349 for (int i = 0; i < nnz; i++) 376 for (octave_idx_type i = 0; i < nnz; i++)
350 A[i] = ridx [i]; 377 A[i] = ridx [i];
351 378
352 // Order the columns (destroys A) 379 // Order the columns (destroys A)
353 OCTAVE_LOCAL_BUFFER (int, stats, COLAMD_STATS); 380 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
354 if (!colamd (n_row, n_col, Alen, A, p, knobs, stats)) 381 if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats))
355 { 382 {
356 colamd_report (stats) ; 383 COLAMD_NAME (_report) (stats) ;
357 error ("colamd: internal error!"); 384 error ("colamd: internal error!");
358 return retval; 385 return retval;
359 } 386 }
360 387
361 // column elimination tree post-ordering (reuse variables) 388 // column elimination tree post-ordering (reuse variables)
362 OCTAVE_LOCAL_BUFFER (int, colbeg, n_col + 1); 389 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1);
363 OCTAVE_LOCAL_BUFFER (int, colend, n_col + 1); 390 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col + 1);
364 OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1); 391 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
365 392
366 for (int i = 0; i < n_col; i++) 393 for (octave_idx_type i = 0; i < n_col; i++)
367 { 394 {
368 colbeg[i] = cidx[p[i]]; 395 colbeg[i] = cidx[p[i]];
369 colend[i] = cidx[p[i]+1]; 396 colend[i] = cidx[p[i]+1];
370 } 397 }
371 398
374 // Calculate the tree post-ordering 401 // Calculate the tree post-ordering
375 TreePostorder (n_col, etree, colbeg); 402 TreePostorder (n_col, etree, colbeg);
376 403
377 // return the permutation vector 404 // return the permutation vector
378 NDArray out_perm (dim_vector (1, n_col)); 405 NDArray out_perm (dim_vector (1, n_col));
379 for (int i = 0; i < n_col; i++) 406 for (octave_idx_type i = 0; i < n_col; i++)
380 out_perm(i) = p [colbeg [i]] + 1; 407 out_perm(i) = p [colbeg [i]] + 1;
381 408
382 retval (0) = out_perm; 409 retval (0) = out_perm;
383 410
384 // print stats if spumoni > 0 411 // print stats if spumoni > 0
385 if (spumoni > 0) 412 if (spumoni > 0)
386 colamd_report (stats) ; 413 COLAMD_NAME (_report) (stats) ;
387 414
388 // Return the stats vector 415 // Return the stats vector
389 if (nargout == 2) 416 if (nargout == 2)
390 { 417 {
391 NDArray out_stats (dim_vector (1, COLAMD_STATS)); 418 NDArray out_stats (dim_vector (1, COLAMD_STATS));
392 for (int i = 0 ; i < COLAMD_STATS ; i++) 419 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
393 out_stats (i) = stats [i] ; 420 out_stats (i) = stats [i] ;
394 retval(1) = out_stats; 421 retval(1) = out_stats;
395 422
396 // fix stats (5) and (6), for 1-based information on 423 // fix stats (5) and (6), for 1-based information on
397 // jumbled matrix. note that this correction doesn't 424 // jumbled matrix. note that this correction doesn't
398 // occur if symamd returns FALSE 425 // occur if symamd returns FALSE
399 out_stats (COLAMD_INFO1) ++ ; 426 out_stats (COLAMD_INFO1) ++ ;
400 out_stats (COLAMD_INFO2) ++ ; 427 out_stats (COLAMD_INFO2) ++ ;
401 } 428 }
402 } 429 }
403
404 #else
405
406 error ("colamd: not available in this version of Octave");
407
408 #endif
409 430
410 return retval; 431 return retval;
411 } 432 }
412 433
413 DEFUN_DLD (symamd, args, nargout, 434 DEFUN_DLD (symamd, args, nargout,
422 sparser Cholesky factor than @var{s}. Sometimes SYMAMD works well for\n\ 443 sparser Cholesky factor than @var{s}. Sometimes SYMAMD works well for\n\
423 symmetric indefinite matrices too. The matrix @var{s} is assumed to be\n\ 444 symmetric indefinite matrices too. The matrix @var{s} is assumed to be\n\
424 symmetric; only the strictly lower triangular part is referenced. @var{s}\n\ 445 symmetric; only the strictly lower triangular part is referenced. @var{s}\n\
425 must be square.\n\ 446 must be square.\n\
426 \n\ 447 \n\
427 @var{knobs} is an optional input argument. If @var{s} is n-by-n, then\n\ 448 @var{knobs} is an optional one- to two-element input vector. If @var{s} is\n\
428 rows and columns with more than @code{@var{knobs} (1) * @var{n}} entries\n\ 449 n-by-n, then rows and columns with more than\n\
429 are removed prior to ordering, and ordered last in the output permutation\n\ 450 @code{max(16,@var{knobs}(1)*sqrt(n))} entries are removed prior to ordering,\n\
430 @var{p}. If the @var{knobs} parameter is not present, then the default of\n\ 451 and ordered last in the output permutation @var{p}. No rows/columns are\n\
431 0.5 is used instead. @code{@var{knobs} (2)} controls the printing of\n\ 452 removed if @code{@var{knobs}(1) < 0}. If @code{@var{knobs} (2)} is nonzero,\n\
432 statistics and error messages.\n\ 453 @code{stats} and @var{knobs} are printed. The default is @code{@var{knobs} \n\
454 = [10 0]}. Note that @var{knobs} differs from earlier versions of symamd.\n\
433 \n\ 455 \n\
434 @var{stats} is an optional 20-element output vector that provides data\n\ 456 @var{stats} is an optional 20-element output vector that provides data\n\
435 about the ordering and the validity of the input matrix @var{s}. Ordering\n\ 457 about the ordering and the validity of the input matrix @var{s}. Ordering\n\
436 statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1) =\n\ 458 statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1) =\n\
437 @var{stats} (2)} is the number of dense or empty rows and columns\n\ 459 @var{stats} (2)} is the number of dense or empty rows and columns\n\
473 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\ 495 @url{http://www.cise.ufl.edu/research/sparse/colamd})\n\
474 @end deftypefn\n\ 496 @end deftypefn\n\
475 @seealso{colperm, colamd}") 497 @seealso{colperm, colamd}")
476 { 498 {
477 octave_value_list retval; 499 octave_value_list retval;
478
479 #if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
480
481 int nargin = args.length (); 500 int nargin = args.length ();
482 int spumoni = 0; 501 int spumoni = 0;
483 502
484 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2) 503 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2)
485 usage ("symamd: incorrect number of input and/or output arguments"); 504 usage ("symamd: incorrect number of input and/or output arguments");
486 else 505 else
487 { 506 {
488 // Get knobs 507 // Get knobs
489 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); 508 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
490 colamd_set_defaults (knobs); 509 COLAMD_NAME (_set_defaults) (knobs);
491 510
492 // Check for user-passed knobs 511 // Check for user-passed knobs
493 if (nargin == 2) 512 if (nargin == 2)
494 { 513 {
495 NDArray User_knobs = args(1).array_value (); 514 NDArray User_knobs = args(1).array_value ();
504 // print knob settings if spumoni is set 523 // print knob settings if spumoni is set
505 if (spumoni > 0) 524 if (spumoni > 0)
506 octave_stdout << "symamd: dense row/col fraction: " 525 octave_stdout << "symamd: dense row/col fraction: "
507 << knobs [COLAMD_DENSE_ROW] << std::endl; 526 << knobs [COLAMD_DENSE_ROW] << std::endl;
508 527
509 int n_row, n_col, nnz; 528 octave_idx_type n_row, n_col, nnz;
510 int *ridx, *cidx; 529 octave_idx_type *ridx, *cidx;
511 SparseMatrix sm; 530 SparseMatrix sm;
512 SparseComplexMatrix scm; 531 SparseComplexMatrix scm;
513 532
514 if (args(0).class_name () == "sparse") 533 if (args(0).class_name () == "sparse")
515 { 534 {
551 error ("symamd: matrix must be square"); 570 error ("symamd: matrix must be square");
552 return retval; 571 return retval;
553 } 572 }
554 573
555 // Allocate workspace for symamd 574 // Allocate workspace for symamd
556 OCTAVE_LOCAL_BUFFER (int, perm, n_col+1); 575 OCTAVE_LOCAL_BUFFER (octave_idx_type, perm, n_col+1);
557 OCTAVE_LOCAL_BUFFER (int, stats, COLAMD_STATS); 576 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
558 if (!symamd (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free)) 577 if (!SYMAMD_NAME () (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free))
559 { 578 {
560 symamd_report (stats) ; 579 SYMAMD_NAME (_report) (stats) ;
561 error ("symamd: internal error!") ; 580 error ("symamd: internal error!") ;
562 return retval; 581 return retval;
563 } 582 }
564 583
565 // column elimination tree post-ordering 584 // column elimination tree post-ordering
566 OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1); 585 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
567 symetree (ridx, cidx, etree, perm, n_col); 586 symetree (ridx, cidx, etree, perm, n_col);
568 587
569 // Calculate the tree post-ordering 588 // Calculate the tree post-ordering
570 OCTAVE_LOCAL_BUFFER (int, post, n_col + 1); 589 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
571 TreePostorder (n_col, etree, post); 590 TreePostorder (n_col, etree, post);
572 591
573 // return the permutation vector 592 // return the permutation vector
574 NDArray out_perm (dim_vector (1, n_col)); 593 NDArray out_perm (dim_vector (1, n_col));
575 for (int i = 0; i < n_col; i++) 594 for (octave_idx_type i = 0; i < n_col; i++)
576 out_perm(i) = perm [post [i]] + 1; 595 out_perm(i) = perm [post [i]] + 1;
577 596
578 retval (0) = out_perm; 597 retval (0) = out_perm;
579 598
580 // print stats if spumoni > 0 599 // print stats if spumoni > 0
581 if (spumoni > 0) 600 if (spumoni > 0)
582 symamd_report (stats) ; 601 SYMAMD_NAME (_report) (stats) ;
583 602
584 // Return the stats vector 603 // Return the stats vector
585 if (nargout == 2) 604 if (nargout == 2)
586 { 605 {
587 NDArray out_stats (dim_vector (1, COLAMD_STATS)); 606 NDArray out_stats (dim_vector (1, COLAMD_STATS));
588 for (int i = 0 ; i < COLAMD_STATS ; i++) 607 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
589 out_stats (i) = stats [i] ; 608 out_stats (i) = stats [i] ;
590 retval(1) = out_stats; 609 retval(1) = out_stats;
591 610
592 // fix stats (5) and (6), for 1-based information on 611 // fix stats (5) and (6), for 1-based information on
593 // jumbled matrix. note that this correction doesn't 612 // jumbled matrix. note that this correction doesn't
595 out_stats (COLAMD_INFO1) ++ ; 614 out_stats (COLAMD_INFO1) ++ ;
596 out_stats (COLAMD_INFO2) ++ ; 615 out_stats (COLAMD_INFO2) ++ ;
597 } 616 }
598 } 617 }
599 618
600 #else
601
602 error ("symamd: not available in this version of Octave");
603
604 #endif
605
606 return retval; 619 return retval;
607 } 620 }
608 621
609 DEFUN_DLD (etree, args, nargout, 622 DEFUN_DLD (etree, args, nargout,
610 "-*- texinfo -*-\n\ 623 "-*- texinfo -*-\n\
622 permutations on the tree.\n\ 635 permutations on the tree.\n\
623 @end deftypefn") 636 @end deftypefn")
624 { 637 {
625 octave_value_list retval; 638 octave_value_list retval;
626 639
627 #if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
628
629 int nargin = args.length (); 640 int nargin = args.length ();
630 641
631 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2) 642 if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2)
632 usage ("etree: incorrect number of input and/or output arguments"); 643 usage ("etree: incorrect number of input and/or output arguments");
633 else 644 else
634 { 645 {
635 int n_row, n_col, nnz; 646 octave_idx_type n_row, n_col, nnz;
636 int *ridx, *cidx; 647 octave_idx_type *ridx, *cidx;
637 bool is_sym = true; 648 bool is_sym = true;
638 SparseMatrix sm; 649 SparseMatrix sm;
639 SparseComplexMatrix scm; 650 SparseComplexMatrix scm;
640 651
641 if (args(0).class_name () == "sparse") 652 if (args(0).class_name () == "sparse")
678 error ("etree: second argument must be a string"); 689 error ("etree: second argument must be a string");
679 return retval; 690 return retval;
680 } 691 }
681 692
682 // column elimination tree post-ordering (reuse variables) 693 // column elimination tree post-ordering (reuse variables)
683 OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1); 694 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
684 695
685 696
686 if (is_sym) 697 if (is_sym)
687 { 698 {
688 if (n_row != n_col) 699 if (n_row != n_col)
692 } 703 }
693 symetree (ridx, cidx, etree, NULL, n_col); 704 symetree (ridx, cidx, etree, NULL, n_col);
694 } 705 }
695 else 706 else
696 { 707 {
697 OCTAVE_LOCAL_BUFFER (int, colbeg, n_col); 708 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
698 OCTAVE_LOCAL_BUFFER (int, colend, n_col); 709 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col);
699 710
700 for (int i = 0; i < n_col; i++) 711 for (octave_idx_type i = 0; i < n_col; i++)
701 { 712 {
702 colbeg[i] = cidx[i]; 713 colbeg[i] = cidx[i];
703 colend[i] = cidx[i+1]; 714 colend[i] = cidx[i+1];
704 } 715 }
705 716
706 coletree (ridx, colbeg, colend, etree, n_row, n_col); 717 coletree (ridx, colbeg, colend, etree, n_row, n_col);
707 } 718 }
708 719
709 NDArray tree (dim_vector (1, n_col)); 720 NDArray tree (dim_vector (1, n_col));
710 for (int i = 0; i < n_col; i++) 721 for (octave_idx_type i = 0; i < n_col; i++)
711 // We flag a root with n_col while Matlab does it with zero 722 // We flag a root with n_col while Matlab does it with zero
712 // Convert for matlab compatiable output 723 // Convert for matlab compatiable output
713 if (etree[i] == n_col) 724 if (etree[i] == n_col)
714 tree (i) = 0; 725 tree (i) = 0;
715 else 726 else
718 retval (0) = tree; 729 retval (0) = tree;
719 730
720 if (nargout == 2) 731 if (nargout == 2)
721 { 732 {
722 // Calculate the tree post-ordering 733 // Calculate the tree post-ordering
723 OCTAVE_LOCAL_BUFFER (int, post, n_col + 1); 734 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
724 TreePostorder (n_col, etree, post); 735 TreePostorder (n_col, etree, post);
725 736
726 NDArray postorder (dim_vector (1, n_col)); 737 NDArray postorder (dim_vector (1, n_col));
727 for (int i = 0; i < n_col; i++) 738 for (octave_idx_type i = 0; i < n_col; i++)
728 postorder (i) = post[i] + 1; 739 postorder (i) = post[i] + 1;
729 740
730 retval (1) = postorder; 741 retval (1) = postorder;
731 } 742 }
732 } 743 }
733
734 #else
735
736 error ("etree: not available in this version of Octave");
737
738 #endif
739 744
740 return retval; 745 return retval;
741 } 746 }
742 747
743 DEFUN_DLD (symbfact, args, nargout, 748 DEFUN_DLD (symbfact, args, nargout,