Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/colamd.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents 01f703952eff
children d0b799dafede
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
50 #else 50 #else
51 #define COLAMD_NAME(name) colamd ## name 51 #define COLAMD_NAME(name) colamd ## name
52 #define SYMAMD_NAME(name) symamd ## name 52 #define SYMAMD_NAME(name) symamd ## name
53 #endif 53 #endif
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)
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 ;
108 108
109 static inline octave_idx_type 109 static inline octave_idx_type
110 find (octave_idx_type i, octave_idx_type *pp) 110 find (octave_idx_type i, octave_idx_type *pp)
111 { 111 {
112 register octave_idx_type p, gp; 112 register octave_idx_type p, gp;
113 113
114 p = pp[i]; 114 p = pp[i];
115 gp = pp[p]; 115 gp = pp[p];
116 116
117 while (gp != p) 117 while (gp != p)
118 { 118 {
124 124
125 return p; 125 return p;
126 } 126 }
127 127
128 static octave_idx_type 128 static octave_idx_type
129 etdfs (octave_idx_type v, octave_idx_type *first_kid, 129 etdfs (octave_idx_type v, octave_idx_type *first_kid,
130 octave_idx_type *next_kid, octave_idx_type *post, 130 octave_idx_type *next_kid, octave_idx_type *post,
131 octave_idx_type postnum) 131 octave_idx_type postnum)
132 { 132 {
133 for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w]) 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); 134 postnum = etdfs (w, first_kid, next_kid, post, postnum);
135 135
148 148
149 // Set up structure describing children 149 // Set up structure describing children
150 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1) 150 for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1)
151 /* do nothing */; 151 /* do nothing */;
152 152
153 for (octave_idx_type v = n-1; v >= 0; v--) 153 for (octave_idx_type v = n-1; v >= 0; v--)
154 { 154 {
155 octave_idx_type dad = parent[v]; 155 octave_idx_type dad = parent[v];
156 next_kid[v] = first_kid[dad]; 156 next_kid[v] = first_kid[dad];
157 first_kid[dad] = v; 157 first_kid[dad] = v;
158 } 158 }
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
173 // Compute firstcol[row] = first nonzero column in row 173 // Compute firstcol[row] = first nonzero column in row
174 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc) 174 for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc)
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
188 // centered at its first vertex, which has the same fill. 188 // centered at its first vertex, which has the same fill.
189 for (octave_idx_type col = 0; col < nc; col++) 189 for (octave_idx_type col = 0; col < nc; col++)
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 }
281 281
282 #ifdef HAVE_COLAMD 282 #ifdef HAVE_COLAMD
283 283
284 int nargin = args.length (); 284 int nargin = args.length ();
285 int spumoni = 0; 285 int spumoni = 0;
286 286
287 if (nargout > 2 || nargin < 1 || nargin > 2) 287 if (nargout > 2 || nargin < 1 || nargin > 2)
288 print_usage (); 288 print_usage ();
289 else 289 else
290 { 290 {
291 // Get knobs 291 // Get knobs
292 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS); 292 OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
293 COLAMD_NAME (_set_defaults) (knobs); 293 COLAMD_NAME (_set_defaults) (knobs);
294 294
295 // Check for user-passed knobs 295 // Check for user-passed knobs
296 if (nargin == 2) 296 if (nargin == 2)
297 { 297 {
298 NDArray User_knobs = args(1).array_value (); 298 NDArray User_knobs = args(1).array_value ();
299 int nel_User_knobs = User_knobs.length (); 299 int nel_User_knobs = User_knobs.length ();
300 300
301 if (nel_User_knobs > 0) 301 if (nel_User_knobs > 0)
302 knobs [COLAMD_DENSE_ROW] = User_knobs (0); 302 knobs [COLAMD_DENSE_ROW] = User_knobs (0);
303 if (nel_User_knobs > 1) 303 if (nel_User_knobs > 1)
304 knobs [COLAMD_DENSE_COL] = User_knobs (1) ; 304 knobs [COLAMD_DENSE_COL] = User_knobs (1) ;
305 if (nel_User_knobs > 2) 305 if (nel_User_knobs > 2)
306 spumoni = static_cast<int> (User_knobs (2)); 306 spumoni = static_cast<int> (User_knobs (2));
307 307
308 // print knob settings if spumoni is set 308 // print knob settings if spumoni is set
309 if (spumoni) 309 if (spumoni)
310 { 310 {
311 311
312 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "." 312 octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "."
313 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n"; 313 << COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n";
314 314
315 if (knobs [COLAMD_DENSE_ROW] >= 0) 315 if (knobs [COLAMD_DENSE_ROW] >= 0)
316 octave_stdout << "knobs(1): " << User_knobs (0) 316 octave_stdout << "knobs(1): " << User_knobs (0)
317 << ", rows with > max(16," 317 << ", rows with > max(16,"
318 << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))" 318 << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))"
319 << " entries removed\n"; 319 << " entries removed\n";
320 else 320 else
321 octave_stdout << "knobs(1): " << User_knobs (0) 321 octave_stdout << "knobs(1): " << User_knobs (0)
322 << ", only completely dense rows removed\n"; 322 << ", only completely dense rows removed\n";
323 323
324 if (knobs [COLAMD_DENSE_COL] >= 0) 324 if (knobs [COLAMD_DENSE_COL] >= 0)
325 octave_stdout << "knobs(2): " << User_knobs (1) 325 octave_stdout << "knobs(2): " << User_knobs (1)
326 << ", cols with > max(16," 326 << ", cols with > max(16,"
327 << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))" 327 << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))"
328 << " entries removed\n"; 328 << " entries removed\n";
329 else 329 else
330 octave_stdout << "knobs(2): " << User_knobs (1) 330 octave_stdout << "knobs(2): " << User_knobs (1)
331 << ", only completely dense columns removed\n"; 331 << ", only completely dense columns removed\n";
332 332
333 octave_stdout << "knobs(3): " << User_knobs (2) 333 octave_stdout << "knobs(3): " << User_knobs (2)
334 << ", statistics and knobs printed\n"; 334 << ", statistics and knobs printed\n";
335 335
336 } 336 }
337 } 337 }
338 338
339 octave_idx_type n_row, n_col, nnz; 339 octave_idx_type n_row, n_col, nnz;
340 octave_idx_type *ridx, *cidx; 340 octave_idx_type *ridx, *cidx;
341 SparseComplexMatrix scm; 341 SparseComplexMatrix scm;
342 SparseMatrix sm; 342 SparseMatrix sm;
343 343
429 NDArray out_stats (dim_vector (1, COLAMD_STATS)); 429 NDArray out_stats (dim_vector (1, COLAMD_STATS));
430 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++) 430 for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
431 out_stats (i) = stats [i] ; 431 out_stats (i) = stats [i] ;
432 retval(1) = out_stats; 432 retval(1) = out_stats;
433 433
434 // fix stats (5) and (6), for 1-based information on 434 // fix stats (5) and (6), for 1-based information on
435 // jumbled matrix. note that this correction doesn't 435 // jumbled matrix. note that this correction doesn't
436 // occur if symamd returns FALSE 436 // occur if symamd returns FALSE
437 out_stats (COLAMD_INFO1) ++ ; 437 out_stats (COLAMD_INFO1) ++ ;
438 out_stats (COLAMD_INFO2) ++ ; 438 out_stats (COLAMD_INFO2) ++ ;
439 } 439 }
440 } 440 }
441 441
442 #else 442 #else
443 443
516 516
517 #ifdef HAVE_COLAMD 517 #ifdef HAVE_COLAMD
518 518
519 int nargin = args.length (); 519 int nargin = args.length ();
520 int spumoni = 0; 520 int spumoni = 0;
521 521
522 if (nargout > 2 || nargin < 1 || nargin > 2) 522 if (nargout > 2 || nargin < 1 || nargin > 2)
523 print_usage (); 523 print_usage ();
524 else 524 else
525 { 525 {
526 // Get knobs 526 // Get knobs
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
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.nnz (); 582 nnz = sm.nnz ();
583 ridx = sm.xridx (); 583 ridx = sm.xridx ();
584 cidx = sm.xcidx (); 584 cidx = sm.xcidx ();
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