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