Eigen_Colamd.h
Go to the documentation of this file.
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>
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>
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>
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>
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 
256 template <typename IndexType>
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 
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 
321 template <typename IndexType>
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  }
352  stats [COLAMD_INFO1] = -1 ;
353  stats [COLAMD_INFO2] = -1 ;
354 
355  if (!A) /* A is not present */
356  {
358  COLAMD_DEBUG0 (("colamd: A not present\n")) ;
359  return (false) ;
360  }
361 
362  if (!p) /* p is not present */
363  {
365  COLAMD_DEBUG0 (("colamd: p not present\n")) ;
366  return (false) ;
367  }
368 
369  if (n_row < 0) /* n_row must be >= 0 */
370  {
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  {
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  {
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  {
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 */
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>
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 */
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  {
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. */
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 
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 
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>
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>
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>
1446 static inline void order_children
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>
1547 static void detect_super_cols
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 ;
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>
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>
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
static void order_children(IndexType n_col, colamd_col< IndexType > Col[], IndexType p[])
static IndexType clear_mark(IndexType n_row, Colamd_Row< IndexType > Row[])
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
union internal::Colamd_Row::@562 shared1
union internal::colamd_col::@560 shared3
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
static IndexType clear_mark(IndexType n_row, Colamd_Row< IndexType > Row[])
if n return
union internal::colamd_col::@558 shared1
bool stats
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)
Definition: Eigen_Colamd.h:699
static void detect_super_cols(colamd_col< IndexType > Col[], IndexType A[], IndexType head[], IndexType row_start, IndexType row_length)
static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS])
Computes a column ordering using the column approximate minimum degree ordering.
Definition: Eigen_Colamd.h:322
m row(1)
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)
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
IndexType colamd_r(IndexType n_row)
Definition: Eigen_Colamd.h:206
static IndexType garbage_collection(IndexType n_row, IndexType n_col, Colamd_Row< IndexType > Row[], colamd_col< IndexType > Col[], IndexType A[], IndexType *pfree)
EIGEN_DEVICE_FUNC SegmentReturnType head(Index n)
This is the const version of head(Index).
Definition: BlockMethods.h:919
static void colamd_set_defaults(double knobs[COLAMD_KNOBS])
set default parameters The use of this routine is optional.
Definition: Eigen_Colamd.h:286
union internal::colamd_col::@561 shared4
IndexType colamd_c(IndexType n_col)
Definition: Eigen_Colamd.h:202
float * p
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)
static void detect_super_cols(colamd_col< IndexType > Col[], IndexType A[], IndexType head[], IndexType row_start, IndexType row_length)
m col(1)
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])
Definition: Eigen_Colamd.h:483
union internal::colamd_col::@559 shared2
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])
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)
Definition: Eigen_Colamd.h:936
IndexType colamd_recommended(IndexType nnz, IndexType n_row, IndexType n_col)
Returns the recommended value of Alen.
Definition: Eigen_Colamd.h:257
std::ptrdiff_t j
static void order_children(IndexType n_col, colamd_col< IndexType > Col[], IndexType p[])
static IndexType garbage_collection(IndexType n_row, IndexType n_col, Colamd_Row< IndexType > Row[], colamd_col< IndexType > Col[], IndexType A[], IndexType *pfree)


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:00