Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/colamd.cc @ 10154:40dfc0c99116

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DLD-FUNCTIONS/*.cc: untabify
author John W. Eaton <jwe@octave.org>
date Wed, 20 Jan 2010 17:33:41 -0500
parents 7c02ec148a3c
children d0ce5e973937
comparison
equal deleted inserted replaced
10153:2c28f9d0360f 10154:40dfc0c99116
54 54
55 // The symmetric column elimination tree code take from the Davis LDL code. 55 // The symmetric column elimination tree code take from the Davis LDL code.
56 // Copyright given elsewhere in this file. 56 // Copyright given elsewhere in this file.
57 static void 57 static void
58 symetree (const octave_idx_type *ridx, const octave_idx_type *cidx, 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) 59 octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n)
60 { 60 {
61 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n); 61 OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n);
62 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0)); 62 OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
63 if (P) 63 if (P)
64 // If P is present then compute Pinv, the inverse of P 64 // If P is present then compute Pinv, the inverse of P
66 Pinv [P [k]] = k ; 66 Pinv [P [k]] = k ;
67 67
68 for (octave_idx_type k = 0 ; k < n ; k++) 68 for (octave_idx_type k = 0 ; k < n ; k++)
69 { 69 {
70 // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) 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 71 Parent [k] = n ; // parent of k is not yet known
72 Flag [k] = k ; // mark node k as visited 72 Flag [k] = k ; // mark node k as visited
73 octave_idx_type kk = (P) ? (P [k]) : (k) ; // kth original, or permuted, column 73 octave_idx_type kk = (P) ? (P [k]) : (k) ; // kth original, or permuted, column
74 octave_idx_type p2 = cidx [kk+1] ; 74 octave_idx_type p2 = cidx [kk+1] ;
75 for (octave_idx_type p = cidx [kk] ; p < p2 ; p++) 75 for (octave_idx_type p = cidx [kk] ; p < p2 ; p++)
76 { 76 {
77 // A (i,k) is nonzero (original or permuted A) 77 // A (i,k) is nonzero (original or permuted A)
78 octave_idx_type i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ; 78 octave_idx_type i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ;
79 if (i < k) 79 if (i < k)
80 { 80 {
81 // follow path from i to root of etree, stop at flagged node 81 // follow path from i to root of etree, stop at flagged node
82 for ( ; Flag [i] != k ; i = Parent [i]) 82 for ( ; Flag [i] != k ; i = Parent [i])
83 { 83 {
84 // find parent of i if not yet determined 84 // find parent of i if not yet determined
85 if (Parent [i] == n) 85 if (Parent [i] == n)
86 Parent [i] = k ; 86 Parent [i] = k ;
87 Flag [i] = k ; // mark i as visited 87 Flag [i] = k ; // mark i as visited
88 } 88 }
89 } 89 }
90 } 90 }
91 } 91 }
92 } 92 }
93 93
94 // The elimination tree post-ordering code below is taken from SuperLU 94 // The elimination tree post-ordering code below is taken from SuperLU
95 static inline octave_idx_type 95 static inline octave_idx_type
138 return postnum; 138 return postnum;
139 } 139 }
140 140
141 static void 141 static void
142 tree_postorder (octave_idx_type n, octave_idx_type *parent, 142 tree_postorder (octave_idx_type n, octave_idx_type *parent,
143 octave_idx_type *post) 143 octave_idx_type *post)
144 { 144 {
145 // Allocate storage for working arrays and results 145 // Allocate storage for working arrays and results
146 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1); 146 OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1);
147 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1); 147 OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1);
148 148
161 etdfs (n, first_kid, next_kid, post, 0); 161 etdfs (n, first_kid, next_kid, post, 0);
162 } 162 }
163 163
164 static void 164 static void
165 coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg, 165 coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg,
166 octave_idx_type *colend, octave_idx_type *parent, 166 octave_idx_type *colend, octave_idx_type *parent,
167 octave_idx_type nr, octave_idx_type nc) 167 octave_idx_type nr, octave_idx_type nc)
168 { 168 {
169 OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc); 169 OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc);
170 OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc); 170 OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc);
171 OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr); 171 OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr);
172 172
175 /* do nothing */; 175 /* do nothing */;
176 176
177 for (octave_idx_type col = 0; col < nc; col++) 177 for (octave_idx_type col = 0; col < nc; col++)
178 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) 178 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
179 { 179 {
180 octave_idx_type row = ridx[p]; 180 octave_idx_type row = ridx[p];
181 if (firstcol[row] > col) 181 if (firstcol[row] > col)
182 firstcol[row] = col; 182 firstcol[row] = col;
183 } 183 }
184 184
185 // Compute etree by Liu's algorithm for symmetric matrices, 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. 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 187 // Thus each row clique in A'*A is replaced by a star
190 { 190 {
191 octave_idx_type cset = make_set (col, pp); 191 octave_idx_type cset = make_set (col, pp);
192 root[cset] = col; 192 root[cset] = col;
193 parent[col] = nc; 193 parent[col] = nc;
194 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) 194 for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
195 { 195 {
196 octave_idx_type row = firstcol[ridx[p]]; 196 octave_idx_type row = firstcol[ridx[p]];
197 if (row >= col) 197 if (row >= col)
198 continue; 198 continue;
199 octave_idx_type rset = find (row, pp); 199 octave_idx_type rset = find (row, pp);
200 octave_idx_type rroot = root[rset]; 200 octave_idx_type rroot = root[rset];
201 if (rroot != col) 201 if (rroot != col)
202 { 202 {
203 parent[rroot] = col; 203 parent[rroot] = col;
204 cset = link (cset, rset, pp); 204 cset = link (cset, rset, pp);
205 root[cset] = col; 205 root[cset] = col;
206 } 206 }
207 } 207 }
208 } 208 }
209 } 209 }
210 210
211 DEFUN_DLD (colamd, args, nargout, 211 DEFUN_DLD (colamd, args, nargout,
212 "-*- texinfo -*-\n\ 212 "-*- texinfo -*-\n\
291 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); 291 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
292 COLAMD_NAME (_set_defaults) (knobs); 292 COLAMD_NAME (_set_defaults) (knobs);
293 293
294 // Check for user-passed knobs 294 // Check for user-passed knobs
295 if (nargin == 2) 295 if (nargin == 2)
296 { 296 {
297 NDArray User_knobs = args(1).array_value (); 297 NDArray User_knobs = args(1).array_value ();
298 int nel_User_knobs = User_knobs.length (); 298 int nel_User_knobs = User_knobs.length ();
299 299
300 if (nel_User_knobs > 0) 300 if (nel_User_knobs > 0)
301 knobs [COLAMD_DENSE_ROW] = User_knobs (0); 301 knobs [COLAMD_DENSE_ROW] = User_knobs (0);
302 if (nel_User_knobs > 1) 302 if (nel_User_knobs > 1)
303 knobs [COLAMD_DENSE_COL] = User_knobs (1) ; 303 knobs [COLAMD_DENSE_COL] = User_knobs (1) ;
304 if (nel_User_knobs > 2) 304 if (nel_User_knobs > 2)
305 spumoni = static_cast<int> (User_knobs (2)); 305 spumoni = static_cast<int> (User_knobs (2));
306 306
307 // print knob settings if spumoni is set 307 // print knob settings if spumoni is set
308 if (spumoni) 308 if (spumoni)
309 { 309 {
310 310
311 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "." 311 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "."
312 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n"; 312 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n";
313 313
314 if (knobs [COLAMD_DENSE_ROW] >= 0) 314 if (knobs [COLAMD_DENSE_ROW] >= 0)
315 octave_stdout << "knobs(1): " << User_knobs (0) 315 octave_stdout << "knobs(1): " << User_knobs (0)
316 << ", rows with > max(16," 316 << ", rows with > max(16,"
317 << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))" 317 << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))"
318 << " entries removed\n"; 318 << " entries removed\n";
319 else 319 else
320 octave_stdout << "knobs(1): " << User_knobs (0) 320 octave_stdout << "knobs(1): " << User_knobs (0)
321 << ", only completely dense rows removed\n"; 321 << ", only completely dense rows removed\n";
322 322
323 if (knobs [COLAMD_DENSE_COL] >= 0) 323 if (knobs [COLAMD_DENSE_COL] >= 0)
324 octave_stdout << "knobs(2): " << User_knobs (1) 324 octave_stdout << "knobs(2): " << User_knobs (1)
325 << ", cols with > max(16," 325 << ", cols with > max(16,"
326 << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))" 326 << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))"
327 << " entries removed\n"; 327 << " entries removed\n";
328 else 328 else
329 octave_stdout << "knobs(2): " << User_knobs (1) 329 octave_stdout << "knobs(2): " << User_knobs (1)
330 << ", only completely dense columns removed\n"; 330 << ", only completely dense columns removed\n";
331 331
332 octave_stdout << "knobs(3): " << User_knobs (2) 332 octave_stdout << "knobs(3): " << User_knobs (2)
333 << ", statistics and knobs printed\n"; 333 << ", statistics and knobs printed\n";
334 334
335 } 335 }
336 } 336 }
337 337
338 octave_idx_type n_row, n_col, nnz; 338 octave_idx_type n_row, n_col, nnz;
339 octave_idx_type *ridx, *cidx; 339 octave_idx_type *ridx, *cidx;
340 SparseComplexMatrix scm; 340 SparseComplexMatrix scm;
341 SparseMatrix sm; 341 SparseMatrix sm;
342 342
343 if (args(0).is_sparse_type ()) 343 if (args(0).is_sparse_type ())
344 { 344 {
345 if (args(0).is_complex_type ()) 345 if (args(0).is_complex_type ())
346 { 346 {
347 scm = args(0). sparse_complex_matrix_value (); 347 scm = args(0). sparse_complex_matrix_value ();
348 n_row = scm.rows (); 348 n_row = scm.rows ();
349 n_col = scm.cols (); 349 n_col = scm.cols ();
350 nnz = scm.nzmax (); 350 nnz = scm.nzmax ();
351 ridx = scm.xridx (); 351 ridx = scm.xridx ();
352 cidx = scm.xcidx (); 352 cidx = scm.xcidx ();
353 } 353 }
354 else 354 else
355 { 355 {
356 sm = args(0).sparse_matrix_value (); 356 sm = args(0).sparse_matrix_value ();
357 357
358 n_row = sm.rows (); 358 n_row = sm.rows ();
359 n_col = sm.cols (); 359 n_col = sm.cols ();
360 nnz = sm.nzmax (); 360 nnz = sm.nzmax ();
361 ridx = sm.xridx (); 361 ridx = sm.xridx ();
362 cidx = sm.xcidx (); 362 cidx = sm.xcidx ();
363 } 363 }
364 } 364 }
365 else 365 else
366 { 366 {
367 if (args(0).is_complex_type ()) 367 if (args(0).is_complex_type ())
368 sm = SparseMatrix (real (args(0).complex_matrix_value ())); 368 sm = SparseMatrix (real (args(0).complex_matrix_value ()));
369 else 369 else
370 sm = SparseMatrix (args(0).matrix_value ()); 370 sm = SparseMatrix (args(0).matrix_value ());
371 371
372 n_row = sm.rows (); 372 n_row = sm.rows ();
373 n_col = sm.cols (); 373 n_col = sm.cols ();
374 nnz = sm.nzmax (); 374 nnz = sm.nzmax ();
375 ridx = sm.xridx (); 375 ridx = sm.xridx ();
376 cidx = sm.xcidx (); 376 cidx = sm.xcidx ();
377 } 377 }
378 378
379 // Allocate workspace for colamd 379 // Allocate workspace for colamd
380 OCTAVE_LOCAL_BUFFER (octave_idx_type, p, n_col+1); 380 OCTAVE_LOCAL_BUFFER (octave_idx_type, p, n_col+1);
381 for (octave_idx_type i = 0; i < n_col+1; i++) 381 for (octave_idx_type i = 0; i < n_col+1; i++)
382 p[i] = cidx [i]; 382 p[i] = cidx [i];
383 383
384 octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col); 384 octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col);
385 OCTAVE_LOCAL_BUFFER (octave_idx_type, A, Alen); 385 OCTAVE_LOCAL_BUFFER (octave_idx_type, A, Alen);
386 for (octave_idx_type i = 0; i < nnz; i++) 386 for (octave_idx_type i = 0; i < nnz; i++)
387 A[i] = ridx [i]; 387 A[i] = ridx [i];
388 388
389 // Order the columns (destroys A) 389 // Order the columns (destroys A)
390 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS); 390 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
391 if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats)) 391 if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats))
392 { 392 {
393 COLAMD_NAME (_report) (stats) ; 393 COLAMD_NAME (_report) (stats) ;
394 error ("colamd: internal error!"); 394 error ("colamd: internal error!");
395 return retval; 395 return retval;
396 } 396 }
397 397
398 // column elimination tree post-ordering (reuse variables) 398 // column elimination tree post-ordering (reuse variables)
399 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1); 399 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1);
400 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, 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); 401 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
402 402
403 for (octave_idx_type i = 0; i < n_col; i++) 403 for (octave_idx_type i = 0; i < n_col; i++)
404 { 404 {
405 colbeg[i] = cidx[p[i]]; 405 colbeg[i] = cidx[p[i]];
406 colend[i] = cidx[p[i]+1]; 406 colend[i] = cidx[p[i]+1];
407 } 407 }
408 408
409 coletree (ridx, colbeg, colend, etree, n_row, n_col); 409 coletree (ridx, colbeg, colend, etree, n_row, n_col);
410 410
411 // Calculate the tree post-ordering 411 // Calculate the tree post-ordering
412 tree_postorder (n_col, etree, colbeg); 412 tree_postorder (n_col, etree, colbeg);
413 413
414 // return the permutation vector 414 // return the permutation vector
415 NDArray out_perm (dim_vector (1, n_col)); 415 NDArray out_perm (dim_vector (1, n_col));
416 for (octave_idx_type i = 0; i < n_col; i++) 416 for (octave_idx_type i = 0; i < n_col; i++)
417 out_perm(i) = p [colbeg [i]] + 1; 417 out_perm(i) = p [colbeg [i]] + 1;
418 418
419 retval(0) = out_perm; 419 retval(0) = out_perm;
420 420
421 // print stats if spumoni > 0 421 // print stats if spumoni > 0
422 if (spumoni > 0) 422 if (spumoni > 0)
423 COLAMD_NAME (_report) (stats) ; 423 COLAMD_NAME (_report) (stats) ;
424 424
425 // Return the stats vector 425 // Return the stats vector
426 if (nargout == 2) 426 if (nargout == 2)
427 { 427 {
428 NDArray out_stats (dim_vector (1, COLAMD_STATS)); 428 NDArray out_stats (dim_vector (1, COLAMD_STATS));
429 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++) 429 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
430 out_stats (i) = stats [i] ; 430 out_stats (i) = stats [i] ;
431 retval(1) = out_stats; 431 retval(1) = out_stats;
432 432
433 // fix stats (5) and (6), for 1-based information on 433 // fix stats (5) and (6), for 1-based information on
434 // jumbled matrix. note that this correction doesn't 434 // jumbled matrix. note that this correction doesn't
435 // occur if symamd returns FALSE 435 // occur if symamd returns FALSE
436 out_stats (COLAMD_INFO1) ++ ; 436 out_stats (COLAMD_INFO1) ++ ;
437 out_stats (COLAMD_INFO2) ++ ; 437 out_stats (COLAMD_INFO2) ++ ;
438 } 438 }
439 } 439 }
440 440
441 #else 441 #else
442 442
443 error ("colamd: not available in this version of Octave"); 443 error ("colamd: not available in this version of Octave");
527 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); 527 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
528 COLAMD_NAME (_set_defaults) (knobs); 528 COLAMD_NAME (_set_defaults) (knobs);
529 529
530 // Check for user-passed knobs 530 // Check for user-passed knobs
531 if (nargin == 2) 531 if (nargin == 2)
532 { 532 {
533 NDArray User_knobs = args(1).array_value (); 533 NDArray User_knobs = args(1).array_value ();
534 int nel_User_knobs = User_knobs.length (); 534 int nel_User_knobs = User_knobs.length ();
535 535
536 if (nel_User_knobs > 0) 536 if (nel_User_knobs > 0)
537 knobs [COLAMD_DENSE_ROW] = User_knobs (COLAMD_DENSE_ROW); 537 knobs [COLAMD_DENSE_ROW] = User_knobs (COLAMD_DENSE_ROW);
538 if (nel_User_knobs > 1) 538 if (nel_User_knobs > 1)
539 spumoni = static_cast<int> (User_knobs (1)); 539 spumoni = static_cast<int> (User_knobs (1));
540 } 540 }
541 541
542 // print knob settings if spumoni is set 542 // print knob settings if spumoni is set
543 if (spumoni > 0) 543 if (spumoni > 0)
544 octave_stdout << "symamd: dense row/col fraction: " 544 octave_stdout << "symamd: dense row/col fraction: "
545 << knobs [COLAMD_DENSE_ROW] << std::endl; 545 << knobs [COLAMD_DENSE_ROW] << std::endl;
546 546
547 octave_idx_type n_row, n_col, nnz; 547 octave_idx_type n_row, n_col, nnz;
548 octave_idx_type *ridx, *cidx; 548 octave_idx_type *ridx, *cidx;
549 SparseMatrix sm; 549 SparseMatrix sm;
550 SparseComplexMatrix scm; 550 SparseComplexMatrix scm;
551 551
552 if (args(0).is_sparse_type ()) 552 if (args(0).is_sparse_type ())
553 { 553 {
554 if (args(0).is_complex_type ()) 554 if (args(0).is_complex_type ())
555 { 555 {
556 scm = args(0).sparse_complex_matrix_value (); 556 scm = args(0).sparse_complex_matrix_value ();
557 n_row = scm.rows (); 557 n_row = scm.rows ();
558 n_col = scm.cols (); 558 n_col = scm.cols ();
559 nnz = scm.nzmax (); 559 nnz = scm.nzmax ();
560 ridx = scm.xridx (); 560 ridx = scm.xridx ();
561 cidx = scm.xcidx (); 561 cidx = scm.xcidx ();
562 } 562 }
563 else 563 else
564 { 564 {
565 sm = args(0).sparse_matrix_value (); 565 sm = args(0).sparse_matrix_value ();
566 n_row = sm.rows (); 566 n_row = sm.rows ();
567 n_col = sm.cols (); 567 n_col = sm.cols ();
568 nnz = sm.nzmax (); 568 nnz = sm.nzmax ();
569 ridx = sm.xridx (); 569 ridx = sm.xridx ();
570 cidx = sm.xcidx (); 570 cidx = sm.xcidx ();
571 } 571 }
572 } 572 }
573 else 573 else
574 { 574 {
575 if (args(0).is_complex_type ()) 575 if (args(0).is_complex_type ())
576 sm = SparseMatrix (real (args(0).complex_matrix_value ())); 576 sm = SparseMatrix (real (args(0).complex_matrix_value ()));
577 else 577 else
578 sm = SparseMatrix (args(0).matrix_value ()); 578 sm = SparseMatrix (args(0).matrix_value ());
579 579
580 n_row = sm.rows (); 580 n_row = sm.rows ();
581 n_col = sm.cols (); 581 n_col = sm.cols ();
582 nnz = sm.nzmax (); 582 nnz = sm.nzmax ();
583 ridx = sm.xridx (); 583 ridx = sm.xridx ();
584 cidx = sm.xcidx (); 584 cidx = sm.xcidx ();
585 } 585 }
586 586
587 if (n_row != n_col) 587 if (n_row != n_col)
588 { 588 {
589 error ("symamd: matrix must be square"); 589 error ("symamd: matrix must be square");
590 return retval; 590 return retval;
591 } 591 }
592 592
593 // Allocate workspace for symamd 593 // Allocate workspace for symamd
594 OCTAVE_LOCAL_BUFFER (octave_idx_type, perm, n_col+1); 594 OCTAVE_LOCAL_BUFFER (octave_idx_type, perm, n_col+1);
595 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS); 595 OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
596 if (!SYMAMD_NAME () (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free)) 596 if (!SYMAMD_NAME () (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free))
597 { 597 {
598 SYMAMD_NAME (_report) (stats) ; 598 SYMAMD_NAME (_report) (stats) ;
599 error ("symamd: internal error!") ; 599 error ("symamd: internal error!") ;
600 return retval; 600 return retval;
601 } 601 }
602 602
603 // column elimination tree post-ordering 603 // column elimination tree post-ordering
604 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); 604 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
605 symetree (ridx, cidx, etree, perm, n_col); 605 symetree (ridx, cidx, etree, perm, n_col);
606 606
609 tree_postorder (n_col, etree, post); 609 tree_postorder (n_col, etree, post);
610 610
611 // return the permutation vector 611 // return the permutation vector
612 NDArray out_perm (dim_vector (1, n_col)); 612 NDArray out_perm (dim_vector (1, n_col));
613 for (octave_idx_type i = 0; i < n_col; i++) 613 for (octave_idx_type i = 0; i < n_col; i++)
614 out_perm(i) = perm [post [i]] + 1; 614 out_perm(i) = perm [post [i]] + 1;
615 615
616 retval(0) = out_perm; 616 retval(0) = out_perm;
617 617
618 // print stats if spumoni > 0 618 // print stats if spumoni > 0
619 if (spumoni > 0) 619 if (spumoni > 0)
620 SYMAMD_NAME (_report) (stats) ; 620 SYMAMD_NAME (_report) (stats) ;
621 621
622 // Return the stats vector 622 // Return the stats vector
623 if (nargout == 2) 623 if (nargout == 2)
624 { 624 {
625 NDArray out_stats (dim_vector (1, COLAMD_STATS)); 625 NDArray out_stats (dim_vector (1, COLAMD_STATS));
626 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++) 626 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
627 out_stats (i) = stats [i] ; 627 out_stats (i) = stats [i] ;
628 retval(1) = out_stats; 628 retval(1) = out_stats;
629 629
630 // fix stats (5) and (6), for 1-based information on 630 // fix stats (5) and (6), for 1-based information on
631 // jumbled matrix. note that this correction doesn't 631 // jumbled matrix. note that this correction doesn't
632 // occur if symamd returns FALSE 632 // occur if symamd returns FALSE
633 out_stats (COLAMD_INFO1) ++ ; 633 out_stats (COLAMD_INFO1) ++ ;
634 out_stats (COLAMD_INFO2) ++ ; 634 out_stats (COLAMD_INFO2) ++ ;
635 } 635 }
636 } 636 }
637 637
638 #else 638 #else
639 639
640 error ("symamd: not available in this version of Octave"); 640 error ("symamd: not available in this version of Octave");
673 bool is_sym = true; 673 bool is_sym = true;
674 SparseMatrix sm; 674 SparseMatrix sm;
675 SparseComplexMatrix scm; 675 SparseComplexMatrix scm;
676 676
677 if (args(0).is_sparse_type ()) 677 if (args(0).is_sparse_type ())
678 { 678 {
679 if (args(0).is_complex_type ()) 679 if (args(0).is_complex_type ())
680 { 680 {
681 scm = args(0).sparse_complex_matrix_value (); 681 scm = args(0).sparse_complex_matrix_value ();
682 n_row = scm.rows (); 682 n_row = scm.rows ();
683 n_col = scm.cols (); 683 n_col = scm.cols ();
684 nnz = scm.nzmax (); 684 nnz = scm.nzmax ();
685 ridx = scm.xridx (); 685 ridx = scm.xridx ();
686 cidx = scm.xcidx (); 686 cidx = scm.xcidx ();
687 } 687 }
688 else 688 else
689 { 689 {
690 sm = args(0).sparse_matrix_value (); 690 sm = args(0).sparse_matrix_value ();
691 n_row = sm.rows (); 691 n_row = sm.rows ();
692 n_col = sm.cols (); 692 n_col = sm.cols ();
693 nnz = sm.nzmax (); 693 nnz = sm.nzmax ();
694 ridx = sm.xridx (); 694 ridx = sm.xridx ();
695 cidx = sm.xcidx (); 695 cidx = sm.xcidx ();
696 } 696 }
697 697
698 } 698 }
699 else 699 else
700 { 700 {
701 error ("etree: must be called with a sparse matrix"); 701 error ("etree: must be called with a sparse matrix");
702 return retval; 702 return retval;
703 } 703 }
704 704
705 if (nargin == 2) 705 if (nargin == 2)
706 { 706 {
707 if (args(1).is_string ()) 707 if (args(1).is_string ())
708 { 708 {
709 std::string str = args(1).string_value (); 709 std::string str = args(1).string_value ();
710 if (str.find ("C") == 0 || str.find ("c") == 0) 710 if (str.find ("C") == 0 || str.find ("c") == 0)
711 is_sym = false; 711 is_sym = false;
712 } 712 }
713 else 713 else
714 { 714 {
715 error ("etree: second argument must be a string"); 715 error ("etree: second argument must be a string");
716 return retval; 716 return retval;
717 } 717 }
718 } 718 }
719 719
720 // column elimination tree post-ordering (reuse variables) 720 // column elimination tree post-ordering (reuse variables)
721 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1); 721 OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
722 722
723 if (is_sym) 723 if (is_sym)
724 { 724 {
725 if (n_row != n_col) 725 if (n_row != n_col)
726 { 726 {
727 error ("etree: matrix is marked as symmetric, but not square"); 727 error ("etree: matrix is marked as symmetric, but not square");
728 return retval; 728 return retval;
729 } 729 }
730 730
731 symetree (ridx, cidx, etree, 0, n_col); 731 symetree (ridx, cidx, etree, 0, n_col);
732 } 732 }
733 else 733 else
734 { 734 {
735 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col); 735 OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
736 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col); 736 OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col);
737 737
738 for (octave_idx_type i = 0; i < n_col; i++) 738 for (octave_idx_type i = 0; i < n_col; i++)
739 { 739 {
740 colbeg[i] = cidx[i]; 740 colbeg[i] = cidx[i];
741 colend[i] = cidx[i+1]; 741 colend[i] = cidx[i+1];
742 } 742 }
743 743
744 coletree (ridx, colbeg, colend, etree, n_row, n_col); 744 coletree (ridx, colbeg, colend, etree, n_row, n_col);
745 } 745 }
746 746
747 NDArray tree (dim_vector (1, n_col)); 747 NDArray tree (dim_vector (1, n_col));
748 for (octave_idx_type i = 0; i < n_col; i++) 748 for (octave_idx_type i = 0; i < n_col; i++)
749 // We flag a root with n_col while Matlab does it with zero 749 // We flag a root with n_col while Matlab does it with zero
750 // Convert for matlab compatiable output 750 // Convert for matlab compatiable output
751 if (etree[i] == n_col) 751 if (etree[i] == n_col)
752 tree(i) = 0; 752 tree(i) = 0;
753 else 753 else
754 tree(i) = etree[i] + 1; 754 tree(i) = etree[i] + 1;
755 755
756 retval(0) = tree; 756 retval(0) = tree;
757 757
758 if (nargout == 2) 758 if (nargout == 2)
759 { 759 {
760 // Calculate the tree post-ordering 760 // Calculate the tree post-ordering
761 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1); 761 OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
762 tree_postorder (n_col, etree, post); 762 tree_postorder (n_col, etree, post);
763 763
764 NDArray postorder (dim_vector (1, n_col)); 764 NDArray postorder (dim_vector (1, n_col));
765 for (octave_idx_type i = 0; i < n_col; i++) 765 for (octave_idx_type i = 0; i < n_col; i++)
766 postorder(i) = post[i] + 1; 766 postorder(i) = post[i] + 1;
767 767
768 retval(1) = postorder; 768 retval(1) = postorder;
769 } 769 }
770 } 770 }
771 771
772 return retval; 772 return retval;
773 } 773 }
774 774