1 // // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10 // This file is modified from the colamd/symamd library. The copyright is below
11
12 // The authors of the code itself are Stefan I. Larimore and Timothy A.
13 // Davis (davis@cise.ufl.edu), University of Florida. The algorithm was
14 // developed in collaboration with John Gilbert, Xerox PARC, and Esmond
15 // Ng, Oak Ridge National Laboratory.
16 //
17 // Date:
18 //
19 // September 8, 2003. Version 2.3.
20 //
21 // Acknowledgements:
22 //
23 // This work was supported by the National Science Foundation, under
24 // grants DMS-9504974 and DMS-9803599.
25 //
26 // Notice:
27 //
28 // Copyright (c) 1998-2003 by the University of Florida.
29 // All Rights Reserved.
30 //
31 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
32 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
33 //
34 // Permission is hereby granted to use, copy, modify, and/or distribute
35 // this program, provided that the Copyright, this License, and the
36 // Availability of the original version is retained on all copies and made
37 // accessible to the end-user of any code or package that includes COLAMD
38 // or any modified version of COLAMD.
39 //
40 // Availability:
41 //
42 // The colamd/symamd library is available at
43 //
44 // http://www.suitesparse.com
45
46
47 #ifndef EIGEN_COLAMD_H
48 #define EIGEN_COLAMD_H
49
50 namespace internal {
51 /* Ensure that debugging is turned off: */
52 #ifndef COLAMD_NDEBUG
53 #define COLAMD_NDEBUG
54 #endif /* NDEBUG */
55 /* ========================================================================== */
56 /* === Knob and statistics definitions ====================================== */
57 /* ========================================================================== */
58
59 /* size of the knobs [ ] array. Only knobs [0..1] are currently used. */
60 #define COLAMD_KNOBS 20
61
62 /* number of output statistics. Only stats [0..6] are currently used. */
63 #define COLAMD_STATS 20
64
65 /* knobs [0] and stats [0]: dense row knob and output statistic. */
66 #define COLAMD_DENSE_ROW 0
67
68 /* knobs [1] and stats [1]: dense column knob and output statistic. */
69 #define COLAMD_DENSE_COL 1
70
71 /* stats [2]: memory defragmentation count output statistic */
72 #define COLAMD_DEFRAG_COUNT 2
73
74 /* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */
75 #define COLAMD_STATUS 3
76
77 /* stats [4..6]: error info, or info on jumbled columns */
78 #define COLAMD_INFO1 4
79 #define COLAMD_INFO2 5
80 #define COLAMD_INFO3 6
81
82 /* error codes returned in stats [3]: */
83 #define COLAMD_OK (0)
84 #define COLAMD_OK_BUT_JUMBLED (1)
85 #define COLAMD_ERROR_A_not_present (-1)
86 #define COLAMD_ERROR_p_not_present (-2)
87 #define COLAMD_ERROR_nrow_negative (-3)
88 #define COLAMD_ERROR_ncol_negative (-4)
89 #define COLAMD_ERROR_nnz_negative (-5)
90 #define COLAMD_ERROR_p0_nonzero (-6)
91 #define COLAMD_ERROR_A_too_small (-7)
92 #define COLAMD_ERROR_col_length_negative (-8)
93 #define COLAMD_ERROR_row_index_out_of_bounds (-9)
94 #define COLAMD_ERROR_out_of_memory (-10)
95 #define COLAMD_ERROR_internal_error (-999)
96
97 /* ========================================================================== */
98 /* === Definitions ========================================================== */
99 /* ========================================================================== */
100
101 #define ONES_COMPLEMENT(r) (-(r)-1)
102
103 /* -------------------------------------------------------------------------- */
104
105 #define COLAMD_EMPTY (-1)
106
107 /* Row and column status */
108 #define ALIVE (0)
109 #define DEAD (-1)
110
111 /* Column status */
112 #define DEAD_PRINCIPAL (-1)
113 #define DEAD_NON_PRINCIPAL (-2)
114
115 /* Macros for row and column status update and checking. */
116 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
117 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
118 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
119 #define COL_IS_DEAD(c) (Col [c].start < ALIVE)
120 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
121 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
122 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
123 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
124 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
125
126 /* ========================================================================== */
127 /* === Colamd reporting mechanism =========================================== */
128 /* ========================================================================== */
129
130 // == Row and Column structures ==
131 template <typename IndexType>
132 struct colamd_col
133 {
134 IndexType start ; /* index for A of first row in this column, or DEAD */
135 /* if column is dead */
136 IndexType length ; /* number of rows in this column */
137 union
138 {
139 IndexType thickness ; /* number of original columns represented by this */
140 /* col, if the column is alive */
141 IndexType parent ; /* parent in parent tree super-column structure, if */
142 /* the column is dead */
143 } shared1 ;
144 union
145 {
146 IndexType score ; /* the score used to maintain heap, if col is alive */
147 IndexType order ; /* pivot ordering of this column, if col is dead */
148 } shared2 ;
149 union
150 {
151 IndexType headhash ; /* head of a hash bucket, if col is at the head of */
152 /* a degree list */
153 IndexType hash ; /* hash value, if col is not in a degree list */
154 IndexType prev ; /* previous column in degree list, if col is in a */
155 /* degree list (but not at the head of a degree list) */
156 } shared3 ;
157 union
158 {
159 IndexType degree_next ; /* next column, if col is in a degree list */
160 IndexType hash_next ; /* next column, if col is in a hash list */
161 } shared4 ;
162
163 };
164
165 template <typename IndexType>
166 struct Colamd_Row
167 {
168 IndexType start ; /* index for A of first col in this row */
169 IndexType length ; /* number of principal columns in this row */
170 union
171 {
172 IndexType degree ; /* number of principal & non-principal columns in row */
173 IndexType p ; /* used as a row pointer in init_rows_cols () */
174 } shared1 ;
175 union
176 {
177 IndexType mark ; /* for computing set differences and marking dead rows*/
178 IndexType first_column ;/* first column in row (used in garbage collection) */
179 } shared2 ;
180
181 };
182
183 /* ========================================================================== */
184 /* === Colamd recommended memory size ======================================= */
185 /* ========================================================================== */
186
187 /*
188 The recommended length Alen of the array A passed to colamd is given by
189 the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any
190 argument is negative. 2*nnz space is required for the row and column
191 indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is
192 required for the Col and Row arrays, respectively, which are internal to
193 colamd. An additional n_col space is the minimal amount of "elbow room",
194 and nnz/5 more space is recommended for run time efficiency.
195
196 This macro is not needed when using symamd.
197
198 Explicit typecast to IndexType added Sept. 23, 2002, COLAMD version 2.2, to avoid
199 gcc -pedantic warning messages.
200 */
201 template <typename IndexType>
colamd_c(IndexType n_col)202 inline IndexType colamd_c(IndexType n_col)
203 { return IndexType( ((n_col) + 1) * sizeof (colamd_col<IndexType>) / sizeof (IndexType) ) ; }
204
205 template <typename IndexType>
colamd_r(IndexType n_row)206 inline IndexType colamd_r(IndexType n_row)
207 { return IndexType(((n_row) + 1) * sizeof (Colamd_Row<IndexType>) / sizeof (IndexType)); }
208
209 // Prototypes of non-user callable routines
210 template <typename IndexType>
211 static IndexType init_rows_cols (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> col [], IndexType A [], IndexType p [], IndexType stats[COLAMD_STATS] );
212
213 template <typename IndexType>
214 static void init_scoring (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType head [], double knobs[COLAMD_KNOBS], IndexType *p_n_row2, IndexType *p_n_col2, IndexType *p_max_deg);
215
216 template <typename IndexType>
217 static IndexType find_ordering (IndexType n_row, IndexType n_col, IndexType Alen, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType n_col2, IndexType max_deg, IndexType pfree);
218
219 template <typename IndexType>
220 static void order_children (IndexType n_col, colamd_col<IndexType> Col [], IndexType p []);
221
222 template <typename IndexType>
223 static void detect_super_cols (colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType row_start, IndexType row_length ) ;
224
225 template <typename IndexType>
226 static IndexType garbage_collection (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType *pfree) ;
227
228 template <typename IndexType>
229 static inline IndexType clear_mark (IndexType n_row, Colamd_Row<IndexType> Row [] ) ;
230
231 /* === No debugging ========================================================= */
232
233 #define COLAMD_DEBUG0(params) ;
234 #define COLAMD_DEBUG1(params) ;
235 #define COLAMD_DEBUG2(params) ;
236 #define COLAMD_DEBUG3(params) ;
237 #define COLAMD_DEBUG4(params) ;
238
239 #define COLAMD_ASSERT(expression) ((void) 0)
240
241
242 /**
243 * \brief Returns the recommended value of Alen
244 *
245 * Returns recommended value of Alen for use by colamd.
246 * Returns -1 if any input argument is negative.
247 * The use of this routine or macro is optional.
248 * Note that the macro uses its arguments more than once,
249 * so be careful for side effects, if you pass expressions as arguments to COLAMD_RECOMMENDED.
250 *
251 * \param nnz nonzeros in A
252 * \param n_row number of rows in A
253 * \param n_col number of columns in A
254 * \return recommended value of Alen for use by colamd
255 */
256 template <typename IndexType>
colamd_recommended(IndexType nnz,IndexType n_row,IndexType n_col)257 inline IndexType colamd_recommended ( IndexType nnz, IndexType n_row, IndexType n_col)
258 {
259 if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0)
260 return (-1);
261 else
262 return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5));
263 }
264
265 /**
266 * \brief set default parameters The use of this routine is optional.
267 *
268 * Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col)
269 * entries are removed prior to ordering. Columns with more than
270 * (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to
271 * ordering, and placed last in the output column ordering.
272 *
273 * COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
274 * respectively, in colamd.h. Default values of these two knobs
275 * are both 0.5. Currently, only knobs [0] and knobs [1] are
276 * used, but future versions may use more knobs. If so, they will
277 * be properly set to their defaults by the future version of
278 * colamd_set_defaults, so that the code that calls colamd will
279 * not need to change, assuming that you either use
280 * colamd_set_defaults, or pass a (double *) NULL pointer as the
281 * knobs array to colamd or symamd.
282 *
283 * \param knobs parameter settings for colamd
284 */
285
colamd_set_defaults(double knobs[COLAMD_KNOBS])286 static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS])
287 {
288 /* === Local variables ================================================== */
289
290 int i ;
291
292 if (!knobs)
293 {
294 return ; /* no knobs to initialize */
295 }
296 for (i = 0 ; i < COLAMD_KNOBS ; i++)
297 {
298 knobs [i] = 0 ;
299 }
300 knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */
301 knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */
302 }
303
304 /**
305 * \brief Computes a column ordering using the column approximate minimum degree ordering
306 *
307 * Computes a column ordering (Q) of A such that P(AQ)=LU or
308 * (AQ)'AQ=LL' have less fill-in and require fewer floating point
309 * operations than factorizing the unpermuted matrix A or A'A,
310 * respectively.
311 *
312 *
313 * \param n_row number of rows in A
314 * \param n_col number of columns in A
315 * \param Alen, size of the array A
316 * \param A row indices of the matrix, of size ALen
317 * \param p column pointers of A, of size n_col+1
318 * \param knobs parameter settings for colamd
319 * \param stats colamd output statistics and error codes
320 */
321 template <typename IndexType>
colamd(IndexType n_row,IndexType n_col,IndexType Alen,IndexType * A,IndexType * p,double knobs[COLAMD_KNOBS],IndexType stats[COLAMD_STATS])322 static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS])
323 {
324 /* === Local variables ================================================== */
325
326 IndexType i ; /* loop index */
327 IndexType nnz ; /* nonzeros in A */
328 IndexType Row_size ; /* size of Row [], in integers */
329 IndexType Col_size ; /* size of Col [], in integers */
330 IndexType need ; /* minimum required length of A */
331 Colamd_Row<IndexType> *Row ; /* pointer into A of Row [0..n_row] array */
332 colamd_col<IndexType> *Col ; /* pointer into A of Col [0..n_col] array */
333 IndexType n_col2 ; /* number of non-dense, non-empty columns */
334 IndexType n_row2 ; /* number of non-dense, non-empty rows */
335 IndexType ngarbage ; /* number of garbage collections performed */
336 IndexType max_deg ; /* maximum row degree */
337 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
338
339
340 /* === Check the input arguments ======================================== */
341
342 if (!stats)
343 {
344 COLAMD_DEBUG0 (("colamd: stats not present\n")) ;
345 return (false) ;
346 }
347 for (i = 0 ; i < COLAMD_STATS ; i++)
348 {
349 stats [i] = 0 ;
350 }
351 stats [COLAMD_STATUS] = COLAMD_OK ;
352 stats [COLAMD_INFO1] = -1 ;
353 stats [COLAMD_INFO2] = -1 ;
354
355 if (!A) /* A is not present */
356 {
357 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
358 COLAMD_DEBUG0 (("colamd: A not present\n")) ;
359 return (false) ;
360 }
361
362 if (!p) /* p is not present */
363 {
364 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
365 COLAMD_DEBUG0 (("colamd: p not present\n")) ;
366 return (false) ;
367 }
368
369 if (n_row < 0) /* n_row must be >= 0 */
370 {
371 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
372 stats [COLAMD_INFO1] = n_row ;
373 COLAMD_DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
374 return (false) ;
375 }
376
377 if (n_col < 0) /* n_col must be >= 0 */
378 {
379 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
380 stats [COLAMD_INFO1] = n_col ;
381 COLAMD_DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
382 return (false) ;
383 }
384
385 nnz = p [n_col] ;
386 if (nnz < 0) /* nnz must be >= 0 */
387 {
388 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
389 stats [COLAMD_INFO1] = nnz ;
390 COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
391 return (false) ;
392 }
393
394 if (p [0] != 0)
395 {
396 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
397 stats [COLAMD_INFO1] = p [0] ;
398 COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
399 return (false) ;
400 }
401
402 /* === If no knobs, set default knobs =================================== */
403
404 if (!knobs)
405 {
406 colamd_set_defaults (default_knobs) ;
407 knobs = default_knobs ;
408 }
409
410 /* === Allocate the Row and Col arrays from array A ===================== */
411
412 Col_size = colamd_c (n_col) ;
413 Row_size = colamd_r (n_row) ;
414 need = 2*nnz + n_col + Col_size + Row_size ;
415
416 if (need > Alen)
417 {
418 /* not enough space in array A to perform the ordering */
419 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
420 stats [COLAMD_INFO1] = need ;
421 stats [COLAMD_INFO2] = Alen ;
422 COLAMD_DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
423 return (false) ;
424 }
425
426 Alen -= Col_size + Row_size ;
427 Col = (colamd_col<IndexType> *) &A [Alen] ;
428 Row = (Colamd_Row<IndexType> *) &A [Alen + Col_size] ;
429
430 /* === Construct the row and column data structures ===================== */
431
432 if (!Eigen::internal::init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
433 {
434 /* input matrix is invalid */
435 COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ;
436 return (false) ;
437 }
438
439 /* === Initialize scores, kill dense rows/columns ======================= */
440
441 Eigen::internal::init_scoring (n_row, n_col, Row, Col, A, p, knobs,
442 &n_row2, &n_col2, &max_deg) ;
443
444 /* === Order the supercolumns =========================================== */
445
446 ngarbage = Eigen::internal::find_ordering (n_row, n_col, Alen, Row, Col, A, p,
447 n_col2, max_deg, 2*nnz) ;
448
449 /* === Order the non-principal columns ================================== */
450
451 Eigen::internal::order_children (n_col, Col, p) ;
452
453 /* === Return statistics in stats ======================================= */
454
455 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
456 stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
457 stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
458 COLAMD_DEBUG0 (("colamd: done.\n")) ;
459 return (true) ;
460 }
461
462 /* ========================================================================== */
463 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
464 /* ========================================================================== */
465
466 /* There are no user-callable routines beyond this point in the file */
467
468
469 /* ========================================================================== */
470 /* === init_rows_cols ======================================================= */
471 /* ========================================================================== */
472
473 /*
474 Takes the column form of the matrix in A and creates the row form of the
475 matrix. Also, row and column attributes are stored in the Col and Row
476 structs. If the columns are un-sorted or contain duplicate row indices,
477 this routine will also sort and remove duplicate row indices from the
478 column form of the matrix. Returns false if the matrix is invalid,
479 true otherwise. Not user-callable.
480 */
481 template <typename IndexType>
init_rows_cols(IndexType n_row,IndexType n_col,Colamd_Row<IndexType> Row[],colamd_col<IndexType> Col[],IndexType A[],IndexType p[],IndexType stats[COLAMD_STATS])482 static IndexType init_rows_cols /* returns true if OK, or false otherwise */
483 (
484 /* === Parameters ======================================================= */
485
486 IndexType n_row, /* number of rows of A */
487 IndexType n_col, /* number of columns of A */
488 Colamd_Row<IndexType> Row [], /* of size n_row+1 */
489 colamd_col<IndexType> Col [], /* of size n_col+1 */
490 IndexType A [], /* row indices of A, of size Alen */
491 IndexType p [], /* pointers to columns in A, of size n_col+1 */
492 IndexType stats [COLAMD_STATS] /* colamd statistics */
493 )
494 {
495 /* === Local variables ================================================== */
496
497 IndexType col ; /* a column index */
498 IndexType row ; /* a row index */
499 IndexType *cp ; /* a column pointer */
500 IndexType *cp_end ; /* a pointer to the end of a column */
501 IndexType *rp ; /* a row pointer */
502 IndexType *rp_end ; /* a pointer to the end of a row */
503 IndexType last_row ; /* previous row */
504
505 /* === Initialize columns, and check column pointers ==================== */
506
507 for (col = 0 ; col < n_col ; col++)
508 {
509 Col [col].start = p [col] ;
510 Col [col].length = p [col+1] - p [col] ;
511
512 if ((Col [col].length) < 0) // extra parentheses to work-around gcc bug 10200
513 {
514 /* column pointers must be non-decreasing */
515 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
516 stats [COLAMD_INFO1] = col ;
517 stats [COLAMD_INFO2] = Col [col].length ;
518 COLAMD_DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
519 return (false) ;
520 }
521
522 Col [col].shared1.thickness = 1 ;
523 Col [col].shared2.score = 0 ;
524 Col [col].shared3.prev = COLAMD_EMPTY ;
525 Col [col].shared4.degree_next = COLAMD_EMPTY ;
526 }
527
528 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
529
530 /* === Scan columns, compute row degrees, and check row indices ========= */
531
532 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
533
534 for (row = 0 ; row < n_row ; row++)
535 {
536 Row [row].length = 0 ;
537 Row [row].shared2.mark = -1 ;
538 }
539
540 for (col = 0 ; col < n_col ; col++)
541 {
542 last_row = -1 ;
543
544 cp = &A [p [col]] ;
545 cp_end = &A [p [col+1]] ;
546
547 while (cp < cp_end)
548 {
549 row = *cp++ ;
550
551 /* make sure row indices within range */
552 if (row < 0 || row >= n_row)
553 {
554 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
555 stats [COLAMD_INFO1] = col ;
556 stats [COLAMD_INFO2] = row ;
557 stats [COLAMD_INFO3] = n_row ;
558 COLAMD_DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
559 return (false) ;
560 }
561
562 if (row <= last_row || Row [row].shared2.mark == col)
563 {
564 /* row index are unsorted or repeated (or both), thus col */
565 /* is jumbled. This is a notice, not an error condition. */
566 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
567 stats [COLAMD_INFO1] = col ;
568 stats [COLAMD_INFO2] = row ;
569 (stats [COLAMD_INFO3]) ++ ;
570 COLAMD_DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
571 }
572
573 if (Row [row].shared2.mark != col)
574 {
575 Row [row].length++ ;
576 }
577 else
578 {
579 /* this is a repeated entry in the column, */
580 /* it will be removed */
581 Col [col].length-- ;
582 }
583
584 /* mark the row as having been seen in this column */
585 Row [row].shared2.mark = col ;
586
587 last_row = row ;
588 }
589 }
590
591 /* === Compute row pointers ============================================= */
592
593 /* row form of the matrix starts directly after the column */
594 /* form of matrix in A */
595 Row [0].start = p [n_col] ;
596 Row [0].shared1.p = Row [0].start ;
597 Row [0].shared2.mark = -1 ;
598 for (row = 1 ; row < n_row ; row++)
599 {
600 Row [row].start = Row [row-1].start + Row [row-1].length ;
601 Row [row].shared1.p = Row [row].start ;
602 Row [row].shared2.mark = -1 ;
603 }
604
605 /* === Create row form ================================================== */
606
607 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
608 {
609 /* if cols jumbled, watch for repeated row indices */
610 for (col = 0 ; col < n_col ; col++)
611 {
612 cp = &A [p [col]] ;
613 cp_end = &A [p [col+1]] ;
614 while (cp < cp_end)
615 {
616 row = *cp++ ;
617 if (Row [row].shared2.mark != col)
618 {
619 A [(Row [row].shared1.p)++] = col ;
620 Row [row].shared2.mark = col ;
621 }
622 }
623 }
624 }
625 else
626 {
627 /* if cols not jumbled, we don't need the mark (this is faster) */
628 for (col = 0 ; col < n_col ; col++)
629 {
630 cp = &A [p [col]] ;
631 cp_end = &A [p [col+1]] ;
632 while (cp < cp_end)
633 {
634 A [(Row [*cp++].shared1.p)++] = col ;
635 }
636 }
637 }
638
639 /* === Clear the row marks and set row degrees ========================== */
640
641 for (row = 0 ; row < n_row ; row++)
642 {
643 Row [row].shared2.mark = 0 ;
644 Row [row].shared1.degree = Row [row].length ;
645 }
646
647 /* === See if we need to re-create columns ============================== */
648
649 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
650 {
651 COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
652
653
654 /* === Compute col pointers ========================================= */
655
656 /* col form of the matrix starts at A [0]. */
657 /* Note, we may have a gap between the col form and the row */
658 /* form if there were duplicate entries, if so, it will be */
659 /* removed upon the first garbage collection */
660 Col [0].start = 0 ;
661 p [0] = Col [0].start ;
662 for (col = 1 ; col < n_col ; col++)
663 {
664 /* note that the lengths here are for pruned columns, i.e. */
665 /* no duplicate row indices will exist for these columns */
666 Col [col].start = Col [col-1].start + Col [col-1].length ;
667 p [col] = Col [col].start ;
668 }
669
670 /* === Re-create col form =========================================== */
671
672 for (row = 0 ; row < n_row ; row++)
673 {
674 rp = &A [Row [row].start] ;
675 rp_end = rp + Row [row].length ;
676 while (rp < rp_end)
677 {
678 A [(p [*rp++])++] = row ;
679 }
680 }
681 }
682
683 /* === Done. Matrix is not (or no longer) jumbled ====================== */
684
685 return (true) ;
686 }
687
688
689 /* ========================================================================== */
690 /* === init_scoring ========================================================= */
691 /* ========================================================================== */
692
693 /*
694 Kills dense or empty columns and rows, calculates an initial score for
695 each column, and places all columns in the degree lists. Not user-callable.
696 */
697 template <typename IndexType>
init_scoring(IndexType n_row,IndexType n_col,Colamd_Row<IndexType> Row[],colamd_col<IndexType> Col[],IndexType A[],IndexType head[],double knobs[COLAMD_KNOBS],IndexType * p_n_row2,IndexType * p_n_col2,IndexType * p_max_deg)698 static void init_scoring
699 (
700 /* === Parameters ======================================================= */
701
702 IndexType n_row, /* number of rows of A */
703 IndexType n_col, /* number of columns of A */
704 Colamd_Row<IndexType> Row [], /* of size n_row+1 */
705 colamd_col<IndexType> Col [], /* of size n_col+1 */
706 IndexType A [], /* column form and row form of A */
707 IndexType head [], /* of size n_col+1 */
708 double knobs [COLAMD_KNOBS],/* parameters */
709 IndexType *p_n_row2, /* number of non-dense, non-empty rows */
710 IndexType *p_n_col2, /* number of non-dense, non-empty columns */
711 IndexType *p_max_deg /* maximum row degree */
712 )
713 {
714 /* === Local variables ================================================== */
715
716 IndexType c ; /* a column index */
717 IndexType r, row ; /* a row index */
718 IndexType *cp ; /* a column pointer */
719 IndexType deg ; /* degree of a row or column */
720 IndexType *cp_end ; /* a pointer to the end of a column */
721 IndexType *new_cp ; /* new column pointer */
722 IndexType col_length ; /* length of pruned column */
723 IndexType score ; /* current column score */
724 IndexType n_col2 ; /* number of non-dense, non-empty columns */
725 IndexType n_row2 ; /* number of non-dense, non-empty rows */
726 IndexType dense_row_count ; /* remove rows with more entries than this */
727 IndexType dense_col_count ; /* remove cols with more entries than this */
728 IndexType min_score ; /* smallest column score */
729 IndexType max_deg ; /* maximum row degree */
730 IndexType next_col ; /* Used to add to degree list.*/
731
732
733 /* === Extract knobs ==================================================== */
734
735 dense_row_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_ROW] * n_col), n_col)) ;
736 dense_col_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_COL] * n_row), n_row)) ;
737 COLAMD_DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
738 max_deg = 0 ;
739 n_col2 = n_col ;
740 n_row2 = n_row ;
741
742 /* === Kill empty columns =============================================== */
743
744 /* Put the empty columns at the end in their natural order, so that LU */
745 /* factorization can proceed as far as possible. */
746 for (c = n_col-1 ; c >= 0 ; c--)
747 {
748 deg = Col [c].length ;
749 if (deg == 0)
750 {
751 /* this is a empty column, kill and order it last */
752 Col [c].shared2.order = --n_col2 ;
753 KILL_PRINCIPAL_COL (c) ;
754 }
755 }
756 COLAMD_DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
757
758 /* === Kill dense columns =============================================== */
759
760 /* Put the dense columns at the end, in their natural order */
761 for (c = n_col-1 ; c >= 0 ; c--)
762 {
763 /* skip any dead columns */
764 if (COL_IS_DEAD (c))
765 {
766 continue ;
767 }
768 deg = Col [c].length ;
769 if (deg > dense_col_count)
770 {
771 /* this is a dense column, kill and order it last */
772 Col [c].shared2.order = --n_col2 ;
773 /* decrement the row degrees */
774 cp = &A [Col [c].start] ;
775 cp_end = cp + Col [c].length ;
776 while (cp < cp_end)
777 {
778 Row [*cp++].shared1.degree-- ;
779 }
780 KILL_PRINCIPAL_COL (c) ;
781 }
782 }
783 COLAMD_DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
784
785 /* === Kill dense and empty rows ======================================== */
786
787 for (r = 0 ; r < n_row ; r++)
788 {
789 deg = Row [r].shared1.degree ;
790 COLAMD_ASSERT (deg >= 0 && deg <= n_col) ;
791 if (deg > dense_row_count || deg == 0)
792 {
793 /* kill a dense or empty row */
794 KILL_ROW (r) ;
795 --n_row2 ;
796 }
797 else
798 {
799 /* keep track of max degree of remaining rows */
800 max_deg = numext::maxi(max_deg, deg) ;
801 }
802 }
803 COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
804
805 /* === Compute initial column scores ==================================== */
806
807 /* At this point the row degrees are accurate. They reflect the number */
808 /* of "live" (non-dense) columns in each row. No empty rows exist. */
809 /* Some "live" columns may contain only dead rows, however. These are */
810 /* pruned in the code below. */
811
812 /* now find the initial matlab score for each column */
813 for (c = n_col-1 ; c >= 0 ; c--)
814 {
815 /* skip dead column */
816 if (COL_IS_DEAD (c))
817 {
818 continue ;
819 }
820 score = 0 ;
821 cp = &A [Col [c].start] ;
822 new_cp = cp ;
823 cp_end = cp + Col [c].length ;
824 while (cp < cp_end)
825 {
826 /* get a row */
827 row = *cp++ ;
828 /* skip if dead */
829 if (ROW_IS_DEAD (row))
830 {
831 continue ;
832 }
833 /* compact the column */
834 *new_cp++ = row ;
835 /* add row's external degree */
836 score += Row [row].shared1.degree - 1 ;
837 /* guard against integer overflow */
838 score = numext::mini(score, n_col) ;
839 }
840 /* determine pruned column length */
841 col_length = (IndexType) (new_cp - &A [Col [c].start]) ;
842 if (col_length == 0)
843 {
844 /* a newly-made null column (all rows in this col are "dense" */
845 /* and have already been killed) */
846 COLAMD_DEBUG2 (("Newly null killed: %d\n", c)) ;
847 Col [c].shared2.order = --n_col2 ;
848 KILL_PRINCIPAL_COL (c) ;
849 }
850 else
851 {
852 /* set column length and set score */
853 COLAMD_ASSERT (score >= 0) ;
854 COLAMD_ASSERT (score <= n_col) ;
855 Col [c].length = col_length ;
856 Col [c].shared2.score = score ;
857 }
858 }
859 COLAMD_DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
860 n_col-n_col2)) ;
861
862 /* At this point, all empty rows and columns are dead. All live columns */
863 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
864 /* yet). Rows may contain dead columns, but all live rows contain at */
865 /* least one live column. */
866
867 /* === Initialize degree lists ========================================== */
868
869
870 /* clear the hash buckets */
871 for (c = 0 ; c <= n_col ; c++)
872 {
873 head [c] = COLAMD_EMPTY ;
874 }
875 min_score = n_col ;
876 /* place in reverse order, so low column indices are at the front */
877 /* of the lists. This is to encourage natural tie-breaking */
878 for (c = n_col-1 ; c >= 0 ; c--)
879 {
880 /* only add principal columns to degree lists */
881 if (COL_IS_ALIVE (c))
882 {
883 COLAMD_DEBUG4 (("place %d score %d minscore %d ncol %d\n",
884 c, Col [c].shared2.score, min_score, n_col)) ;
885
886 /* === Add columns score to DList =============================== */
887
888 score = Col [c].shared2.score ;
889
890 COLAMD_ASSERT (min_score >= 0) ;
891 COLAMD_ASSERT (min_score <= n_col) ;
892 COLAMD_ASSERT (score >= 0) ;
893 COLAMD_ASSERT (score <= n_col) ;
894 COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ;
895
896 /* now add this column to dList at proper score location */
897 next_col = head [score] ;
898 Col [c].shared3.prev = COLAMD_EMPTY ;
899 Col [c].shared4.degree_next = next_col ;
900
901 /* if there already was a column with the same score, set its */
902 /* previous pointer to this new column */
903 if (next_col != COLAMD_EMPTY)
904 {
905 Col [next_col].shared3.prev = c ;
906 }
907 head [score] = c ;
908
909 /* see if this score is less than current min */
910 min_score = numext::mini(min_score, score) ;
911
912
913 }
914 }
915
916
917 /* === Return number of remaining columns, and max row degree =========== */
918
919 *p_n_col2 = n_col2 ;
920 *p_n_row2 = n_row2 ;
921 *p_max_deg = max_deg ;
922 }
923
924
925 /* ========================================================================== */
926 /* === find_ordering ======================================================== */
927 /* ========================================================================== */
928
929 /*
930 Order the principal columns of the supercolumn form of the matrix
931 (no supercolumns on input). Uses a minimum approximate column minimum
932 degree ordering method. Not user-callable.
933 */
934 template <typename IndexType>
find_ordering(IndexType n_row,IndexType n_col,IndexType Alen,Colamd_Row<IndexType> Row[],colamd_col<IndexType> Col[],IndexType A[],IndexType head[],IndexType n_col2,IndexType max_deg,IndexType pfree)935 static IndexType find_ordering /* return the number of garbage collections */
936 (
937 /* === Parameters ======================================================= */
938
939 IndexType n_row, /* number of rows of A */
940 IndexType n_col, /* number of columns of A */
941 IndexType Alen, /* size of A, 2*nnz + n_col or larger */
942 Colamd_Row<IndexType> Row [], /* of size n_row+1 */
943 colamd_col<IndexType> Col [], /* of size n_col+1 */
944 IndexType A [], /* column form and row form of A */
945 IndexType head [], /* of size n_col+1 */
946 IndexType n_col2, /* Remaining columns to order */
947 IndexType max_deg, /* Maximum row degree */
948 IndexType pfree /* index of first free slot (2*nnz on entry) */
949 )
950 {
951 /* === Local variables ================================================== */
952
953 IndexType k ; /* current pivot ordering step */
954 IndexType pivot_col ; /* current pivot column */
955 IndexType *cp ; /* a column pointer */
956 IndexType *rp ; /* a row pointer */
957 IndexType pivot_row ; /* current pivot row */
958 IndexType *new_cp ; /* modified column pointer */
959 IndexType *new_rp ; /* modified row pointer */
960 IndexType pivot_row_start ; /* pointer to start of pivot row */
961 IndexType pivot_row_degree ; /* number of columns in pivot row */
962 IndexType pivot_row_length ; /* number of supercolumns in pivot row */
963 IndexType pivot_col_score ; /* score of pivot column */
964 IndexType needed_memory ; /* free space needed for pivot row */
965 IndexType *cp_end ; /* pointer to the end of a column */
966 IndexType *rp_end ; /* pointer to the end of a row */
967 IndexType row ; /* a row index */
968 IndexType col ; /* a column index */
969 IndexType max_score ; /* maximum possible score */
970 IndexType cur_score ; /* score of current column */
971 unsigned int hash ; /* hash value for supernode detection */
972 IndexType head_column ; /* head of hash bucket */
973 IndexType first_col ; /* first column in hash bucket */
974 IndexType tag_mark ; /* marker value for mark array */
975 IndexType row_mark ; /* Row [row].shared2.mark */
976 IndexType set_difference ; /* set difference size of row with pivot row */
977 IndexType min_score ; /* smallest column score */
978 IndexType col_thickness ; /* "thickness" (no. of columns in a supercol) */
979 IndexType max_mark ; /* maximum value of tag_mark */
980 IndexType pivot_col_thickness ; /* number of columns represented by pivot col */
981 IndexType prev_col ; /* Used by Dlist operations. */
982 IndexType next_col ; /* Used by Dlist operations. */
983 IndexType ngarbage ; /* number of garbage collections performed */
984
985
986 /* === Initialization and clear mark ==================================== */
987
988 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */
989 tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
990 min_score = 0 ;
991 ngarbage = 0 ;
992 COLAMD_DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
993
994 /* === Order the columns ================================================ */
995
996 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
997 {
998
999 /* === Select pivot column, and order it ============================ */
1000
1001 /* make sure degree list isn't empty */
1002 COLAMD_ASSERT (min_score >= 0) ;
1003 COLAMD_ASSERT (min_score <= n_col) ;
1004 COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ;
1005
1006 /* get pivot column from head of minimum degree list */
1007 while (min_score < n_col && head [min_score] == COLAMD_EMPTY)
1008 {
1009 min_score++ ;
1010 }
1011 pivot_col = head [min_score] ;
1012 COLAMD_ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
1013 next_col = Col [pivot_col].shared4.degree_next ;
1014 head [min_score] = next_col ;
1015 if (next_col != COLAMD_EMPTY)
1016 {
1017 Col [next_col].shared3.prev = COLAMD_EMPTY ;
1018 }
1019
1020 COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ;
1021 COLAMD_DEBUG3 (("Pivot col: %d\n", pivot_col)) ;
1022
1023 /* remember score for defrag check */
1024 pivot_col_score = Col [pivot_col].shared2.score ;
1025
1026 /* the pivot column is the kth column in the pivot order */
1027 Col [pivot_col].shared2.order = k ;
1028
1029 /* increment order count by column thickness */
1030 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
1031 k += pivot_col_thickness ;
1032 COLAMD_ASSERT (pivot_col_thickness > 0) ;
1033
1034 /* === Garbage_collection, if necessary ============================= */
1035
1036 needed_memory = numext::mini(pivot_col_score, n_col - k) ;
1037 if (pfree + needed_memory >= Alen)
1038 {
1039 pfree = Eigen::internal::garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
1040 ngarbage++ ;
1041 /* after garbage collection we will have enough */
1042 COLAMD_ASSERT (pfree + needed_memory < Alen) ;
1043 /* garbage collection has wiped out the Row[].shared2.mark array */
1044 tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
1045
1046 }
1047
1048 /* === Compute pivot row pattern ==================================== */
1049
1050 /* get starting location for this new merged row */
1051 pivot_row_start = pfree ;
1052
1053 /* initialize new row counts to zero */
1054 pivot_row_degree = 0 ;
1055
1056 /* tag pivot column as having been visited so it isn't included */
1057 /* in merged pivot row */
1058 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
1059
1060 /* pivot row is the union of all rows in the pivot column pattern */
1061 cp = &A [Col [pivot_col].start] ;
1062 cp_end = cp + Col [pivot_col].length ;
1063 while (cp < cp_end)
1064 {
1065 /* get a row */
1066 row = *cp++ ;
1067 COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
1068 /* skip if row is dead */
1069 if (ROW_IS_DEAD (row))
1070 {
1071 continue ;
1072 }
1073 rp = &A [Row [row].start] ;
1074 rp_end = rp + Row [row].length ;
1075 while (rp < rp_end)
1076 {
1077 /* get a column */
1078 col = *rp++ ;
1079 /* add the column, if alive and untagged */
1080 col_thickness = Col [col].shared1.thickness ;
1081 if (col_thickness > 0 && COL_IS_ALIVE (col))
1082 {
1083 /* tag column in pivot row */
1084 Col [col].shared1.thickness = -col_thickness ;
1085 COLAMD_ASSERT (pfree < Alen) ;
1086 /* place column in pivot row */
1087 A [pfree++] = col ;
1088 pivot_row_degree += col_thickness ;
1089 }
1090 }
1091 }
1092
1093 /* clear tag on pivot column */
1094 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
1095 max_deg = numext::maxi(max_deg, pivot_row_degree) ;
1096
1097
1098 /* === Kill all rows used to construct pivot row ==================== */
1099
1100 /* also kill pivot row, temporarily */
1101 cp = &A [Col [pivot_col].start] ;
1102 cp_end = cp + Col [pivot_col].length ;
1103 while (cp < cp_end)
1104 {
1105 /* may be killing an already dead row */
1106 row = *cp++ ;
1107 COLAMD_DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
1108 KILL_ROW (row) ;
1109 }
1110
1111 /* === Select a row index to use as the new pivot row =============== */
1112
1113 pivot_row_length = pfree - pivot_row_start ;
1114 if (pivot_row_length > 0)
1115 {
1116 /* pick the "pivot" row arbitrarily (first row in col) */
1117 pivot_row = A [Col [pivot_col].start] ;
1118 COLAMD_DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
1119 }
1120 else
1121 {
1122 /* there is no pivot row, since it is of zero length */
1123 pivot_row = COLAMD_EMPTY ;
1124 COLAMD_ASSERT (pivot_row_length == 0) ;
1125 }
1126 COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
1127
1128 /* === Approximate degree computation =============================== */
1129
1130 /* Here begins the computation of the approximate degree. The column */
1131 /* score is the sum of the pivot row "length", plus the size of the */
1132 /* set differences of each row in the column minus the pattern of the */
1133 /* pivot row itself. The column ("thickness") itself is also */
1134 /* excluded from the column score (we thus use an approximate */
1135 /* external degree). */
1136
1137 /* The time taken by the following code (compute set differences, and */
1138 /* add them up) is proportional to the size of the data structure */
1139 /* being scanned - that is, the sum of the sizes of each column in */
1140 /* the pivot row. Thus, the amortized time to compute a column score */
1141 /* is proportional to the size of that column (where size, in this */
1142 /* context, is the column "length", or the number of row indices */
1143 /* in that column). The number of row indices in a column is */
1144 /* monotonically non-decreasing, from the length of the original */
1145 /* column on input to colamd. */
1146
1147 /* === Compute set differences ====================================== */
1148
1149 COLAMD_DEBUG3 (("** Computing set differences phase. **\n")) ;
1150
1151 /* pivot row is currently dead - it will be revived later. */
1152
1153 COLAMD_DEBUG3 (("Pivot row: ")) ;
1154 /* for each column in pivot row */
1155 rp = &A [pivot_row_start] ;
1156 rp_end = rp + pivot_row_length ;
1157 while (rp < rp_end)
1158 {
1159 col = *rp++ ;
1160 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1161 COLAMD_DEBUG3 (("Col: %d\n", col)) ;
1162
1163 /* clear tags used to construct pivot row pattern */
1164 col_thickness = -Col [col].shared1.thickness ;
1165 COLAMD_ASSERT (col_thickness > 0) ;
1166 Col [col].shared1.thickness = col_thickness ;
1167
1168 /* === Remove column from degree list =========================== */
1169
1170 cur_score = Col [col].shared2.score ;
1171 prev_col = Col [col].shared3.prev ;
1172 next_col = Col [col].shared4.degree_next ;
1173 COLAMD_ASSERT (cur_score >= 0) ;
1174 COLAMD_ASSERT (cur_score <= n_col) ;
1175 COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ;
1176 if (prev_col == COLAMD_EMPTY)
1177 {
1178 head [cur_score] = next_col ;
1179 }
1180 else
1181 {
1182 Col [prev_col].shared4.degree_next = next_col ;
1183 }
1184 if (next_col != COLAMD_EMPTY)
1185 {
1186 Col [next_col].shared3.prev = prev_col ;
1187 }
1188
1189 /* === Scan the column ========================================== */
1190
1191 cp = &A [Col [col].start] ;
1192 cp_end = cp + Col [col].length ;
1193 while (cp < cp_end)
1194 {
1195 /* get a row */
1196 row = *cp++ ;
1197 row_mark = Row [row].shared2.mark ;
1198 /* skip if dead */
1199 if (ROW_IS_MARKED_DEAD (row_mark))
1200 {
1201 continue ;
1202 }
1203 COLAMD_ASSERT (row != pivot_row) ;
1204 set_difference = row_mark - tag_mark ;
1205 /* check if the row has been seen yet */
1206 if (set_difference < 0)
1207 {
1208 COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ;
1209 set_difference = Row [row].shared1.degree ;
1210 }
1211 /* subtract column thickness from this row's set difference */
1212 set_difference -= col_thickness ;
1213 COLAMD_ASSERT (set_difference >= 0) ;
1214 /* absorb this row if the set difference becomes zero */
1215 if (set_difference == 0)
1216 {
1217 COLAMD_DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
1218 KILL_ROW (row) ;
1219 }
1220 else
1221 {
1222 /* save the new mark */
1223 Row [row].shared2.mark = set_difference + tag_mark ;
1224 }
1225 }
1226 }
1227
1228
1229 /* === Add up set differences for each column ======================= */
1230
1231 COLAMD_DEBUG3 (("** Adding set differences phase. **\n")) ;
1232
1233 /* for each column in pivot row */
1234 rp = &A [pivot_row_start] ;
1235 rp_end = rp + pivot_row_length ;
1236 while (rp < rp_end)
1237 {
1238 /* get a column */
1239 col = *rp++ ;
1240 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1241 hash = 0 ;
1242 cur_score = 0 ;
1243 cp = &A [Col [col].start] ;
1244 /* compact the column */
1245 new_cp = cp ;
1246 cp_end = cp + Col [col].length ;
1247
1248 COLAMD_DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
1249
1250 while (cp < cp_end)
1251 {
1252 /* get a row */
1253 row = *cp++ ;
1254 COLAMD_ASSERT(row >= 0 && row < n_row) ;
1255 row_mark = Row [row].shared2.mark ;
1256 /* skip if dead */
1257 if (ROW_IS_MARKED_DEAD (row_mark))
1258 {
1259 continue ;
1260 }
1261 COLAMD_ASSERT (row_mark > tag_mark) ;
1262 /* compact the column */
1263 *new_cp++ = row ;
1264 /* compute hash function */
1265 hash += row ;
1266 /* add set difference */
1267 cur_score += row_mark - tag_mark ;
1268 /* integer overflow... */
1269 cur_score = numext::mini(cur_score, n_col) ;
1270 }
1271
1272 /* recompute the column's length */
1273 Col [col].length = (IndexType) (new_cp - &A [Col [col].start]) ;
1274
1275 /* === Further mass elimination ================================= */
1276
1277 if (Col [col].length == 0)
1278 {
1279 COLAMD_DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
1280 /* nothing left but the pivot row in this column */
1281 KILL_PRINCIPAL_COL (col) ;
1282 pivot_row_degree -= Col [col].shared1.thickness ;
1283 COLAMD_ASSERT (pivot_row_degree >= 0) ;
1284 /* order it */
1285 Col [col].shared2.order = k ;
1286 /* increment order count by column thickness */
1287 k += Col [col].shared1.thickness ;
1288 }
1289 else
1290 {
1291 /* === Prepare for supercolumn detection ==================== */
1292
1293 COLAMD_DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
1294
1295 /* save score so far */
1296 Col [col].shared2.score = cur_score ;
1297
1298 /* add column to hash table, for supercolumn detection */
1299 hash %= n_col + 1 ;
1300
1301 COLAMD_DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
1302 COLAMD_ASSERT (hash <= n_col) ;
1303
1304 head_column = head [hash] ;
1305 if (head_column > COLAMD_EMPTY)
1306 {
1307 /* degree list "hash" is non-empty, use prev (shared3) of */
1308 /* first column in degree list as head of hash bucket */
1309 first_col = Col [head_column].shared3.headhash ;
1310 Col [head_column].shared3.headhash = col ;
1311 }
1312 else
1313 {
1314 /* degree list "hash" is empty, use head as hash bucket */
1315 first_col = - (head_column + 2) ;
1316 head [hash] = - (col + 2) ;
1317 }
1318 Col [col].shared4.hash_next = first_col ;
1319
1320 /* save hash function in Col [col].shared3.hash */
1321 Col [col].shared3.hash = (IndexType) hash ;
1322 COLAMD_ASSERT (COL_IS_ALIVE (col)) ;
1323 }
1324 }
1325
1326 /* The approximate external column degree is now computed. */
1327
1328 /* === Supercolumn detection ======================================== */
1329
1330 COLAMD_DEBUG3 (("** Supercolumn detection phase. **\n")) ;
1331
1332 Eigen::internal::detect_super_cols (Col, A, head, pivot_row_start, pivot_row_length) ;
1333
1334 /* === Kill the pivotal column ====================================== */
1335
1336 KILL_PRINCIPAL_COL (pivot_col) ;
1337
1338 /* === Clear mark =================================================== */
1339
1340 tag_mark += (max_deg + 1) ;
1341 if (tag_mark >= max_mark)
1342 {
1343 COLAMD_DEBUG2 (("clearing tag_mark\n")) ;
1344 tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
1345 }
1346
1347 /* === Finalize the new pivot row, and column scores ================ */
1348
1349 COLAMD_DEBUG3 (("** Finalize scores phase. **\n")) ;
1350
1351 /* for each column in pivot row */
1352 rp = &A [pivot_row_start] ;
1353 /* compact the pivot row */
1354 new_rp = rp ;
1355 rp_end = rp + pivot_row_length ;
1356 while (rp < rp_end)
1357 {
1358 col = *rp++ ;
1359 /* skip dead columns */
1360 if (COL_IS_DEAD (col))
1361 {
1362 continue ;
1363 }
1364 *new_rp++ = col ;
1365 /* add new pivot row to column */
1366 A [Col [col].start + (Col [col].length++)] = pivot_row ;
1367
1368 /* retrieve score so far and add on pivot row's degree. */
1369 /* (we wait until here for this in case the pivot */
1370 /* row's degree was reduced due to mass elimination). */
1371 cur_score = Col [col].shared2.score + pivot_row_degree ;
1372
1373 /* calculate the max possible score as the number of */
1374 /* external columns minus the 'k' value minus the */
1375 /* columns thickness */
1376 max_score = n_col - k - Col [col].shared1.thickness ;
1377
1378 /* make the score the external degree of the union-of-rows */
1379 cur_score -= Col [col].shared1.thickness ;
1380
1381 /* make sure score is less or equal than the max score */
1382 cur_score = numext::mini(cur_score, max_score) ;
1383 COLAMD_ASSERT (cur_score >= 0) ;
1384
1385 /* store updated score */
1386 Col [col].shared2.score = cur_score ;
1387
1388 /* === Place column back in degree list ========================= */
1389
1390 COLAMD_ASSERT (min_score >= 0) ;
1391 COLAMD_ASSERT (min_score <= n_col) ;
1392 COLAMD_ASSERT (cur_score >= 0) ;
1393 COLAMD_ASSERT (cur_score <= n_col) ;
1394 COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ;
1395 next_col = head [cur_score] ;
1396 Col [col].shared4.degree_next = next_col ;
1397 Col [col].shared3.prev = COLAMD_EMPTY ;
1398 if (next_col != COLAMD_EMPTY)
1399 {
1400 Col [next_col].shared3.prev = col ;
1401 }
1402 head [cur_score] = col ;
1403
1404 /* see if this score is less than current min */
1405 min_score = numext::mini(min_score, cur_score) ;
1406
1407 }
1408
1409 /* === Resurrect the new pivot row ================================== */
1410
1411 if (pivot_row_degree > 0)
1412 {
1413 /* update pivot row length to reflect any cols that were killed */
1414 /* during super-col detection and mass elimination */
1415 Row [pivot_row].start = pivot_row_start ;
1416 Row [pivot_row].length = (IndexType) (new_rp - &A[pivot_row_start]) ;
1417 Row [pivot_row].shared1.degree = pivot_row_degree ;
1418 Row [pivot_row].shared2.mark = 0 ;
1419 /* pivot row is no longer dead */
1420 }
1421 }
1422
1423 /* === All principal columns have now been ordered ====================== */
1424
1425 return (ngarbage) ;
1426 }
1427
1428
1429 /* ========================================================================== */
1430 /* === order_children ======================================================= */
1431 /* ========================================================================== */
1432
1433 /*
1434 The find_ordering routine has ordered all of the principal columns (the
1435 representatives of the supercolumns). The non-principal columns have not
1436 yet been ordered. This routine orders those columns by walking up the
1437 parent tree (a column is a child of the column which absorbed it). The
1438 final permutation vector is then placed in p [0 ... n_col-1], with p [0]
1439 being the first column, and p [n_col-1] being the last. It doesn't look
1440 like it at first glance, but be assured that this routine takes time linear
1441 in the number of columns. Although not immediately obvious, the time
1442 taken by this routine is O (n_col), that is, linear in the number of
1443 columns. Not user-callable.
1444 */
1445 template <typename IndexType>
order_children(IndexType n_col,colamd_col<IndexType> Col[],IndexType p[])1446 static inline void order_children
1447 (
1448 /* === Parameters ======================================================= */
1449
1450 IndexType n_col, /* number of columns of A */
1451 colamd_col<IndexType> Col [], /* of size n_col+1 */
1452 IndexType p [] /* p [0 ... n_col-1] is the column permutation*/
1453 )
1454 {
1455 /* === Local variables ================================================== */
1456
1457 IndexType i ; /* loop counter for all columns */
1458 IndexType c ; /* column index */
1459 IndexType parent ; /* index of column's parent */
1460 IndexType order ; /* column's order */
1461
1462 /* === Order each non-principal column ================================== */
1463
1464 for (i = 0 ; i < n_col ; i++)
1465 {
1466 /* find an un-ordered non-principal column */
1467 COLAMD_ASSERT (COL_IS_DEAD (i)) ;
1468 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY)
1469 {
1470 parent = i ;
1471 /* once found, find its principal parent */
1472 do
1473 {
1474 parent = Col [parent].shared1.parent ;
1475 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
1476
1477 /* now, order all un-ordered non-principal columns along path */
1478 /* to this parent. collapse tree at the same time */
1479 c = i ;
1480 /* get order of parent */
1481 order = Col [parent].shared2.order ;
1482
1483 do
1484 {
1485 COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ;
1486
1487 /* order this column */
1488 Col [c].shared2.order = order++ ;
1489 /* collaps tree */
1490 Col [c].shared1.parent = parent ;
1491
1492 /* get immediate parent of this column */
1493 c = Col [c].shared1.parent ;
1494
1495 /* continue until we hit an ordered column. There are */
1496 /* guarranteed not to be anymore unordered columns */
1497 /* above an ordered column */
1498 } while (Col [c].shared2.order == COLAMD_EMPTY) ;
1499
1500 /* re-order the super_col parent to largest order for this group */
1501 Col [parent].shared2.order = order ;
1502 }
1503 }
1504
1505 /* === Generate the permutation ========================================= */
1506
1507 for (c = 0 ; c < n_col ; c++)
1508 {
1509 p [Col [c].shared2.order] = c ;
1510 }
1511 }
1512
1513
1514 /* ========================================================================== */
1515 /* === detect_super_cols ==================================================== */
1516 /* ========================================================================== */
1517
1518 /*
1519 Detects supercolumns by finding matches between columns in the hash buckets.
1520 Check amongst columns in the set A [row_start ... row_start + row_length-1].
1521 The columns under consideration are currently *not* in the degree lists,
1522 and have already been placed in the hash buckets.
1523
1524 The hash bucket for columns whose hash function is equal to h is stored
1525 as follows:
1526
1527 if head [h] is >= 0, then head [h] contains a degree list, so:
1528
1529 head [h] is the first column in degree bucket h.
1530 Col [head [h]].headhash gives the first column in hash bucket h.
1531
1532 otherwise, the degree list is empty, and:
1533
1534 -(head [h] + 2) is the first column in hash bucket h.
1535
1536 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
1537 column" pointer. Col [c].shared3.hash is used instead as the hash number
1538 for that column. The value of Col [c].shared4.hash_next is the next column
1539 in the same hash bucket.
1540
1541 Assuming no, or "few" hash collisions, the time taken by this routine is
1542 linear in the sum of the sizes (lengths) of each column whose score has
1543 just been computed in the approximate degree computation.
1544 Not user-callable.
1545 */
1546 template <typename IndexType>
detect_super_cols(colamd_col<IndexType> Col[],IndexType A[],IndexType head[],IndexType row_start,IndexType row_length)1547 static void detect_super_cols
1548 (
1549 /* === Parameters ======================================================= */
1550
1551 colamd_col<IndexType> Col [], /* of size n_col+1 */
1552 IndexType A [], /* row indices of A */
1553 IndexType head [], /* head of degree lists and hash buckets */
1554 IndexType row_start, /* pointer to set of columns to check */
1555 IndexType row_length /* number of columns to check */
1556 )
1557 {
1558 /* === Local variables ================================================== */
1559
1560 IndexType hash ; /* hash value for a column */
1561 IndexType *rp ; /* pointer to a row */
1562 IndexType c ; /* a column index */
1563 IndexType super_c ; /* column index of the column to absorb into */
1564 IndexType *cp1 ; /* column pointer for column super_c */
1565 IndexType *cp2 ; /* column pointer for column c */
1566 IndexType length ; /* length of column super_c */
1567 IndexType prev_c ; /* column preceding c in hash bucket */
1568 IndexType i ; /* loop counter */
1569 IndexType *rp_end ; /* pointer to the end of the row */
1570 IndexType col ; /* a column index in the row to check */
1571 IndexType head_column ; /* first column in hash bucket or degree list */
1572 IndexType first_col ; /* first column in hash bucket */
1573
1574 /* === Consider each column in the row ================================== */
1575
1576 rp = &A [row_start] ;
1577 rp_end = rp + row_length ;
1578 while (rp < rp_end)
1579 {
1580 col = *rp++ ;
1581 if (COL_IS_DEAD (col))
1582 {
1583 continue ;
1584 }
1585
1586 /* get hash number for this column */
1587 hash = Col [col].shared3.hash ;
1588 COLAMD_ASSERT (hash <= n_col) ;
1589
1590 /* === Get the first column in this hash bucket ===================== */
1591
1592 head_column = head [hash] ;
1593 if (head_column > COLAMD_EMPTY)
1594 {
1595 first_col = Col [head_column].shared3.headhash ;
1596 }
1597 else
1598 {
1599 first_col = - (head_column + 2) ;
1600 }
1601
1602 /* === Consider each column in the hash bucket ====================== */
1603
1604 for (super_c = first_col ; super_c != COLAMD_EMPTY ;
1605 super_c = Col [super_c].shared4.hash_next)
1606 {
1607 COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ;
1608 COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ;
1609 length = Col [super_c].length ;
1610
1611 /* prev_c is the column preceding column c in the hash bucket */
1612 prev_c = super_c ;
1613
1614 /* === Compare super_c with all columns after it ================ */
1615
1616 for (c = Col [super_c].shared4.hash_next ;
1617 c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next)
1618 {
1619 COLAMD_ASSERT (c != super_c) ;
1620 COLAMD_ASSERT (COL_IS_ALIVE (c)) ;
1621 COLAMD_ASSERT (Col [c].shared3.hash == hash) ;
1622
1623 /* not identical if lengths or scores are different */
1624 if (Col [c].length != length ||
1625 Col [c].shared2.score != Col [super_c].shared2.score)
1626 {
1627 prev_c = c ;
1628 continue ;
1629 }
1630
1631 /* compare the two columns */
1632 cp1 = &A [Col [super_c].start] ;
1633 cp2 = &A [Col [c].start] ;
1634
1635 for (i = 0 ; i < length ; i++)
1636 {
1637 /* the columns are "clean" (no dead rows) */
1638 COLAMD_ASSERT (ROW_IS_ALIVE (*cp1)) ;
1639 COLAMD_ASSERT (ROW_IS_ALIVE (*cp2)) ;
1640 /* row indices will same order for both supercols, */
1641 /* no gather scatter nessasary */
1642 if (*cp1++ != *cp2++)
1643 {
1644 break ;
1645 }
1646 }
1647
1648 /* the two columns are different if the for-loop "broke" */
1649 if (i != length)
1650 {
1651 prev_c = c ;
1652 continue ;
1653 }
1654
1655 /* === Got it! two columns are identical =================== */
1656
1657 COLAMD_ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
1658
1659 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
1660 Col [c].shared1.parent = super_c ;
1661 KILL_NON_PRINCIPAL_COL (c) ;
1662 /* order c later, in order_children() */
1663 Col [c].shared2.order = COLAMD_EMPTY ;
1664 /* remove c from hash bucket */
1665 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
1666 }
1667 }
1668
1669 /* === Empty this hash bucket ======================================= */
1670
1671 if (head_column > COLAMD_EMPTY)
1672 {
1673 /* corresponding degree list "hash" is not empty */
1674 Col [head_column].shared3.headhash = COLAMD_EMPTY ;
1675 }
1676 else
1677 {
1678 /* corresponding degree list "hash" is empty */
1679 head [hash] = COLAMD_EMPTY ;
1680 }
1681 }
1682 }
1683
1684
1685 /* ========================================================================== */
1686 /* === garbage_collection =================================================== */
1687 /* ========================================================================== */
1688
1689 /*
1690 Defragments and compacts columns and rows in the workspace A. Used when
1691 all avaliable memory has been used while performing row merging. Returns
1692 the index of the first free position in A, after garbage collection. The
1693 time taken by this routine is linear is the size of the array A, which is
1694 itself linear in the number of nonzeros in the input matrix.
1695 Not user-callable.
1696 */
1697 template <typename IndexType>
garbage_collection(IndexType n_row,IndexType n_col,Colamd_Row<IndexType> Row[],colamd_col<IndexType> Col[],IndexType A[],IndexType * pfree)1698 static IndexType garbage_collection /* returns the new value of pfree */
1699 (
1700 /* === Parameters ======================================================= */
1701
1702 IndexType n_row, /* number of rows */
1703 IndexType n_col, /* number of columns */
1704 Colamd_Row<IndexType> Row [], /* row info */
1705 colamd_col<IndexType> Col [], /* column info */
1706 IndexType A [], /* A [0 ... Alen-1] holds the matrix */
1707 IndexType *pfree /* &A [0] ... pfree is in use */
1708 )
1709 {
1710 /* === Local variables ================================================== */
1711
1712 IndexType *psrc ; /* source pointer */
1713 IndexType *pdest ; /* destination pointer */
1714 IndexType j ; /* counter */
1715 IndexType r ; /* a row index */
1716 IndexType c ; /* a column index */
1717 IndexType length ; /* length of a row or column */
1718
1719 /* === Defragment the columns =========================================== */
1720
1721 pdest = &A[0] ;
1722 for (c = 0 ; c < n_col ; c++)
1723 {
1724 if (COL_IS_ALIVE (c))
1725 {
1726 psrc = &A [Col [c].start] ;
1727
1728 /* move and compact the column */
1729 COLAMD_ASSERT (pdest <= psrc) ;
1730 Col [c].start = (IndexType) (pdest - &A [0]) ;
1731 length = Col [c].length ;
1732 for (j = 0 ; j < length ; j++)
1733 {
1734 r = *psrc++ ;
1735 if (ROW_IS_ALIVE (r))
1736 {
1737 *pdest++ = r ;
1738 }
1739 }
1740 Col [c].length = (IndexType) (pdest - &A [Col [c].start]) ;
1741 }
1742 }
1743
1744 /* === Prepare to defragment the rows =================================== */
1745
1746 for (r = 0 ; r < n_row ; r++)
1747 {
1748 if (ROW_IS_ALIVE (r))
1749 {
1750 if (Row [r].length == 0)
1751 {
1752 /* this row is of zero length. cannot compact it, so kill it */
1753 COLAMD_DEBUG3 (("Defrag row kill\n")) ;
1754 KILL_ROW (r) ;
1755 }
1756 else
1757 {
1758 /* save first column index in Row [r].shared2.first_column */
1759 psrc = &A [Row [r].start] ;
1760 Row [r].shared2.first_column = *psrc ;
1761 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ;
1762 /* flag the start of the row with the one's complement of row */
1763 *psrc = ONES_COMPLEMENT (r) ;
1764
1765 }
1766 }
1767 }
1768
1769 /* === Defragment the rows ============================================== */
1770
1771 psrc = pdest ;
1772 while (psrc < pfree)
1773 {
1774 /* find a negative number ... the start of a row */
1775 if (*psrc++ < 0)
1776 {
1777 psrc-- ;
1778 /* get the row index */
1779 r = ONES_COMPLEMENT (*psrc) ;
1780 COLAMD_ASSERT (r >= 0 && r < n_row) ;
1781 /* restore first column index */
1782 *psrc = Row [r].shared2.first_column ;
1783 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ;
1784
1785 /* move and compact the row */
1786 COLAMD_ASSERT (pdest <= psrc) ;
1787 Row [r].start = (IndexType) (pdest - &A [0]) ;
1788 length = Row [r].length ;
1789 for (j = 0 ; j < length ; j++)
1790 {
1791 c = *psrc++ ;
1792 if (COL_IS_ALIVE (c))
1793 {
1794 *pdest++ = c ;
1795 }
1796 }
1797 Row [r].length = (IndexType) (pdest - &A [Row [r].start]) ;
1798
1799 }
1800 }
1801 /* ensure we found all the rows */
1802 COLAMD_ASSERT (debug_rows == 0) ;
1803
1804 /* === Return the new value of pfree ==================================== */
1805
1806 return ((IndexType) (pdest - &A [0])) ;
1807 }
1808
1809
1810 /* ========================================================================== */
1811 /* === clear_mark =========================================================== */
1812 /* ========================================================================== */
1813
1814 /*
1815 Clears the Row [].shared2.mark array, and returns the new tag_mark.
1816 Return value is the new tag_mark. Not user-callable.
1817 */
1818 template <typename IndexType>
clear_mark(IndexType n_row,Colamd_Row<IndexType> Row[])1819 static inline IndexType clear_mark /* return the new value for tag_mark */
1820 (
1821 /* === Parameters ======================================================= */
1822
1823 IndexType n_row, /* number of rows in A */
1824 Colamd_Row<IndexType> Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
1825 )
1826 {
1827 /* === Local variables ================================================== */
1828
1829 IndexType r ;
1830
1831 for (r = 0 ; r < n_row ; r++)
1832 {
1833 if (ROW_IS_ALIVE (r))
1834 {
1835 Row [r].shared2.mark = 0 ;
1836 }
1837 }
1838 return (1) ;
1839 }
1840
1841
1842 } // namespace internal
1843 #endif
1844