Home | History | Annotate | Download | only in OrderingMethods
      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 (at) 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 (at) cise.ufl.edu), University of Florida.  The algorithm was
     14 //   developed in collaboration with John Gilbert, Xerox PARC, and Esmond
     15 //   Ng, Oak Ridge National Laboratory.
     16 //
     17 //     Date:
     18 //
     19 //   September 8, 2003.  Version 2.3.
     20 //
     21 //     Acknowledgements:
     22 //
     23 //   This work was supported by the National Science Foundation, under
     24 //   grants DMS-9504974 and DMS-9803599.
     25 //
     26 //     Notice:
     27 //
     28 //   Copyright (c) 1998-2003 by the University of Florida.
     29 //   All Rights Reserved.
     30 //
     31 //   THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
     32 //   EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
     33 //
     34 //   Permission is hereby granted to use, copy, modify, and/or distribute
     35 //   this program, provided that the Copyright, this License, and the
     36 //   Availability of the original version is retained on all copies and made
     37 //   accessible to the end-user of any code or package that includes COLAMD
     38 //   or any modified version of COLAMD.
     39 //
     40 //     Availability:
     41 //
     42 //   The colamd/symamd library is available at
     43 //
     44 //       http://www.suitesparse.com
     45 
     46 
     47 #ifndef EIGEN_COLAMD_H
     48 #define EIGEN_COLAMD_H
     49 
     50 namespace internal {
     51 /* Ensure that debugging is turned off: */
     52 #ifndef COLAMD_NDEBUG
     53 #define COLAMD_NDEBUG
     54 #endif /* NDEBUG */
     55 /* ========================================================================== */
     56 /* === Knob and statistics definitions ====================================== */
     57 /* ========================================================================== */
     58 
     59 /* size of the knobs [ ] array.  Only knobs [0..1] are currently used. */
     60 #define COLAMD_KNOBS 20
     61 
     62 /* number of output statistics.  Only stats [0..6] are currently used. */
     63 #define COLAMD_STATS 20
     64 
     65 /* knobs [0] and stats [0]: dense row knob and output statistic. */
     66 #define COLAMD_DENSE_ROW 0
     67 
     68 /* knobs [1] and stats [1]: dense column knob and output statistic. */
     69 #define COLAMD_DENSE_COL 1
     70 
     71 /* stats [2]: memory defragmentation count output statistic */
     72 #define COLAMD_DEFRAG_COUNT 2
     73 
     74 /* stats [3]: colamd status:  zero OK, > 0 warning or notice, < 0 error */
     75 #define COLAMD_STATUS 3
     76 
     77 /* stats [4..6]: error info, or info on jumbled columns */
     78 #define COLAMD_INFO1 4
     79 #define COLAMD_INFO2 5
     80 #define COLAMD_INFO3 6
     81 
     82 /* error codes returned in stats [3]: */
     83 #define COLAMD_OK       (0)
     84 #define COLAMD_OK_BUT_JUMBLED     (1)
     85 #define COLAMD_ERROR_A_not_present    (-1)
     86 #define COLAMD_ERROR_p_not_present    (-2)
     87 #define COLAMD_ERROR_nrow_negative    (-3)
     88 #define COLAMD_ERROR_ncol_negative    (-4)
     89 #define COLAMD_ERROR_nnz_negative   (-5)
     90 #define COLAMD_ERROR_p0_nonzero     (-6)
     91 #define COLAMD_ERROR_A_too_small    (-7)
     92 #define COLAMD_ERROR_col_length_negative  (-8)
     93 #define COLAMD_ERROR_row_index_out_of_bounds  (-9)
     94 #define COLAMD_ERROR_out_of_memory    (-10)
     95 #define COLAMD_ERROR_internal_error   (-999)
     96 
     97 /* ========================================================================== */
     98 /* === Definitions ========================================================== */
     99 /* ========================================================================== */
    100 
    101 #define ONES_COMPLEMENT(r) (-(r)-1)
    102 
    103 /* -------------------------------------------------------------------------- */
    104 
    105 #define COLAMD_EMPTY (-1)
    106 
    107 /* Row and column status */
    108 #define ALIVE (0)
    109 #define DEAD  (-1)
    110 
    111 /* Column status */
    112 #define DEAD_PRINCIPAL    (-1)
    113 #define DEAD_NON_PRINCIPAL  (-2)
    114 
    115 /* Macros for row and column status update and checking. */
    116 #define ROW_IS_DEAD(r)      ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
    117 #define ROW_IS_MARKED_DEAD(row_mark)  (row_mark < ALIVE)
    118 #define ROW_IS_ALIVE(r)     (Row [r].shared2.mark >= ALIVE)
    119 #define COL_IS_DEAD(c)      (Col [c].start < ALIVE)
    120 #define COL_IS_ALIVE(c)     (Col [c].start >= ALIVE)
    121 #define COL_IS_DEAD_PRINCIPAL(c)  (Col [c].start == DEAD_PRINCIPAL)
    122 #define KILL_ROW(r)     { Row [r].shared2.mark = DEAD ; }
    123 #define KILL_PRINCIPAL_COL(c)   { Col [c].start = DEAD_PRINCIPAL ; }
    124 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
    125 
    126 /* ========================================================================== */
    127 /* === Colamd reporting mechanism =========================================== */
    128 /* ========================================================================== */
    129 
    130 // == Row and Column structures ==
    131 template <typename IndexType>
    132 struct colamd_col
    133 {
    134   IndexType start ;   /* index for A of first row in this column, or DEAD */
    135   /* if column is dead */
    136   IndexType length ;  /* number of rows in this column */
    137   union
    138   {
    139     IndexType thickness ; /* number of original columns represented by this */
    140     /* col, if the column is alive */
    141     IndexType parent ;  /* parent in parent tree super-column structure, if */
    142     /* the column is dead */
    143   } shared1 ;
    144   union
    145   {
    146     IndexType score ; /* the score used to maintain heap, if col is alive */
    147     IndexType order ; /* pivot ordering of this column, if col is dead */
    148   } shared2 ;
    149   union
    150   {
    151     IndexType headhash ;  /* head of a hash bucket, if col is at the head of */
    152     /* a degree list */
    153     IndexType hash ;  /* hash value, if col is not in a degree list */
    154     IndexType prev ;  /* previous column in degree list, if col is in a */
    155     /* degree list (but not at the head of a degree list) */
    156   } shared3 ;
    157   union
    158   {
    159     IndexType degree_next ; /* next column, if col is in a degree list */
    160     IndexType hash_next ;   /* next column, if col is in a hash list */
    161   } shared4 ;
    162 
    163 };
    164 
    165 template <typename IndexType>
    166 struct Colamd_Row
    167 {
    168   IndexType start ;   /* index for A of first col in this row */
    169   IndexType length ;  /* number of principal columns in this row */
    170   union
    171   {
    172     IndexType degree ;  /* number of principal & non-principal columns in row */
    173     IndexType p ;   /* used as a row pointer in init_rows_cols () */
    174   } shared1 ;
    175   union
    176   {
    177     IndexType mark ;  /* for computing set differences and marking dead rows*/
    178     IndexType first_column ;/* first column in row (used in garbage collection) */
    179   } shared2 ;
    180 
    181 };
    182 
    183 /* ========================================================================== */
    184 /* === Colamd recommended memory size ======================================= */
    185 /* ========================================================================== */
    186 
    187 /*
    188   The recommended length Alen of the array A passed to colamd is given by
    189   the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro.  It returns -1 if any
    190   argument is negative.  2*nnz space is required for the row and column
    191   indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is
    192   required for the Col and Row arrays, respectively, which are internal to
    193   colamd.  An additional n_col space is the minimal amount of "elbow room",
    194   and nnz/5 more space is recommended for run time efficiency.
    195 
    196   This macro is not needed when using symamd.
    197 
    198   Explicit typecast to IndexType added Sept. 23, 2002, COLAMD version 2.2, to avoid
    199   gcc -pedantic warning messages.
    200 */
    201 template <typename IndexType>
    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 
    242 /**
    243  * \brief Returns the recommended value of Alen
    244  *
    245  * Returns recommended value of Alen for use by colamd.
    246  * Returns -1 if any input argument is negative.
    247  * The use of this routine or macro is optional.
    248  * Note that the macro uses its arguments   more than once,
    249  * so be careful for side effects, if you pass expressions as arguments to COLAMD_RECOMMENDED.
    250  *
    251  * \param nnz nonzeros in A
    252  * \param n_row number of rows in A
    253  * \param n_col number of columns in A
    254  * \return recommended value of Alen for use by colamd
    255  */
    256 template <typename IndexType>
    257 inline IndexType colamd_recommended ( IndexType nnz, IndexType n_row, IndexType n_col)
    258 {
    259   if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0)
    260     return (-1);
    261   else
    262     return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5));
    263 }
    264 
    265 /**
    266  * \brief set default parameters  The use of this routine is optional.
    267  *
    268  * Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col)
    269  * entries are removed prior to ordering.  Columns with more than
    270  * (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to
    271  * ordering, and placed last in the output column ordering.
    272  *
    273  * COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
    274  * respectively, in colamd.h.  Default values of these two knobs
    275  * are both 0.5.  Currently, only knobs [0] and knobs [1] are
    276  * used, but future versions may use more knobs.  If so, they will
    277  * be properly set to their defaults by the future version of
    278  * colamd_set_defaults, so that the code that calls colamd will
    279  * not need to change, assuming that you either use
    280  * colamd_set_defaults, or pass a (double *) NULL pointer as the
    281  * knobs array to colamd or symamd.
    282  *
    283  * \param knobs parameter settings for colamd
    284  */
    285 
    286 static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS])
    287 {
    288   /* === Local variables ================================================== */
    289 
    290   int i ;
    291 
    292   if (!knobs)
    293   {
    294     return ;      /* no knobs to initialize */
    295   }
    296   for (i = 0 ; i < COLAMD_KNOBS ; i++)
    297   {
    298     knobs [i] = 0 ;
    299   }
    300   knobs [COLAMD_DENSE_ROW] = 0.5 ;  /* ignore rows over 50% dense */
    301   knobs [COLAMD_DENSE_COL] = 0.5 ;  /* ignore columns over 50% dense */
    302 }
    303 
    304 /**
    305  * \brief  Computes a column ordering using the column approximate minimum degree ordering
    306  *
    307  * Computes a column ordering (Q) of A such that P(AQ)=LU or
    308  * (AQ)'AQ=LL' have less fill-in and require fewer floating point
    309  * operations than factorizing the unpermuted matrix A or A'A,
    310  * respectively.
    311  *
    312  *
    313  * \param n_row number of rows in A
    314  * \param n_col number of columns in A
    315  * \param Alen, size of the array A
    316  * \param A row indices of the matrix, of size ALen
    317  * \param p column pointers of A, of size n_col+1
    318  * \param knobs parameter settings for colamd
    319  * \param stats colamd output statistics and error codes
    320  */
    321 template <typename IndexType>
    322 static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS])
    323 {
    324   /* === Local variables ================================================== */
    325 
    326   IndexType i ;     /* loop index */
    327   IndexType nnz ;     /* nonzeros in A */
    328   IndexType Row_size ;    /* size of Row [], in integers */
    329   IndexType Col_size ;    /* size of Col [], in integers */
    330   IndexType need ;      /* minimum required length of A */
    331   Colamd_Row<IndexType> *Row ;   /* pointer into A of Row [0..n_row] array */
    332   colamd_col<IndexType> *Col ;   /* pointer into A of Col [0..n_col] array */
    333   IndexType n_col2 ;    /* number of non-dense, non-empty columns */
    334   IndexType n_row2 ;    /* number of non-dense, non-empty rows */
    335   IndexType ngarbage ;    /* number of garbage collections performed */
    336   IndexType max_deg ;   /* maximum row degree */
    337   double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
    338 
    339 
    340   /* === Check the input arguments ======================================== */
    341 
    342   if (!stats)
    343   {
    344     COLAMD_DEBUG0 (("colamd: stats not present\n")) ;
    345     return (false) ;
    346   }
    347   for (i = 0 ; i < COLAMD_STATS ; i++)
    348   {
    349     stats [i] = 0 ;
    350   }
    351   stats [COLAMD_STATUS] = COLAMD_OK ;
    352   stats [COLAMD_INFO1] = -1 ;
    353   stats [COLAMD_INFO2] = -1 ;
    354 
    355   if (!A)   /* A is not present */
    356   {
    357     stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
    358     COLAMD_DEBUG0 (("colamd: A not present\n")) ;
    359     return (false) ;
    360   }
    361 
    362   if (!p)   /* p is not present */
    363   {
    364     stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
    365     COLAMD_DEBUG0 (("colamd: p not present\n")) ;
    366     return (false) ;
    367   }
    368 
    369   if (n_row < 0)  /* n_row must be >= 0 */
    370   {
    371     stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
    372     stats [COLAMD_INFO1] = n_row ;
    373     COLAMD_DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
    374     return (false) ;
    375   }
    376 
    377   if (n_col < 0)  /* n_col must be >= 0 */
    378   {
    379     stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
    380     stats [COLAMD_INFO1] = n_col ;
    381     COLAMD_DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
    382     return (false) ;
    383   }
    384 
    385   nnz = p [n_col] ;
    386   if (nnz < 0)  /* nnz must be >= 0 */
    387   {
    388     stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
    389     stats [COLAMD_INFO1] = nnz ;
    390     COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
    391     return (false) ;
    392   }
    393 
    394   if (p [0] != 0)
    395   {
    396     stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
    397     stats [COLAMD_INFO1] = p [0] ;
    398     COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
    399     return (false) ;
    400   }
    401 
    402   /* === If no knobs, set default knobs =================================== */
    403 
    404   if (!knobs)
    405   {
    406     colamd_set_defaults (default_knobs) ;
    407     knobs = default_knobs ;
    408   }
    409 
    410   /* === Allocate the Row and Col arrays from array A ===================== */
    411 
    412   Col_size = colamd_c (n_col) ;
    413   Row_size = colamd_r (n_row) ;
    414   need = 2*nnz + n_col + Col_size + Row_size ;
    415 
    416   if (need > Alen)
    417   {
    418     /* not enough space in array A to perform the ordering */
    419     stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
    420     stats [COLAMD_INFO1] = need ;
    421     stats [COLAMD_INFO2] = Alen ;
    422     COLAMD_DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
    423     return (false) ;
    424   }
    425 
    426   Alen -= Col_size + Row_size ;
    427   Col = (colamd_col<IndexType> *) &A [Alen] ;
    428   Row = (Colamd_Row<IndexType> *) &A [Alen + Col_size] ;
    429 
    430   /* === Construct the row and column data structures ===================== */
    431 
    432   if (!Eigen::internal::init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
    433   {
    434     /* input matrix is invalid */
    435     COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ;
    436     return (false) ;
    437   }
    438 
    439   /* === Initialize scores, kill dense rows/columns ======================= */
    440 
    441   Eigen::internal::init_scoring (n_row, n_col, Row, Col, A, p, knobs,
    442 		&n_row2, &n_col2, &max_deg) ;
    443 
    444   /* === Order the supercolumns =========================================== */
    445 
    446   ngarbage = Eigen::internal::find_ordering (n_row, n_col, Alen, Row, Col, A, p,
    447 			    n_col2, max_deg, 2*nnz) ;
    448 
    449   /* === Order the non-principal columns ================================== */
    450 
    451   Eigen::internal::order_children (n_col, Col, p) ;
    452 
    453   /* === Return statistics in stats ======================================= */
    454 
    455   stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
    456   stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
    457   stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
    458   COLAMD_DEBUG0 (("colamd: done.\n")) ;
    459   return (true) ;
    460 }
    461 
    462 /* ========================================================================== */
    463 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
    464 /* ========================================================================== */
    465 
    466 /* There are no user-callable routines beyond this point in the file */
    467 
    468 
    469 /* ========================================================================== */
    470 /* === init_rows_cols ======================================================= */
    471 /* ========================================================================== */
    472 
    473 /*
    474   Takes the column form of the matrix in A and creates the row form of the
    475   matrix.  Also, row and column attributes are stored in the Col and Row
    476   structs.  If the columns are un-sorted or contain duplicate row indices,
    477   this routine will also sort and remove duplicate row indices from the
    478   column form of the matrix.  Returns false if the matrix is invalid,
    479   true otherwise.  Not user-callable.
    480 */
    481 template <typename IndexType>
    482 static IndexType init_rows_cols  /* returns true if OK, or false otherwise */
    483   (
    484     /* === Parameters ======================================================= */
    485 
    486     IndexType n_row,      /* number of rows of A */
    487     IndexType n_col,      /* number of columns of A */
    488     Colamd_Row<IndexType> Row [],    /* of size n_row+1 */
    489     colamd_col<IndexType> Col [],    /* of size n_col+1 */
    490     IndexType A [],     /* row indices of A, of size Alen */
    491     IndexType p [],     /* pointers to columns in A, of size n_col+1 */
    492     IndexType stats [COLAMD_STATS]  /* colamd statistics */
    493     )
    494 {
    495   /* === Local variables ================================================== */
    496 
    497   IndexType col ;     /* a column index */
    498   IndexType row ;     /* a row index */
    499   IndexType *cp ;     /* a column pointer */
    500   IndexType *cp_end ;   /* a pointer to the end of a column */
    501   IndexType *rp ;     /* a row pointer */
    502   IndexType *rp_end ;   /* a pointer to the end of a row */
    503   IndexType last_row ;    /* previous row */
    504 
    505   /* === Initialize columns, and check column pointers ==================== */
    506 
    507   for (col = 0 ; col < n_col ; col++)
    508   {
    509     Col [col].start = p [col] ;
    510     Col [col].length = p [col+1] - p [col] ;
    511 
    512     if ((Col [col].length) < 0) // extra parentheses to work-around gcc bug 10200
    513     {
    514       /* column pointers must be non-decreasing */
    515       stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
    516       stats [COLAMD_INFO1] = col ;
    517       stats [COLAMD_INFO2] = Col [col].length ;
    518       COLAMD_DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
    519       return (false) ;
    520     }
    521 
    522     Col [col].shared1.thickness = 1 ;
    523     Col [col].shared2.score = 0 ;
    524     Col [col].shared3.prev = COLAMD_EMPTY ;
    525     Col [col].shared4.degree_next = COLAMD_EMPTY ;
    526   }
    527 
    528   /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
    529 
    530   /* === Scan columns, compute row degrees, and check row indices ========= */
    531 
    532   stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
    533 
    534   for (row = 0 ; row < n_row ; row++)
    535   {
    536     Row [row].length = 0 ;
    537     Row [row].shared2.mark = -1 ;
    538   }
    539 
    540   for (col = 0 ; col < n_col ; col++)
    541   {
    542     last_row = -1 ;
    543 
    544     cp = &A [p [col]] ;
    545     cp_end = &A [p [col+1]] ;
    546 
    547     while (cp < cp_end)
    548     {
    549       row = *cp++ ;
    550 
    551       /* make sure row indices within range */
    552       if (row < 0 || row >= n_row)
    553       {
    554 	stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
    555 	stats [COLAMD_INFO1] = col ;
    556 	stats [COLAMD_INFO2] = row ;
    557 	stats [COLAMD_INFO3] = n_row ;
    558 	COLAMD_DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
    559 	return (false) ;
    560       }
    561 
    562       if (row <= last_row || Row [row].shared2.mark == col)
    563       {
    564 	/* row index are unsorted or repeated (or both), thus col */
    565 	/* is jumbled.  This is a notice, not an error condition. */
    566 	stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
    567 	stats [COLAMD_INFO1] = col ;
    568 	stats [COLAMD_INFO2] = row ;
    569 	(stats [COLAMD_INFO3]) ++ ;
    570 	COLAMD_DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
    571       }
    572 
    573       if (Row [row].shared2.mark != col)
    574       {
    575 	Row [row].length++ ;
    576       }
    577       else
    578       {
    579 	/* this is a repeated entry in the column, */
    580 	/* it will be removed */
    581 	Col [col].length-- ;
    582       }
    583 
    584       /* mark the row as having been seen in this column */
    585       Row [row].shared2.mark = col ;
    586 
    587       last_row = row ;
    588     }
    589   }
    590 
    591   /* === Compute row pointers ============================================= */
    592 
    593   /* row form of the matrix starts directly after the column */
    594   /* form of matrix in A */
    595   Row [0].start = p [n_col] ;
    596   Row [0].shared1.p = Row [0].start ;
    597   Row [0].shared2.mark = -1 ;
    598   for (row = 1 ; row < n_row ; row++)
    599   {
    600     Row [row].start = Row [row-1].start + Row [row-1].length ;
    601     Row [row].shared1.p = Row [row].start ;
    602     Row [row].shared2.mark = -1 ;
    603   }
    604 
    605   /* === Create row form ================================================== */
    606 
    607   if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
    608   {
    609     /* if cols jumbled, watch for repeated row indices */
    610     for (col = 0 ; col < n_col ; col++)
    611     {
    612       cp = &A [p [col]] ;
    613       cp_end = &A [p [col+1]] ;
    614       while (cp < cp_end)
    615       {
    616 	row = *cp++ ;
    617 	if (Row [row].shared2.mark != col)
    618 	{
    619 	  A [(Row [row].shared1.p)++] = col ;
    620 	  Row [row].shared2.mark = col ;
    621 	}
    622       }
    623     }
    624   }
    625   else
    626   {
    627     /* if cols not jumbled, we don't need the mark (this is faster) */
    628     for (col = 0 ; col < n_col ; col++)
    629     {
    630       cp = &A [p [col]] ;
    631       cp_end = &A [p [col+1]] ;
    632       while (cp < cp_end)
    633       {
    634 	A [(Row [*cp++].shared1.p)++] = col ;
    635       }
    636     }
    637   }
    638 
    639   /* === Clear the row marks and set row degrees ========================== */
    640 
    641   for (row = 0 ; row < n_row ; row++)
    642   {
    643     Row [row].shared2.mark = 0 ;
    644     Row [row].shared1.degree = Row [row].length ;
    645   }
    646 
    647   /* === See if we need to re-create columns ============================== */
    648 
    649   if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
    650   {
    651     COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
    652 
    653 
    654     /* === Compute col pointers ========================================= */
    655 
    656     /* col form of the matrix starts at A [0]. */
    657     /* Note, we may have a gap between the col form and the row */
    658     /* form if there were duplicate entries, if so, it will be */
    659     /* removed upon the first garbage collection */
    660     Col [0].start = 0 ;
    661     p [0] = Col [0].start ;
    662     for (col = 1 ; col < n_col ; col++)
    663     {
    664       /* note that the lengths here are for pruned columns, i.e. */
    665       /* no duplicate row indices will exist for these columns */
    666       Col [col].start = Col [col-1].start + Col [col-1].length ;
    667       p [col] = Col [col].start ;
    668     }
    669 
    670     /* === Re-create col form =========================================== */
    671 
    672     for (row = 0 ; row < n_row ; row++)
    673     {
    674       rp = &A [Row [row].start] ;
    675       rp_end = rp + Row [row].length ;
    676       while (rp < rp_end)
    677       {
    678 	A [(p [*rp++])++] = row ;
    679       }
    680     }
    681   }
    682 
    683   /* === Done.  Matrix is not (or no longer) jumbled ====================== */
    684 
    685   return (true) ;
    686 }
    687 
    688 
    689 /* ========================================================================== */
    690 /* === init_scoring ========================================================= */
    691 /* ========================================================================== */
    692 
    693 /*
    694   Kills dense or empty columns and rows, calculates an initial score for
    695   each column, and places all columns in the degree lists.  Not user-callable.
    696 */
    697 template <typename IndexType>
    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 (head [min_score] == COLAMD_EMPTY && min_score < n_col)
   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
   1447 (
   1448   /* === Parameters ======================================================= */
   1449 
   1450   IndexType n_col,      /* number of columns of A */
   1451   colamd_col<IndexType> Col [],    /* of size n_col+1 */
   1452   IndexType p []      /* p [0 ... n_col-1] is the column permutation*/
   1453   )
   1454 {
   1455   /* === Local variables ================================================== */
   1456 
   1457   IndexType i ;     /* loop counter for all columns */
   1458   IndexType c ;     /* column index */
   1459   IndexType parent ;    /* index of column's parent */
   1460   IndexType order ;     /* column's order */
   1461 
   1462   /* === Order each non-principal column ================================== */
   1463 
   1464   for (i = 0 ; i < n_col ; i++)
   1465   {
   1466     /* find an un-ordered non-principal column */
   1467     COLAMD_ASSERT (COL_IS_DEAD (i)) ;
   1468     if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY)
   1469     {
   1470       parent = i ;
   1471       /* once found, find its principal parent */
   1472       do
   1473       {
   1474 	parent = Col [parent].shared1.parent ;
   1475       } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
   1476 
   1477       /* now, order all un-ordered non-principal columns along path */
   1478       /* to this parent.  collapse tree at the same time */
   1479       c = i ;
   1480       /* get order of parent */
   1481       order = Col [parent].shared2.order ;
   1482 
   1483       do
   1484       {
   1485 	COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ;
   1486 
   1487 	/* order this column */
   1488 	Col [c].shared2.order = order++ ;
   1489 	/* collaps tree */
   1490 	Col [c].shared1.parent = parent ;
   1491 
   1492 	/* get immediate parent of this column */
   1493 	c = Col [c].shared1.parent ;
   1494 
   1495 	/* continue until we hit an ordered column.  There are */
   1496 	/* guarranteed not to be anymore unordered columns */
   1497 	/* above an ordered column */
   1498       } while (Col [c].shared2.order == COLAMD_EMPTY) ;
   1499 
   1500       /* re-order the super_col parent to largest order for this group */
   1501       Col [parent].shared2.order = order ;
   1502     }
   1503   }
   1504 
   1505   /* === Generate the permutation ========================================= */
   1506 
   1507   for (c = 0 ; c < n_col ; c++)
   1508   {
   1509     p [Col [c].shared2.order] = c ;
   1510   }
   1511 }
   1512 
   1513 
   1514 /* ========================================================================== */
   1515 /* === detect_super_cols ==================================================== */
   1516 /* ========================================================================== */
   1517 
   1518 /*
   1519   Detects supercolumns by finding matches between columns in the hash buckets.
   1520   Check amongst columns in the set A [row_start ... row_start + row_length-1].
   1521   The columns under consideration are currently *not* in the degree lists,
   1522   and have already been placed in the hash buckets.
   1523 
   1524   The hash bucket for columns whose hash function is equal to h is stored
   1525   as follows:
   1526 
   1527   if head [h] is >= 0, then head [h] contains a degree list, so:
   1528 
   1529   head [h] is the first column in degree bucket h.
   1530   Col [head [h]].headhash gives the first column in hash bucket h.
   1531 
   1532   otherwise, the degree list is empty, and:
   1533 
   1534   -(head [h] + 2) is the first column in hash bucket h.
   1535 
   1536   For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
   1537   column" pointer.  Col [c].shared3.hash is used instead as the hash number
   1538   for that column.  The value of Col [c].shared4.hash_next is the next column
   1539   in the same hash bucket.
   1540 
   1541   Assuming no, or "few" hash collisions, the time taken by this routine is
   1542   linear in the sum of the sizes (lengths) of each column whose score has
   1543   just been computed in the approximate degree computation.
   1544   Not user-callable.
   1545 */
   1546 template <typename IndexType>
   1547 static void detect_super_cols
   1548 (
   1549   /* === Parameters ======================================================= */
   1550 
   1551   colamd_col<IndexType> Col [],    /* of size n_col+1 */
   1552   IndexType A [],     /* row indices of A */
   1553   IndexType head [],    /* head of degree lists and hash buckets */
   1554   IndexType row_start,    /* pointer to set of columns to check */
   1555   IndexType row_length    /* number of columns to check */
   1556 )
   1557 {
   1558   /* === Local variables ================================================== */
   1559 
   1560   IndexType hash ;      /* hash value for a column */
   1561   IndexType *rp ;     /* pointer to a row */
   1562   IndexType c ;     /* a column index */
   1563   IndexType super_c ;   /* column index of the column to absorb into */
   1564   IndexType *cp1 ;      /* column pointer for column super_c */
   1565   IndexType *cp2 ;      /* column pointer for column c */
   1566   IndexType length ;    /* length of column super_c */
   1567   IndexType prev_c ;    /* column preceding c in hash bucket */
   1568   IndexType i ;     /* loop counter */
   1569   IndexType *rp_end ;   /* pointer to the end of the row */
   1570   IndexType col ;     /* a column index in the row to check */
   1571   IndexType head_column ;   /* first column in hash bucket or degree list */
   1572   IndexType first_col ;   /* first column in hash bucket */
   1573 
   1574   /* === Consider each column in the row ================================== */
   1575 
   1576   rp = &A [row_start] ;
   1577   rp_end = rp + row_length ;
   1578   while (rp < rp_end)
   1579   {
   1580     col = *rp++ ;
   1581     if (COL_IS_DEAD (col))
   1582     {
   1583       continue ;
   1584     }
   1585 
   1586     /* get hash number for this column */
   1587     hash = Col [col].shared3.hash ;
   1588     COLAMD_ASSERT (hash <= n_col) ;
   1589 
   1590     /* === Get the first column in this hash bucket ===================== */
   1591 
   1592     head_column = head [hash] ;
   1593     if (head_column > COLAMD_EMPTY)
   1594     {
   1595       first_col = Col [head_column].shared3.headhash ;
   1596     }
   1597     else
   1598     {
   1599       first_col = - (head_column + 2) ;
   1600     }
   1601 
   1602     /* === Consider each column in the hash bucket ====================== */
   1603 
   1604     for (super_c = first_col ; super_c != COLAMD_EMPTY ;
   1605 	 super_c = Col [super_c].shared4.hash_next)
   1606     {
   1607       COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ;
   1608       COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ;
   1609       length = Col [super_c].length ;
   1610 
   1611       /* prev_c is the column preceding column c in the hash bucket */
   1612       prev_c = super_c ;
   1613 
   1614       /* === Compare super_c with all columns after it ================ */
   1615 
   1616       for (c = Col [super_c].shared4.hash_next ;
   1617 	   c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next)
   1618       {
   1619 	COLAMD_ASSERT (c != super_c) ;
   1620 	COLAMD_ASSERT (COL_IS_ALIVE (c)) ;
   1621 	COLAMD_ASSERT (Col [c].shared3.hash == hash) ;
   1622 
   1623 	/* not identical if lengths or scores are different */
   1624 	if (Col [c].length != length ||
   1625 	    Col [c].shared2.score != Col [super_c].shared2.score)
   1626 	{
   1627 	  prev_c = c ;
   1628 	  continue ;
   1629 	}
   1630 
   1631 	/* compare the two columns */
   1632 	cp1 = &A [Col [super_c].start] ;
   1633 	cp2 = &A [Col [c].start] ;
   1634 
   1635 	for (i = 0 ; i < length ; i++)
   1636 	{
   1637 	  /* the columns are "clean" (no dead rows) */
   1638 	  COLAMD_ASSERT (ROW_IS_ALIVE (*cp1))  ;
   1639 	  COLAMD_ASSERT (ROW_IS_ALIVE (*cp2))  ;
   1640 	  /* row indices will same order for both supercols, */
   1641 	  /* no gather scatter nessasary */
   1642 	  if (*cp1++ != *cp2++)
   1643 	  {
   1644 	    break ;
   1645 	  }
   1646 	}
   1647 
   1648 	/* the two columns are different if the for-loop "broke" */
   1649 	if (i != length)
   1650 	{
   1651 	  prev_c = c ;
   1652 	  continue ;
   1653 	}
   1654 
   1655 	/* === Got it!  two columns are identical =================== */
   1656 
   1657 	COLAMD_ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
   1658 
   1659 	Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
   1660 	Col [c].shared1.parent = super_c ;
   1661 	KILL_NON_PRINCIPAL_COL (c) ;
   1662 	/* order c later, in order_children() */
   1663 	Col [c].shared2.order = COLAMD_EMPTY ;
   1664 	/* remove c from hash bucket */
   1665 	Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
   1666       }
   1667     }
   1668 
   1669     /* === Empty this hash bucket ======================================= */
   1670 
   1671     if (head_column > COLAMD_EMPTY)
   1672     {
   1673       /* corresponding degree list "hash" is not empty */
   1674       Col [head_column].shared3.headhash = COLAMD_EMPTY ;
   1675     }
   1676     else
   1677     {
   1678       /* corresponding degree list "hash" is empty */
   1679       head [hash] = COLAMD_EMPTY ;
   1680     }
   1681   }
   1682 }
   1683 
   1684 
   1685 /* ========================================================================== */
   1686 /* === garbage_collection =================================================== */
   1687 /* ========================================================================== */
   1688 
   1689 /*
   1690   Defragments and compacts columns and rows in the workspace A.  Used when
   1691   all avaliable memory has been used while performing row merging.  Returns
   1692   the index of the first free position in A, after garbage collection.  The
   1693   time taken by this routine is linear is the size of the array A, which is
   1694   itself linear in the number of nonzeros in the input matrix.
   1695   Not user-callable.
   1696 */
   1697 template <typename IndexType>
   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
   1844