Eigen_Colamd.h
Go to the documentation of this file.
00001 // // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 // This file is modified from the colamd/symamd library. The copyright is below
00011 
00012 //   The authors of the code itself are Stefan I. Larimore and Timothy A.
00013 //   Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was
00014 //   developed in collaboration with John Gilbert, Xerox PARC, and Esmond
00015 //   Ng, Oak Ridge National Laboratory.
00016 // 
00017 //     Date:
00018 // 
00019 //   September 8, 2003.  Version 2.3.
00020 // 
00021 //     Acknowledgements:
00022 // 
00023 //   This work was supported by the National Science Foundation, under
00024 //   grants DMS-9504974 and DMS-9803599.
00025 // 
00026 //     Notice:
00027 // 
00028 //   Copyright (c) 1998-2003 by the University of Florida.
00029 //   All Rights Reserved.
00030 // 
00031 //   THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
00032 //   EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
00033 // 
00034 //   Permission is hereby granted to use, copy, modify, and/or distribute
00035 //   this program, provided that the Copyright, this License, and the
00036 //   Availability of the original version is retained on all copies and made
00037 //   accessible to the end-user of any code or package that includes COLAMD
00038 //   or any modified version of COLAMD. 
00039 // 
00040 //     Availability:
00041 // 
00042 //   The colamd/symamd library is available at
00043 // 
00044 //       http://www.cise.ufl.edu/research/sparse/colamd/
00045 
00046 //   This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h
00047 //   file.  It is required by the colamd.c, colamdmex.c, and symamdmex.c
00048 //   files, and by any C code that calls the routines whose prototypes are
00049 //   listed below, or that uses the colamd/symamd definitions listed below.
00050   
00051 #ifndef EIGEN_COLAMD_H
00052 #define EIGEN_COLAMD_H
00053 
00054 namespace internal {
00055 /* Ensure that debugging is turned off: */
00056 #ifndef COLAMD_NDEBUG
00057 #define COLAMD_NDEBUG
00058 #endif /* NDEBUG */
00059 /* ========================================================================== */
00060 /* === Knob and statistics definitions ====================================== */
00061 /* ========================================================================== */
00062 
00063 /* size of the knobs [ ] array.  Only knobs [0..1] are currently used. */
00064 #define COLAMD_KNOBS 20
00065 
00066 /* number of output statistics.  Only stats [0..6] are currently used. */
00067 #define COLAMD_STATS 20 
00068 
00069 /* knobs [0] and stats [0]: dense row knob and output statistic. */
00070 #define COLAMD_DENSE_ROW 0
00071 
00072 /* knobs [1] and stats [1]: dense column knob and output statistic. */
00073 #define COLAMD_DENSE_COL 1
00074 
00075 /* stats [2]: memory defragmentation count output statistic */
00076 #define COLAMD_DEFRAG_COUNT 2
00077 
00078 /* stats [3]: colamd status:  zero OK, > 0 warning or notice, < 0 error */
00079 #define COLAMD_STATUS 3
00080 
00081 /* stats [4..6]: error info, or info on jumbled columns */ 
00082 #define COLAMD_INFO1 4
00083 #define COLAMD_INFO2 5
00084 #define COLAMD_INFO3 6
00085 
00086 /* error codes returned in stats [3]: */
00087 #define COLAMD_OK       (0)
00088 #define COLAMD_OK_BUT_JUMBLED     (1)
00089 #define COLAMD_ERROR_A_not_present    (-1)
00090 #define COLAMD_ERROR_p_not_present    (-2)
00091 #define COLAMD_ERROR_nrow_negative    (-3)
00092 #define COLAMD_ERROR_ncol_negative    (-4)
00093 #define COLAMD_ERROR_nnz_negative   (-5)
00094 #define COLAMD_ERROR_p0_nonzero     (-6)
00095 #define COLAMD_ERROR_A_too_small    (-7)
00096 #define COLAMD_ERROR_col_length_negative  (-8)
00097 #define COLAMD_ERROR_row_index_out_of_bounds  (-9)
00098 #define COLAMD_ERROR_out_of_memory    (-10)
00099 #define COLAMD_ERROR_internal_error   (-999)
00100 
00101 /* ========================================================================== */
00102 /* === Definitions ========================================================== */
00103 /* ========================================================================== */
00104 
00105 #define COLAMD_MAX(a,b) (((a) > (b)) ? (a) : (b))
00106 #define COLAMD_MIN(a,b) (((a) < (b)) ? (a) : (b))
00107 
00108 #define ONES_COMPLEMENT(r) (-(r)-1)
00109 
00110 /* -------------------------------------------------------------------------- */
00111 
00112 #define COLAMD_EMPTY (-1)
00113 
00114 /* Row and column status */
00115 #define ALIVE (0)
00116 #define DEAD  (-1)
00117 
00118 /* Column status */
00119 #define DEAD_PRINCIPAL    (-1)
00120 #define DEAD_NON_PRINCIPAL  (-2)
00121 
00122 /* Macros for row and column status update and checking. */
00123 #define ROW_IS_DEAD(r)      ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
00124 #define ROW_IS_MARKED_DEAD(row_mark)  (row_mark < ALIVE)
00125 #define ROW_IS_ALIVE(r)     (Row [r].shared2.mark >= ALIVE)
00126 #define COL_IS_DEAD(c)      (Col [c].start < ALIVE)
00127 #define COL_IS_ALIVE(c)     (Col [c].start >= ALIVE)
00128 #define COL_IS_DEAD_PRINCIPAL(c)  (Col [c].start == DEAD_PRINCIPAL)
00129 #define KILL_ROW(r)     { Row [r].shared2.mark = DEAD ; }
00130 #define KILL_PRINCIPAL_COL(c)   { Col [c].start = DEAD_PRINCIPAL ; }
00131 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
00132 
00133 /* ========================================================================== */
00134 /* === Colamd reporting mechanism =========================================== */
00135 /* ========================================================================== */
00136 
00137 // == Row and Column structures ==
00138 template <typename Index>
00139 struct colamd_col
00140 {
00141   Index start ;   /* index for A of first row in this column, or DEAD */
00142   /* if column is dead */
00143   Index length ;  /* number of rows in this column */
00144   union
00145   {
00146     Index thickness ; /* number of original columns represented by this */
00147     /* col, if the column is alive */
00148     Index parent ;  /* parent in parent tree super-column structure, if */
00149     /* the column is dead */
00150   } shared1 ;
00151   union
00152   {
00153     Index score ; /* the score used to maintain heap, if col is alive */
00154     Index order ; /* pivot ordering of this column, if col is dead */
00155   } shared2 ;
00156   union
00157   {
00158     Index headhash ;  /* head of a hash bucket, if col is at the head of */
00159     /* a degree list */
00160     Index hash ;  /* hash value, if col is not in a degree list */
00161     Index prev ;  /* previous column in degree list, if col is in a */
00162     /* degree list (but not at the head of a degree list) */
00163   } shared3 ;
00164   union
00165   {
00166     Index degree_next ; /* next column, if col is in a degree list */
00167     Index hash_next ;   /* next column, if col is in a hash list */
00168   } shared4 ;
00169   
00170 };
00171  
00172 template <typename Index>
00173 struct Colamd_Row
00174 {
00175   Index start ;   /* index for A of first col in this row */
00176   Index length ;  /* number of principal columns in this row */
00177   union
00178   {
00179     Index degree ;  /* number of principal & non-principal columns in row */
00180     Index p ;   /* used as a row pointer in init_rows_cols () */
00181   } shared1 ;
00182   union
00183   {
00184     Index mark ;  /* for computing set differences and marking dead rows*/
00185     Index first_column ;/* first column in row (used in garbage collection) */
00186   } shared2 ;
00187   
00188 };
00189  
00190 /* ========================================================================== */
00191 /* === Colamd recommended memory size ======================================= */
00192 /* ========================================================================== */
00193  
00194 /*
00195   The recommended length Alen of the array A passed to colamd is given by
00196   the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro.  It returns -1 if any
00197   argument is negative.  2*nnz space is required for the row and column
00198   indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is
00199   required for the Col and Row arrays, respectively, which are internal to
00200   colamd.  An additional n_col space is the minimal amount of "elbow room",
00201   and nnz/5 more space is recommended for run time efficiency.
00202   
00203   This macro is not needed when using symamd.
00204   
00205   Explicit typecast to Index added Sept. 23, 2002, COLAMD version 2.2, to avoid
00206   gcc -pedantic warning messages.
00207 */
00208 template <typename Index>
00209 inline Index colamd_c(Index n_col) 
00210 { return Index( ((n_col) + 1) * sizeof (colamd_col<Index>) / sizeof (Index) ) ; }
00211 
00212 template <typename Index>
00213 inline Index  colamd_r(Index n_row)
00214 { return Index(((n_row) + 1) * sizeof (Colamd_Row<Index>) / sizeof (Index)); }
00215 
00216 // Prototypes of non-user callable routines
00217 template <typename Index>
00218 static Index init_rows_cols (Index n_row, Index n_col, Colamd_Row<Index> Row [], colamd_col<Index> col [], Index A [], Index p [], Index stats[COLAMD_STATS] ); 
00219 
00220 template <typename Index>
00221 static void init_scoring (Index n_row, Index n_col, Colamd_Row<Index> Row [], colamd_col<Index> Col [], Index A [], Index head [], double knobs[COLAMD_KNOBS], Index *p_n_row2, Index *p_n_col2, Index *p_max_deg);
00222 
00223 template <typename Index>
00224 static Index find_ordering (Index n_row, Index n_col, Index Alen, Colamd_Row<Index> Row [], colamd_col<Index> Col [], Index A [], Index head [], Index n_col2, Index max_deg, Index pfree);
00225 
00226 template <typename Index>
00227 static void order_children (Index n_col, colamd_col<Index> Col [], Index p []);
00228 
00229 template <typename Index>
00230 static void detect_super_cols (colamd_col<Index> Col [], Index A [], Index head [], Index row_start, Index row_length ) ;
00231 
00232 template <typename Index>
00233 static Index garbage_collection (Index n_row, Index n_col, Colamd_Row<Index> Row [], colamd_col<Index> Col [], Index A [], Index *pfree) ;
00234 
00235 template <typename Index>
00236 static inline  Index clear_mark (Index n_row, Colamd_Row<Index> Row [] ) ;
00237 
00238 /* === No debugging ========================================================= */
00239 
00240 #define COLAMD_DEBUG0(params) ;
00241 #define COLAMD_DEBUG1(params) ;
00242 #define COLAMD_DEBUG2(params) ;
00243 #define COLAMD_DEBUG3(params) ;
00244 #define COLAMD_DEBUG4(params) ;
00245 
00246 #define COLAMD_ASSERT(expression) ((void) 0)
00247 
00248 
00263 template <typename Index>
00264 inline Index colamd_recommended ( Index nnz, Index n_row, Index n_col)
00265 {
00266   if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0)
00267     return (-1);
00268   else
00269     return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5)); 
00270 }
00271 
00293 static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS])
00294 {
00295   /* === Local variables ================================================== */
00296   
00297   int i ;
00298 
00299   if (!knobs)
00300   {
00301     return ;      /* no knobs to initialize */
00302   }
00303   for (i = 0 ; i < COLAMD_KNOBS ; i++)
00304   {
00305     knobs [i] = 0 ;
00306   }
00307   knobs [COLAMD_DENSE_ROW] = 0.5 ;  /* ignore rows over 50% dense */
00308   knobs [COLAMD_DENSE_COL] = 0.5 ;  /* ignore columns over 50% dense */
00309 }
00310 
00328 template <typename Index>
00329 static bool colamd(Index n_row, Index n_col, Index Alen, Index *A, Index *p, double knobs[COLAMD_KNOBS], Index stats[COLAMD_STATS])
00330 {
00331   /* === Local variables ================================================== */
00332   
00333   Index i ;     /* loop index */
00334   Index nnz ;     /* nonzeros in A */
00335   Index Row_size ;    /* size of Row [], in integers */
00336   Index Col_size ;    /* size of Col [], in integers */
00337   Index need ;      /* minimum required length of A */
00338   Colamd_Row<Index> *Row ;   /* pointer into A of Row [0..n_row] array */
00339   colamd_col<Index> *Col ;   /* pointer into A of Col [0..n_col] array */
00340   Index n_col2 ;    /* number of non-dense, non-empty columns */
00341   Index n_row2 ;    /* number of non-dense, non-empty rows */
00342   Index ngarbage ;    /* number of garbage collections performed */
00343   Index max_deg ;   /* maximum row degree */
00344   double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
00345   
00346   
00347   /* === Check the input arguments ======================================== */
00348   
00349   if (!stats)
00350   {
00351     COLAMD_DEBUG0 (("colamd: stats not present\n")) ;
00352     return (false) ;
00353   }
00354   for (i = 0 ; i < COLAMD_STATS ; i++)
00355   {
00356     stats [i] = 0 ;
00357   }
00358   stats [COLAMD_STATUS] = COLAMD_OK ;
00359   stats [COLAMD_INFO1] = -1 ;
00360   stats [COLAMD_INFO2] = -1 ;
00361   
00362   if (!A)   /* A is not present */
00363   {
00364     stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
00365     COLAMD_DEBUG0 (("colamd: A not present\n")) ;
00366     return (false) ;
00367   }
00368   
00369   if (!p)   /* p is not present */
00370   {
00371     stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
00372     COLAMD_DEBUG0 (("colamd: p not present\n")) ;
00373     return (false) ;
00374   }
00375   
00376   if (n_row < 0)  /* n_row must be >= 0 */
00377   {
00378     stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
00379     stats [COLAMD_INFO1] = n_row ;
00380     COLAMD_DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
00381     return (false) ;
00382   }
00383   
00384   if (n_col < 0)  /* n_col must be >= 0 */
00385   {
00386     stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
00387     stats [COLAMD_INFO1] = n_col ;
00388     COLAMD_DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
00389     return (false) ;
00390   }
00391   
00392   nnz = p [n_col] ;
00393   if (nnz < 0)  /* nnz must be >= 0 */
00394   {
00395     stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
00396     stats [COLAMD_INFO1] = nnz ;
00397     COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
00398     return (false) ;
00399   }
00400   
00401   if (p [0] != 0)
00402   {
00403     stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
00404     stats [COLAMD_INFO1] = p [0] ;
00405     COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
00406     return (false) ;
00407   }
00408   
00409   /* === If no knobs, set default knobs =================================== */
00410   
00411   if (!knobs)
00412   {
00413     colamd_set_defaults (default_knobs) ;
00414     knobs = default_knobs ;
00415   }
00416   
00417   /* === Allocate the Row and Col arrays from array A ===================== */
00418   
00419   Col_size = colamd_c (n_col) ;
00420   Row_size = colamd_r (n_row) ;
00421   need = 2*nnz + n_col + Col_size + Row_size ;
00422   
00423   if (need > Alen)
00424   {
00425     /* not enough space in array A to perform the ordering */
00426     stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
00427     stats [COLAMD_INFO1] = need ;
00428     stats [COLAMD_INFO2] = Alen ;
00429     COLAMD_DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
00430     return (false) ;
00431   }
00432   
00433   Alen -= Col_size + Row_size ;
00434   Col = (colamd_col<Index> *) &A [Alen] ;
00435   Row = (Colamd_Row<Index> *) &A [Alen + Col_size] ;
00436 
00437   /* === Construct the row and column data structures ===================== */
00438   
00439   if (!Eigen::internal::init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
00440   {
00441     /* input matrix is invalid */
00442     COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ;
00443     return (false) ;
00444   }
00445   
00446   /* === Initialize scores, kill dense rows/columns ======================= */
00447 
00448   Eigen::internal::init_scoring (n_row, n_col, Row, Col, A, p, knobs,
00449                 &n_row2, &n_col2, &max_deg) ;
00450   
00451   /* === Order the supercolumns =========================================== */
00452   
00453   ngarbage = Eigen::internal::find_ordering (n_row, n_col, Alen, Row, Col, A, p,
00454                             n_col2, max_deg, 2*nnz) ;
00455   
00456   /* === Order the non-principal columns ================================== */
00457   
00458   Eigen::internal::order_children (n_col, Col, p) ;
00459   
00460   /* === Return statistics in stats ======================================= */
00461   
00462   stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
00463   stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
00464   stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
00465   COLAMD_DEBUG0 (("colamd: done.\n")) ; 
00466   return (true) ;
00467 }
00468 
00469 /* ========================================================================== */
00470 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
00471 /* ========================================================================== */
00472 
00473 /* There are no user-callable routines beyond this point in the file */
00474 
00475 
00476 /* ========================================================================== */
00477 /* === init_rows_cols ======================================================= */
00478 /* ========================================================================== */
00479 
00480 /*
00481   Takes the column form of the matrix in A and creates the row form of the
00482   matrix.  Also, row and column attributes are stored in the Col and Row
00483   structs.  If the columns are un-sorted or contain duplicate row indices,
00484   this routine will also sort and remove duplicate row indices from the
00485   column form of the matrix.  Returns false if the matrix is invalid,
00486   true otherwise.  Not user-callable.
00487 */
00488 template <typename Index>
00489 static Index init_rows_cols  /* returns true if OK, or false otherwise */
00490   (
00491     /* === Parameters ======================================================= */
00492 
00493     Index n_row,      /* number of rows of A */
00494     Index n_col,      /* number of columns of A */
00495     Colamd_Row<Index> Row [],    /* of size n_row+1 */
00496     colamd_col<Index> Col [],    /* of size n_col+1 */
00497     Index A [],     /* row indices of A, of size Alen */
00498     Index p [],     /* pointers to columns in A, of size n_col+1 */
00499     Index stats [COLAMD_STATS]  /* colamd statistics */ 
00500     )
00501 {
00502   /* === Local variables ================================================== */
00503 
00504   Index col ;     /* a column index */
00505   Index row ;     /* a row index */
00506   Index *cp ;     /* a column pointer */
00507   Index *cp_end ;   /* a pointer to the end of a column */
00508   Index *rp ;     /* a row pointer */
00509   Index *rp_end ;   /* a pointer to the end of a row */
00510   Index last_row ;    /* previous row */
00511 
00512   /* === Initialize columns, and check column pointers ==================== */
00513 
00514   for (col = 0 ; col < n_col ; col++)
00515   {
00516     Col [col].start = p [col] ;
00517     Col [col].length = p [col+1] - p [col] ;
00518 
00519     if (Col [col].length < 0)
00520     {
00521       /* column pointers must be non-decreasing */
00522       stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
00523       stats [COLAMD_INFO1] = col ;
00524       stats [COLAMD_INFO2] = Col [col].length ;
00525       COLAMD_DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
00526       return (false) ;
00527     }
00528 
00529     Col [col].shared1.thickness = 1 ;
00530     Col [col].shared2.score = 0 ;
00531     Col [col].shared3.prev = COLAMD_EMPTY ;
00532     Col [col].shared4.degree_next = COLAMD_EMPTY ;
00533   }
00534 
00535   /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
00536 
00537   /* === Scan columns, compute row degrees, and check row indices ========= */
00538 
00539   stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
00540 
00541   for (row = 0 ; row < n_row ; row++)
00542   {
00543     Row [row].length = 0 ;
00544     Row [row].shared2.mark = -1 ;
00545   }
00546 
00547   for (col = 0 ; col < n_col ; col++)
00548   {
00549     last_row = -1 ;
00550 
00551     cp = &A [p [col]] ;
00552     cp_end = &A [p [col+1]] ;
00553 
00554     while (cp < cp_end)
00555     {
00556       row = *cp++ ;
00557 
00558       /* make sure row indices within range */
00559       if (row < 0 || row >= n_row)
00560       {
00561         stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
00562         stats [COLAMD_INFO1] = col ;
00563         stats [COLAMD_INFO2] = row ;
00564         stats [COLAMD_INFO3] = n_row ;
00565         COLAMD_DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
00566         return (false) ;
00567       }
00568 
00569       if (row <= last_row || Row [row].shared2.mark == col)
00570       {
00571         /* row index are unsorted or repeated (or both), thus col */
00572         /* is jumbled.  This is a notice, not an error condition. */
00573         stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
00574         stats [COLAMD_INFO1] = col ;
00575         stats [COLAMD_INFO2] = row ;
00576         (stats [COLAMD_INFO3]) ++ ;
00577         COLAMD_DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
00578       }
00579 
00580       if (Row [row].shared2.mark != col)
00581       {
00582         Row [row].length++ ;
00583       }
00584       else
00585       {
00586         /* this is a repeated entry in the column, */
00587         /* it will be removed */
00588         Col [col].length-- ;
00589       }
00590 
00591       /* mark the row as having been seen in this column */
00592       Row [row].shared2.mark = col ;
00593 
00594       last_row = row ;
00595     }
00596   }
00597 
00598   /* === Compute row pointers ============================================= */
00599 
00600   /* row form of the matrix starts directly after the column */
00601   /* form of matrix in A */
00602   Row [0].start = p [n_col] ;
00603   Row [0].shared1.p = Row [0].start ;
00604   Row [0].shared2.mark = -1 ;
00605   for (row = 1 ; row < n_row ; row++)
00606   {
00607     Row [row].start = Row [row-1].start + Row [row-1].length ;
00608     Row [row].shared1.p = Row [row].start ;
00609     Row [row].shared2.mark = -1 ;
00610   }
00611 
00612   /* === Create row form ================================================== */
00613 
00614   if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
00615   {
00616     /* if cols jumbled, watch for repeated row indices */
00617     for (col = 0 ; col < n_col ; col++)
00618     {
00619       cp = &A [p [col]] ;
00620       cp_end = &A [p [col+1]] ;
00621       while (cp < cp_end)
00622       {
00623         row = *cp++ ;
00624         if (Row [row].shared2.mark != col)
00625         {
00626           A [(Row [row].shared1.p)++] = col ;
00627           Row [row].shared2.mark = col ;
00628         }
00629       }
00630     }
00631   }
00632   else
00633   {
00634     /* if cols not jumbled, we don't need the mark (this is faster) */
00635     for (col = 0 ; col < n_col ; col++)
00636     {
00637       cp = &A [p [col]] ;
00638       cp_end = &A [p [col+1]] ;
00639       while (cp < cp_end)
00640       {
00641         A [(Row [*cp++].shared1.p)++] = col ;
00642       }
00643     }
00644   }
00645 
00646   /* === Clear the row marks and set row degrees ========================== */
00647 
00648   for (row = 0 ; row < n_row ; row++)
00649   {
00650     Row [row].shared2.mark = 0 ;
00651     Row [row].shared1.degree = Row [row].length ;
00652   }
00653 
00654   /* === See if we need to re-create columns ============================== */
00655 
00656   if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
00657   {
00658     COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
00659 
00660 
00661     /* === Compute col pointers ========================================= */
00662 
00663     /* col form of the matrix starts at A [0]. */
00664     /* Note, we may have a gap between the col form and the row */
00665     /* form if there were duplicate entries, if so, it will be */
00666     /* removed upon the first garbage collection */
00667     Col [0].start = 0 ;
00668     p [0] = Col [0].start ;
00669     for (col = 1 ; col < n_col ; col++)
00670     {
00671       /* note that the lengths here are for pruned columns, i.e. */
00672       /* no duplicate row indices will exist for these columns */
00673       Col [col].start = Col [col-1].start + Col [col-1].length ;
00674       p [col] = Col [col].start ;
00675     }
00676 
00677     /* === Re-create col form =========================================== */
00678 
00679     for (row = 0 ; row < n_row ; row++)
00680     {
00681       rp = &A [Row [row].start] ;
00682       rp_end = rp + Row [row].length ;
00683       while (rp < rp_end)
00684       {
00685         A [(p [*rp++])++] = row ;
00686       }
00687     }
00688   }
00689 
00690   /* === Done.  Matrix is not (or no longer) jumbled ====================== */
00691 
00692   return (true) ;
00693 }
00694 
00695 
00696 /* ========================================================================== */
00697 /* === init_scoring ========================================================= */
00698 /* ========================================================================== */
00699 
00700 /*
00701   Kills dense or empty columns and rows, calculates an initial score for
00702   each column, and places all columns in the degree lists.  Not user-callable.
00703 */
00704 template <typename Index>
00705 static void init_scoring
00706   (
00707     /* === Parameters ======================================================= */
00708 
00709     Index n_row,      /* number of rows of A */
00710     Index n_col,      /* number of columns of A */
00711     Colamd_Row<Index> Row [],    /* of size n_row+1 */
00712     colamd_col<Index> Col [],    /* of size n_col+1 */
00713     Index A [],     /* column form and row form of A */
00714     Index head [],    /* of size n_col+1 */
00715     double knobs [COLAMD_KNOBS],/* parameters */
00716     Index *p_n_row2,    /* number of non-dense, non-empty rows */
00717     Index *p_n_col2,    /* number of non-dense, non-empty columns */
00718     Index *p_max_deg    /* maximum row degree */
00719     )
00720 {
00721   /* === Local variables ================================================== */
00722 
00723   Index c ;     /* a column index */
00724   Index r, row ;    /* a row index */
00725   Index *cp ;     /* a column pointer */
00726   Index deg ;     /* degree of a row or column */
00727   Index *cp_end ;   /* a pointer to the end of a column */
00728   Index *new_cp ;   /* new column pointer */
00729   Index col_length ;    /* length of pruned column */
00730   Index score ;     /* current column score */
00731   Index n_col2 ;    /* number of non-dense, non-empty columns */
00732   Index n_row2 ;    /* number of non-dense, non-empty rows */
00733   Index dense_row_count ; /* remove rows with more entries than this */
00734   Index dense_col_count ; /* remove cols with more entries than this */
00735   Index min_score ;   /* smallest column score */
00736   Index max_deg ;   /* maximum row degree */
00737   Index next_col ;    /* Used to add to degree list.*/
00738 
00739 
00740   /* === Extract knobs ==================================================== */
00741 
00742   dense_row_count = COLAMD_MAX (0, COLAMD_MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ;
00743   dense_col_count = COLAMD_MAX (0, COLAMD_MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ;
00744   COLAMD_DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
00745   max_deg = 0 ;
00746   n_col2 = n_col ;
00747   n_row2 = n_row ;
00748 
00749   /* === Kill empty columns =============================================== */
00750 
00751   /* Put the empty columns at the end in their natural order, so that LU */
00752   /* factorization can proceed as far as possible. */
00753   for (c = n_col-1 ; c >= 0 ; c--)
00754   {
00755     deg = Col [c].length ;
00756     if (deg == 0)
00757     {
00758       /* this is a empty column, kill and order it last */
00759       Col [c].shared2.order = --n_col2 ;
00760       KILL_PRINCIPAL_COL (c) ;
00761     }
00762   }
00763   COLAMD_DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
00764 
00765   /* === Kill dense columns =============================================== */
00766 
00767   /* Put the dense columns at the end, in their natural order */
00768   for (c = n_col-1 ; c >= 0 ; c--)
00769   {
00770     /* skip any dead columns */
00771     if (COL_IS_DEAD (c))
00772     {
00773       continue ;
00774     }
00775     deg = Col [c].length ;
00776     if (deg > dense_col_count)
00777     {
00778       /* this is a dense column, kill and order it last */
00779       Col [c].shared2.order = --n_col2 ;
00780       /* decrement the row degrees */
00781       cp = &A [Col [c].start] ;
00782       cp_end = cp + Col [c].length ;
00783       while (cp < cp_end)
00784       {
00785         Row [*cp++].shared1.degree-- ;
00786       }
00787       KILL_PRINCIPAL_COL (c) ;
00788     }
00789   }
00790   COLAMD_DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
00791 
00792   /* === Kill dense and empty rows ======================================== */
00793 
00794   for (r = 0 ; r < n_row ; r++)
00795   {
00796     deg = Row [r].shared1.degree ;
00797     COLAMD_ASSERT (deg >= 0 && deg <= n_col) ;
00798     if (deg > dense_row_count || deg == 0)
00799     {
00800       /* kill a dense or empty row */
00801       KILL_ROW (r) ;
00802       --n_row2 ;
00803     }
00804     else
00805     {
00806       /* keep track of max degree of remaining rows */
00807       max_deg = COLAMD_MAX (max_deg, deg) ;
00808     }
00809   }
00810   COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
00811 
00812   /* === Compute initial column scores ==================================== */
00813 
00814   /* At this point the row degrees are accurate.  They reflect the number */
00815   /* of "live" (non-dense) columns in each row.  No empty rows exist. */
00816   /* Some "live" columns may contain only dead rows, however.  These are */
00817   /* pruned in the code below. */
00818 
00819   /* now find the initial matlab score for each column */
00820   for (c = n_col-1 ; c >= 0 ; c--)
00821   {
00822     /* skip dead column */
00823     if (COL_IS_DEAD (c))
00824     {
00825       continue ;
00826     }
00827     score = 0 ;
00828     cp = &A [Col [c].start] ;
00829     new_cp = cp ;
00830     cp_end = cp + Col [c].length ;
00831     while (cp < cp_end)
00832     {
00833       /* get a row */
00834       row = *cp++ ;
00835       /* skip if dead */
00836       if (ROW_IS_DEAD (row))
00837       {
00838         continue ;
00839       }
00840       /* compact the column */
00841       *new_cp++ = row ;
00842       /* add row's external degree */
00843       score += Row [row].shared1.degree - 1 ;
00844       /* guard against integer overflow */
00845       score = COLAMD_MIN (score, n_col) ;
00846     }
00847     /* determine pruned column length */
00848     col_length = (Index) (new_cp - &A [Col [c].start]) ;
00849     if (col_length == 0)
00850     {
00851       /* a newly-made null column (all rows in this col are "dense" */
00852       /* and have already been killed) */
00853       COLAMD_DEBUG2 (("Newly null killed: %d\n", c)) ;
00854       Col [c].shared2.order = --n_col2 ;
00855       KILL_PRINCIPAL_COL (c) ;
00856     }
00857     else
00858     {
00859       /* set column length and set score */
00860       COLAMD_ASSERT (score >= 0) ;
00861       COLAMD_ASSERT (score <= n_col) ;
00862       Col [c].length = col_length ;
00863       Col [c].shared2.score = score ;
00864     }
00865   }
00866   COLAMD_DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
00867                   n_col-n_col2)) ;
00868 
00869   /* At this point, all empty rows and columns are dead.  All live columns */
00870   /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
00871   /* yet).  Rows may contain dead columns, but all live rows contain at */
00872   /* least one live column. */
00873 
00874   /* === Initialize degree lists ========================================== */
00875 
00876 
00877   /* clear the hash buckets */
00878   for (c = 0 ; c <= n_col ; c++)
00879   {
00880     head [c] = COLAMD_EMPTY ;
00881   }
00882   min_score = n_col ;
00883   /* place in reverse order, so low column indices are at the front */
00884   /* of the lists.  This is to encourage natural tie-breaking */
00885   for (c = n_col-1 ; c >= 0 ; c--)
00886   {
00887     /* only add principal columns to degree lists */
00888     if (COL_IS_ALIVE (c))
00889     {
00890       COLAMD_DEBUG4 (("place %d score %d minscore %d ncol %d\n",
00891                       c, Col [c].shared2.score, min_score, n_col)) ;
00892 
00893       /* === Add columns score to DList =============================== */
00894 
00895       score = Col [c].shared2.score ;
00896 
00897       COLAMD_ASSERT (min_score >= 0) ;
00898       COLAMD_ASSERT (min_score <= n_col) ;
00899       COLAMD_ASSERT (score >= 0) ;
00900       COLAMD_ASSERT (score <= n_col) ;
00901       COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ;
00902 
00903       /* now add this column to dList at proper score location */
00904       next_col = head [score] ;
00905       Col [c].shared3.prev = COLAMD_EMPTY ;
00906       Col [c].shared4.degree_next = next_col ;
00907 
00908       /* if there already was a column with the same score, set its */
00909       /* previous pointer to this new column */
00910       if (next_col != COLAMD_EMPTY)
00911       {
00912         Col [next_col].shared3.prev = c ;
00913       }
00914       head [score] = c ;
00915 
00916       /* see if this score is less than current min */
00917       min_score = COLAMD_MIN (min_score, score) ;
00918 
00919 
00920     }
00921   }
00922 
00923 
00924   /* === Return number of remaining columns, and max row degree =========== */
00925 
00926   *p_n_col2 = n_col2 ;
00927   *p_n_row2 = n_row2 ;
00928   *p_max_deg = max_deg ;
00929 }
00930 
00931 
00932 /* ========================================================================== */
00933 /* === find_ordering ======================================================== */
00934 /* ========================================================================== */
00935 
00936 /*
00937   Order the principal columns of the supercolumn form of the matrix
00938   (no supercolumns on input).  Uses a minimum approximate column minimum
00939   degree ordering method.  Not user-callable.
00940 */
00941 template <typename Index>
00942 static Index find_ordering /* return the number of garbage collections */
00943   (
00944     /* === Parameters ======================================================= */
00945 
00946     Index n_row,      /* number of rows of A */
00947     Index n_col,      /* number of columns of A */
00948     Index Alen,     /* size of A, 2*nnz + n_col or larger */
00949     Colamd_Row<Index> Row [],    /* of size n_row+1 */
00950     colamd_col<Index> Col [],    /* of size n_col+1 */
00951     Index A [],     /* column form and row form of A */
00952     Index head [],    /* of size n_col+1 */
00953     Index n_col2,     /* Remaining columns to order */
00954     Index max_deg,    /* Maximum row degree */
00955     Index pfree     /* index of first free slot (2*nnz on entry) */
00956     )
00957 {
00958   /* === Local variables ================================================== */
00959 
00960   Index k ;     /* current pivot ordering step */
00961   Index pivot_col ;   /* current pivot column */
00962   Index *cp ;     /* a column pointer */
00963   Index *rp ;     /* a row pointer */
00964   Index pivot_row ;   /* current pivot row */
00965   Index *new_cp ;   /* modified column pointer */
00966   Index *new_rp ;   /* modified row pointer */
00967   Index pivot_row_start ; /* pointer to start of pivot row */
00968   Index pivot_row_degree ;  /* number of columns in pivot row */
00969   Index pivot_row_length ;  /* number of supercolumns in pivot row */
00970   Index pivot_col_score ; /* score of pivot column */
00971   Index needed_memory ;   /* free space needed for pivot row */
00972   Index *cp_end ;   /* pointer to the end of a column */
00973   Index *rp_end ;   /* pointer to the end of a row */
00974   Index row ;     /* a row index */
00975   Index col ;     /* a column index */
00976   Index max_score ;   /* maximum possible score */
00977   Index cur_score ;   /* score of current column */
00978   unsigned int hash ;   /* hash value for supernode detection */
00979   Index head_column ;   /* head of hash bucket */
00980   Index first_col ;   /* first column in hash bucket */
00981   Index tag_mark ;    /* marker value for mark array */
00982   Index row_mark ;    /* Row [row].shared2.mark */
00983   Index set_difference ;  /* set difference size of row with pivot row */
00984   Index min_score ;   /* smallest column score */
00985   Index col_thickness ;   /* "thickness" (no. of columns in a supercol) */
00986   Index max_mark ;    /* maximum value of tag_mark */
00987   Index pivot_col_thickness ; /* number of columns represented by pivot col */
00988   Index prev_col ;    /* Used by Dlist operations. */
00989   Index next_col ;    /* Used by Dlist operations. */
00990   Index ngarbage ;    /* number of garbage collections performed */
00991 
00992 
00993   /* === Initialization and clear mark ==================================== */
00994 
00995   max_mark = INT_MAX - n_col ;  /* INT_MAX defined in <limits.h> */
00996   tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
00997   min_score = 0 ;
00998   ngarbage = 0 ;
00999   COLAMD_DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
01000 
01001   /* === Order the columns ================================================ */
01002 
01003   for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
01004   {
01005 
01006     /* === Select pivot column, and order it ============================ */
01007 
01008     /* make sure degree list isn't empty */
01009     COLAMD_ASSERT (min_score >= 0) ;
01010     COLAMD_ASSERT (min_score <= n_col) ;
01011     COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ;
01012 
01013     /* get pivot column from head of minimum degree list */
01014     while (head [min_score] == COLAMD_EMPTY && min_score < n_col)
01015     {
01016       min_score++ ;
01017     }
01018     pivot_col = head [min_score] ;
01019     COLAMD_ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
01020     next_col = Col [pivot_col].shared4.degree_next ;
01021     head [min_score] = next_col ;
01022     if (next_col != COLAMD_EMPTY)
01023     {
01024       Col [next_col].shared3.prev = COLAMD_EMPTY ;
01025     }
01026 
01027     COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ;
01028     COLAMD_DEBUG3 (("Pivot col: %d\n", pivot_col)) ;
01029 
01030     /* remember score for defrag check */
01031     pivot_col_score = Col [pivot_col].shared2.score ;
01032 
01033     /* the pivot column is the kth column in the pivot order */
01034     Col [pivot_col].shared2.order = k ;
01035 
01036     /* increment order count by column thickness */
01037     pivot_col_thickness = Col [pivot_col].shared1.thickness ;
01038     k += pivot_col_thickness ;
01039     COLAMD_ASSERT (pivot_col_thickness > 0) ;
01040 
01041     /* === Garbage_collection, if necessary ============================= */
01042 
01043     needed_memory = COLAMD_MIN (pivot_col_score, n_col - k) ;
01044     if (pfree + needed_memory >= Alen)
01045     {
01046       pfree = Eigen::internal::garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
01047       ngarbage++ ;
01048       /* after garbage collection we will have enough */
01049       COLAMD_ASSERT (pfree + needed_memory < Alen) ;
01050       /* garbage collection has wiped out the Row[].shared2.mark array */
01051       tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
01052 
01053     }
01054 
01055     /* === Compute pivot row pattern ==================================== */
01056 
01057     /* get starting location for this new merged row */
01058     pivot_row_start = pfree ;
01059 
01060     /* initialize new row counts to zero */
01061     pivot_row_degree = 0 ;
01062 
01063     /* tag pivot column as having been visited so it isn't included */
01064     /* in merged pivot row */
01065     Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
01066 
01067     /* pivot row is the union of all rows in the pivot column pattern */
01068     cp = &A [Col [pivot_col].start] ;
01069     cp_end = cp + Col [pivot_col].length ;
01070     while (cp < cp_end)
01071     {
01072       /* get a row */
01073       row = *cp++ ;
01074       COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
01075       /* skip if row is dead */
01076       if (ROW_IS_DEAD (row))
01077       {
01078         continue ;
01079       }
01080       rp = &A [Row [row].start] ;
01081       rp_end = rp + Row [row].length ;
01082       while (rp < rp_end)
01083       {
01084         /* get a column */
01085         col = *rp++ ;
01086         /* add the column, if alive and untagged */
01087         col_thickness = Col [col].shared1.thickness ;
01088         if (col_thickness > 0 && COL_IS_ALIVE (col))
01089         {
01090           /* tag column in pivot row */
01091           Col [col].shared1.thickness = -col_thickness ;
01092           COLAMD_ASSERT (pfree < Alen) ;
01093           /* place column in pivot row */
01094           A [pfree++] = col ;
01095           pivot_row_degree += col_thickness ;
01096         }
01097       }
01098     }
01099 
01100     /* clear tag on pivot column */
01101     Col [pivot_col].shared1.thickness = pivot_col_thickness ;
01102     max_deg = COLAMD_MAX (max_deg, pivot_row_degree) ;
01103 
01104 
01105     /* === Kill all rows used to construct pivot row ==================== */
01106 
01107     /* also kill pivot row, temporarily */
01108     cp = &A [Col [pivot_col].start] ;
01109     cp_end = cp + Col [pivot_col].length ;
01110     while (cp < cp_end)
01111     {
01112       /* may be killing an already dead row */
01113       row = *cp++ ;
01114       COLAMD_DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
01115       KILL_ROW (row) ;
01116     }
01117 
01118     /* === Select a row index to use as the new pivot row =============== */
01119 
01120     pivot_row_length = pfree - pivot_row_start ;
01121     if (pivot_row_length > 0)
01122     {
01123       /* pick the "pivot" row arbitrarily (first row in col) */
01124       pivot_row = A [Col [pivot_col].start] ;
01125       COLAMD_DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
01126     }
01127     else
01128     {
01129       /* there is no pivot row, since it is of zero length */
01130       pivot_row = COLAMD_EMPTY ;
01131       COLAMD_ASSERT (pivot_row_length == 0) ;
01132     }
01133     COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
01134 
01135     /* === Approximate degree computation =============================== */
01136 
01137     /* Here begins the computation of the approximate degree.  The column */
01138     /* score is the sum of the pivot row "length", plus the size of the */
01139     /* set differences of each row in the column minus the pattern of the */
01140     /* pivot row itself.  The column ("thickness") itself is also */
01141     /* excluded from the column score (we thus use an approximate */
01142     /* external degree). */
01143 
01144     /* The time taken by the following code (compute set differences, and */
01145     /* add them up) is proportional to the size of the data structure */
01146     /* being scanned - that is, the sum of the sizes of each column in */
01147     /* the pivot row.  Thus, the amortized time to compute a column score */
01148     /* is proportional to the size of that column (where size, in this */
01149     /* context, is the column "length", or the number of row indices */
01150     /* in that column).  The number of row indices in a column is */
01151     /* monotonically non-decreasing, from the length of the original */
01152     /* column on input to colamd. */
01153 
01154     /* === Compute set differences ====================================== */
01155 
01156     COLAMD_DEBUG3 (("** Computing set differences phase. **\n")) ;
01157 
01158     /* pivot row is currently dead - it will be revived later. */
01159 
01160     COLAMD_DEBUG3 (("Pivot row: ")) ;
01161     /* for each column in pivot row */
01162     rp = &A [pivot_row_start] ;
01163     rp_end = rp + pivot_row_length ;
01164     while (rp < rp_end)
01165     {
01166       col = *rp++ ;
01167       COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
01168       COLAMD_DEBUG3 (("Col: %d\n", col)) ;
01169 
01170       /* clear tags used to construct pivot row pattern */
01171       col_thickness = -Col [col].shared1.thickness ;
01172       COLAMD_ASSERT (col_thickness > 0) ;
01173       Col [col].shared1.thickness = col_thickness ;
01174 
01175       /* === Remove column from degree list =========================== */
01176 
01177       cur_score = Col [col].shared2.score ;
01178       prev_col = Col [col].shared3.prev ;
01179       next_col = Col [col].shared4.degree_next ;
01180       COLAMD_ASSERT (cur_score >= 0) ;
01181       COLAMD_ASSERT (cur_score <= n_col) ;
01182       COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ;
01183       if (prev_col == COLAMD_EMPTY)
01184       {
01185         head [cur_score] = next_col ;
01186       }
01187       else
01188       {
01189         Col [prev_col].shared4.degree_next = next_col ;
01190       }
01191       if (next_col != COLAMD_EMPTY)
01192       {
01193         Col [next_col].shared3.prev = prev_col ;
01194       }
01195 
01196       /* === Scan the column ========================================== */
01197 
01198       cp = &A [Col [col].start] ;
01199       cp_end = cp + Col [col].length ;
01200       while (cp < cp_end)
01201       {
01202         /* get a row */
01203         row = *cp++ ;
01204         row_mark = Row [row].shared2.mark ;
01205         /* skip if dead */
01206         if (ROW_IS_MARKED_DEAD (row_mark))
01207         {
01208           continue ;
01209         }
01210         COLAMD_ASSERT (row != pivot_row) ;
01211         set_difference = row_mark - tag_mark ;
01212         /* check if the row has been seen yet */
01213         if (set_difference < 0)
01214         {
01215           COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ;
01216           set_difference = Row [row].shared1.degree ;
01217         }
01218         /* subtract column thickness from this row's set difference */
01219         set_difference -= col_thickness ;
01220         COLAMD_ASSERT (set_difference >= 0) ;
01221         /* absorb this row if the set difference becomes zero */
01222         if (set_difference == 0)
01223         {
01224           COLAMD_DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
01225           KILL_ROW (row) ;
01226         }
01227         else
01228         {
01229           /* save the new mark */
01230           Row [row].shared2.mark = set_difference + tag_mark ;
01231         }
01232       }
01233     }
01234 
01235 
01236     /* === Add up set differences for each column ======================= */
01237 
01238     COLAMD_DEBUG3 (("** Adding set differences phase. **\n")) ;
01239 
01240     /* for each column in pivot row */
01241     rp = &A [pivot_row_start] ;
01242     rp_end = rp + pivot_row_length ;
01243     while (rp < rp_end)
01244     {
01245       /* get a column */
01246       col = *rp++ ;
01247       COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
01248       hash = 0 ;
01249       cur_score = 0 ;
01250       cp = &A [Col [col].start] ;
01251       /* compact the column */
01252       new_cp = cp ;
01253       cp_end = cp + Col [col].length ;
01254 
01255       COLAMD_DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
01256 
01257       while (cp < cp_end)
01258       {
01259         /* get a row */
01260         row = *cp++ ;
01261         COLAMD_ASSERT(row >= 0 && row < n_row) ;
01262         row_mark = Row [row].shared2.mark ;
01263         /* skip if dead */
01264         if (ROW_IS_MARKED_DEAD (row_mark))
01265         {
01266           continue ;
01267         }
01268         COLAMD_ASSERT (row_mark > tag_mark) ;
01269         /* compact the column */
01270         *new_cp++ = row ;
01271         /* compute hash function */
01272         hash += row ;
01273         /* add set difference */
01274         cur_score += row_mark - tag_mark ;
01275         /* integer overflow... */
01276         cur_score = COLAMD_MIN (cur_score, n_col) ;
01277       }
01278 
01279       /* recompute the column's length */
01280       Col [col].length = (Index) (new_cp - &A [Col [col].start]) ;
01281 
01282       /* === Further mass elimination ================================= */
01283 
01284       if (Col [col].length == 0)
01285       {
01286         COLAMD_DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
01287         /* nothing left but the pivot row in this column */
01288         KILL_PRINCIPAL_COL (col) ;
01289         pivot_row_degree -= Col [col].shared1.thickness ;
01290         COLAMD_ASSERT (pivot_row_degree >= 0) ;
01291         /* order it */
01292         Col [col].shared2.order = k ;
01293         /* increment order count by column thickness */
01294         k += Col [col].shared1.thickness ;
01295       }
01296       else
01297       {
01298         /* === Prepare for supercolumn detection ==================== */
01299 
01300         COLAMD_DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
01301 
01302         /* save score so far */
01303         Col [col].shared2.score = cur_score ;
01304 
01305         /* add column to hash table, for supercolumn detection */
01306         hash %= n_col + 1 ;
01307 
01308         COLAMD_DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
01309         COLAMD_ASSERT (hash <= n_col) ;
01310 
01311         head_column = head [hash] ;
01312         if (head_column > COLAMD_EMPTY)
01313         {
01314           /* degree list "hash" is non-empty, use prev (shared3) of */
01315           /* first column in degree list as head of hash bucket */
01316           first_col = Col [head_column].shared3.headhash ;
01317           Col [head_column].shared3.headhash = col ;
01318         }
01319         else
01320         {
01321           /* degree list "hash" is empty, use head as hash bucket */
01322           first_col = - (head_column + 2) ;
01323           head [hash] = - (col + 2) ;
01324         }
01325         Col [col].shared4.hash_next = first_col ;
01326 
01327         /* save hash function in Col [col].shared3.hash */
01328         Col [col].shared3.hash = (Index) hash ;
01329         COLAMD_ASSERT (COL_IS_ALIVE (col)) ;
01330       }
01331     }
01332 
01333     /* The approximate external column degree is now computed.  */
01334 
01335     /* === Supercolumn detection ======================================== */
01336 
01337     COLAMD_DEBUG3 (("** Supercolumn detection phase. **\n")) ;
01338 
01339     Eigen::internal::detect_super_cols (Col, A, head, pivot_row_start, pivot_row_length) ;
01340 
01341     /* === Kill the pivotal column ====================================== */
01342 
01343     KILL_PRINCIPAL_COL (pivot_col) ;
01344 
01345     /* === Clear mark =================================================== */
01346 
01347     tag_mark += (max_deg + 1) ;
01348     if (tag_mark >= max_mark)
01349     {
01350       COLAMD_DEBUG2 (("clearing tag_mark\n")) ;
01351       tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
01352     }
01353 
01354     /* === Finalize the new pivot row, and column scores ================ */
01355 
01356     COLAMD_DEBUG3 (("** Finalize scores phase. **\n")) ;
01357 
01358     /* for each column in pivot row */
01359     rp = &A [pivot_row_start] ;
01360     /* compact the pivot row */
01361     new_rp = rp ;
01362     rp_end = rp + pivot_row_length ;
01363     while (rp < rp_end)
01364     {
01365       col = *rp++ ;
01366       /* skip dead columns */
01367       if (COL_IS_DEAD (col))
01368       {
01369         continue ;
01370       }
01371       *new_rp++ = col ;
01372       /* add new pivot row to column */
01373       A [Col [col].start + (Col [col].length++)] = pivot_row ;
01374 
01375       /* retrieve score so far and add on pivot row's degree. */
01376       /* (we wait until here for this in case the pivot */
01377       /* row's degree was reduced due to mass elimination). */
01378       cur_score = Col [col].shared2.score + pivot_row_degree ;
01379 
01380       /* calculate the max possible score as the number of */
01381       /* external columns minus the 'k' value minus the */
01382       /* columns thickness */
01383       max_score = n_col - k - Col [col].shared1.thickness ;
01384 
01385       /* make the score the external degree of the union-of-rows */
01386       cur_score -= Col [col].shared1.thickness ;
01387 
01388       /* make sure score is less or equal than the max score */
01389       cur_score = COLAMD_MIN (cur_score, max_score) ;
01390       COLAMD_ASSERT (cur_score >= 0) ;
01391 
01392       /* store updated score */
01393       Col [col].shared2.score = cur_score ;
01394 
01395       /* === Place column back in degree list ========================= */
01396 
01397       COLAMD_ASSERT (min_score >= 0) ;
01398       COLAMD_ASSERT (min_score <= n_col) ;
01399       COLAMD_ASSERT (cur_score >= 0) ;
01400       COLAMD_ASSERT (cur_score <= n_col) ;
01401       COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ;
01402       next_col = head [cur_score] ;
01403       Col [col].shared4.degree_next = next_col ;
01404       Col [col].shared3.prev = COLAMD_EMPTY ;
01405       if (next_col != COLAMD_EMPTY)
01406       {
01407         Col [next_col].shared3.prev = col ;
01408       }
01409       head [cur_score] = col ;
01410 
01411       /* see if this score is less than current min */
01412       min_score = COLAMD_MIN (min_score, cur_score) ;
01413 
01414     }
01415 
01416     /* === Resurrect the new pivot row ================================== */
01417 
01418     if (pivot_row_degree > 0)
01419     {
01420       /* update pivot row length to reflect any cols that were killed */
01421       /* during super-col detection and mass elimination */
01422       Row [pivot_row].start  = pivot_row_start ;
01423       Row [pivot_row].length = (Index) (new_rp - &A[pivot_row_start]) ;
01424       Row [pivot_row].shared1.degree = pivot_row_degree ;
01425       Row [pivot_row].shared2.mark = 0 ;
01426       /* pivot row is no longer dead */
01427     }
01428   }
01429 
01430   /* === All principal columns have now been ordered ====================== */
01431 
01432   return (ngarbage) ;
01433 }
01434 
01435 
01436 /* ========================================================================== */
01437 /* === order_children ======================================================= */
01438 /* ========================================================================== */
01439 
01440 /*
01441   The find_ordering routine has ordered all of the principal columns (the
01442   representatives of the supercolumns).  The non-principal columns have not
01443   yet been ordered.  This routine orders those columns by walking up the
01444   parent tree (a column is a child of the column which absorbed it).  The
01445   final permutation vector is then placed in p [0 ... n_col-1], with p [0]
01446   being the first column, and p [n_col-1] being the last.  It doesn't look
01447   like it at first glance, but be assured that this routine takes time linear
01448   in the number of columns.  Although not immediately obvious, the time
01449   taken by this routine is O (n_col), that is, linear in the number of
01450   columns.  Not user-callable.
01451 */
01452 template <typename Index>
01453 static inline  void order_children
01454 (
01455   /* === Parameters ======================================================= */
01456 
01457   Index n_col,      /* number of columns of A */
01458   colamd_col<Index> Col [],    /* of size n_col+1 */
01459   Index p []      /* p [0 ... n_col-1] is the column permutation*/
01460   )
01461 {
01462   /* === Local variables ================================================== */
01463 
01464   Index i ;     /* loop counter for all columns */
01465   Index c ;     /* column index */
01466   Index parent ;    /* index of column's parent */
01467   Index order ;     /* column's order */
01468 
01469   /* === Order each non-principal column ================================== */
01470 
01471   for (i = 0 ; i < n_col ; i++)
01472   {
01473     /* find an un-ordered non-principal column */
01474     COLAMD_ASSERT (COL_IS_DEAD (i)) ;
01475     if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY)
01476     {
01477       parent = i ;
01478       /* once found, find its principal parent */
01479       do
01480       {
01481         parent = Col [parent].shared1.parent ;
01482       } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
01483 
01484       /* now, order all un-ordered non-principal columns along path */
01485       /* to this parent.  collapse tree at the same time */
01486       c = i ;
01487       /* get order of parent */
01488       order = Col [parent].shared2.order ;
01489 
01490       do
01491       {
01492         COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ;
01493 
01494         /* order this column */
01495         Col [c].shared2.order = order++ ;
01496         /* collaps tree */
01497         Col [c].shared1.parent = parent ;
01498 
01499         /* get immediate parent of this column */
01500         c = Col [c].shared1.parent ;
01501 
01502         /* continue until we hit an ordered column.  There are */
01503         /* guarranteed not to be anymore unordered columns */
01504         /* above an ordered column */
01505       } while (Col [c].shared2.order == COLAMD_EMPTY) ;
01506 
01507       /* re-order the super_col parent to largest order for this group */
01508       Col [parent].shared2.order = order ;
01509     }
01510   }
01511 
01512   /* === Generate the permutation ========================================= */
01513 
01514   for (c = 0 ; c < n_col ; c++)
01515   {
01516     p [Col [c].shared2.order] = c ;
01517   }
01518 }
01519 
01520 
01521 /* ========================================================================== */
01522 /* === detect_super_cols ==================================================== */
01523 /* ========================================================================== */
01524 
01525 /*
01526   Detects supercolumns by finding matches between columns in the hash buckets.
01527   Check amongst columns in the set A [row_start ... row_start + row_length-1].
01528   The columns under consideration are currently *not* in the degree lists,
01529   and have already been placed in the hash buckets.
01530 
01531   The hash bucket for columns whose hash function is equal to h is stored
01532   as follows:
01533 
01534   if head [h] is >= 0, then head [h] contains a degree list, so:
01535 
01536   head [h] is the first column in degree bucket h.
01537   Col [head [h]].headhash gives the first column in hash bucket h.
01538 
01539   otherwise, the degree list is empty, and:
01540 
01541   -(head [h] + 2) is the first column in hash bucket h.
01542 
01543   For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
01544   column" pointer.  Col [c].shared3.hash is used instead as the hash number
01545   for that column.  The value of Col [c].shared4.hash_next is the next column
01546   in the same hash bucket.
01547 
01548   Assuming no, or "few" hash collisions, the time taken by this routine is
01549   linear in the sum of the sizes (lengths) of each column whose score has
01550   just been computed in the approximate degree computation.
01551   Not user-callable.
01552 */
01553 template <typename Index>
01554 static void detect_super_cols
01555 (
01556   /* === Parameters ======================================================= */
01557   
01558   colamd_col<Index> Col [],    /* of size n_col+1 */
01559   Index A [],     /* row indices of A */
01560   Index head [],    /* head of degree lists and hash buckets */
01561   Index row_start,    /* pointer to set of columns to check */
01562   Index row_length    /* number of columns to check */
01563 )
01564 {
01565   /* === Local variables ================================================== */
01566 
01567   Index hash ;      /* hash value for a column */
01568   Index *rp ;     /* pointer to a row */
01569   Index c ;     /* a column index */
01570   Index super_c ;   /* column index of the column to absorb into */
01571   Index *cp1 ;      /* column pointer for column super_c */
01572   Index *cp2 ;      /* column pointer for column c */
01573   Index length ;    /* length of column super_c */
01574   Index prev_c ;    /* column preceding c in hash bucket */
01575   Index i ;     /* loop counter */
01576   Index *rp_end ;   /* pointer to the end of the row */
01577   Index col ;     /* a column index in the row to check */
01578   Index head_column ;   /* first column in hash bucket or degree list */
01579   Index first_col ;   /* first column in hash bucket */
01580 
01581   /* === Consider each column in the row ================================== */
01582 
01583   rp = &A [row_start] ;
01584   rp_end = rp + row_length ;
01585   while (rp < rp_end)
01586   {
01587     col = *rp++ ;
01588     if (COL_IS_DEAD (col))
01589     {
01590       continue ;
01591     }
01592 
01593     /* get hash number for this column */
01594     hash = Col [col].shared3.hash ;
01595     COLAMD_ASSERT (hash <= n_col) ;
01596 
01597     /* === Get the first column in this hash bucket ===================== */
01598 
01599     head_column = head [hash] ;
01600     if (head_column > COLAMD_EMPTY)
01601     {
01602       first_col = Col [head_column].shared3.headhash ;
01603     }
01604     else
01605     {
01606       first_col = - (head_column + 2) ;
01607     }
01608 
01609     /* === Consider each column in the hash bucket ====================== */
01610 
01611     for (super_c = first_col ; super_c != COLAMD_EMPTY ;
01612          super_c = Col [super_c].shared4.hash_next)
01613     {
01614       COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ;
01615       COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ;
01616       length = Col [super_c].length ;
01617 
01618       /* prev_c is the column preceding column c in the hash bucket */
01619       prev_c = super_c ;
01620 
01621       /* === Compare super_c with all columns after it ================ */
01622 
01623       for (c = Col [super_c].shared4.hash_next ;
01624            c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next)
01625       {
01626         COLAMD_ASSERT (c != super_c) ;
01627         COLAMD_ASSERT (COL_IS_ALIVE (c)) ;
01628         COLAMD_ASSERT (Col [c].shared3.hash == hash) ;
01629 
01630         /* not identical if lengths or scores are different */
01631         if (Col [c].length != length ||
01632             Col [c].shared2.score != Col [super_c].shared2.score)
01633         {
01634           prev_c = c ;
01635           continue ;
01636         }
01637 
01638         /* compare the two columns */
01639         cp1 = &A [Col [super_c].start] ;
01640         cp2 = &A [Col [c].start] ;
01641 
01642         for (i = 0 ; i < length ; i++)
01643         {
01644           /* the columns are "clean" (no dead rows) */
01645           COLAMD_ASSERT (ROW_IS_ALIVE (*cp1))  ;
01646           COLAMD_ASSERT (ROW_IS_ALIVE (*cp2))  ;
01647           /* row indices will same order for both supercols, */
01648           /* no gather scatter nessasary */
01649           if (*cp1++ != *cp2++)
01650           {
01651             break ;
01652           }
01653         }
01654 
01655         /* the two columns are different if the for-loop "broke" */
01656         if (i != length)
01657         {
01658           prev_c = c ;
01659           continue ;
01660         }
01661 
01662         /* === Got it!  two columns are identical =================== */
01663 
01664         COLAMD_ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
01665 
01666         Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
01667         Col [c].shared1.parent = super_c ;
01668         KILL_NON_PRINCIPAL_COL (c) ;
01669         /* order c later, in order_children() */
01670         Col [c].shared2.order = COLAMD_EMPTY ;
01671         /* remove c from hash bucket */
01672         Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
01673       }
01674     }
01675 
01676     /* === Empty this hash bucket ======================================= */
01677 
01678     if (head_column > COLAMD_EMPTY)
01679     {
01680       /* corresponding degree list "hash" is not empty */
01681       Col [head_column].shared3.headhash = COLAMD_EMPTY ;
01682     }
01683     else
01684     {
01685       /* corresponding degree list "hash" is empty */
01686       head [hash] = COLAMD_EMPTY ;
01687     }
01688   }
01689 }
01690 
01691 
01692 /* ========================================================================== */
01693 /* === garbage_collection =================================================== */
01694 /* ========================================================================== */
01695 
01696 /*
01697   Defragments and compacts columns and rows in the workspace A.  Used when
01698   all avaliable memory has been used while performing row merging.  Returns
01699   the index of the first free position in A, after garbage collection.  The
01700   time taken by this routine is linear is the size of the array A, which is
01701   itself linear in the number of nonzeros in the input matrix.
01702   Not user-callable.
01703 */
01704 template <typename Index>
01705 static Index garbage_collection  /* returns the new value of pfree */
01706   (
01707     /* === Parameters ======================================================= */
01708     
01709     Index n_row,      /* number of rows */
01710     Index n_col,      /* number of columns */
01711     Colamd_Row<Index> Row [],    /* row info */
01712     colamd_col<Index> Col [],    /* column info */
01713     Index A [],     /* A [0 ... Alen-1] holds the matrix */
01714     Index *pfree      /* &A [0] ... pfree is in use */
01715     )
01716 {
01717   /* === Local variables ================================================== */
01718 
01719   Index *psrc ;     /* source pointer */
01720   Index *pdest ;    /* destination pointer */
01721   Index j ;     /* counter */
01722   Index r ;     /* a row index */
01723   Index c ;     /* a column index */
01724   Index length ;    /* length of a row or column */
01725 
01726   /* === Defragment the columns =========================================== */
01727 
01728   pdest = &A[0] ;
01729   for (c = 0 ; c < n_col ; c++)
01730   {
01731     if (COL_IS_ALIVE (c))
01732     {
01733       psrc = &A [Col [c].start] ;
01734 
01735       /* move and compact the column */
01736       COLAMD_ASSERT (pdest <= psrc) ;
01737       Col [c].start = (Index) (pdest - &A [0]) ;
01738       length = Col [c].length ;
01739       for (j = 0 ; j < length ; j++)
01740       {
01741         r = *psrc++ ;
01742         if (ROW_IS_ALIVE (r))
01743         {
01744           *pdest++ = r ;
01745         }
01746       }
01747       Col [c].length = (Index) (pdest - &A [Col [c].start]) ;
01748     }
01749   }
01750 
01751   /* === Prepare to defragment the rows =================================== */
01752 
01753   for (r = 0 ; r < n_row ; r++)
01754   {
01755     if (ROW_IS_ALIVE (r))
01756     {
01757       if (Row [r].length == 0)
01758       {
01759         /* this row is of zero length.  cannot compact it, so kill it */
01760         COLAMD_DEBUG3 (("Defrag row kill\n")) ;
01761         KILL_ROW (r) ;
01762       }
01763       else
01764       {
01765         /* save first column index in Row [r].shared2.first_column */
01766         psrc = &A [Row [r].start] ;
01767         Row [r].shared2.first_column = *psrc ;
01768         COLAMD_ASSERT (ROW_IS_ALIVE (r)) ;
01769         /* flag the start of the row with the one's complement of row */
01770         *psrc = ONES_COMPLEMENT (r) ;
01771 
01772       }
01773     }
01774   }
01775 
01776   /* === Defragment the rows ============================================== */
01777 
01778   psrc = pdest ;
01779   while (psrc < pfree)
01780   {
01781     /* find a negative number ... the start of a row */
01782     if (*psrc++ < 0)
01783     {
01784       psrc-- ;
01785       /* get the row index */
01786       r = ONES_COMPLEMENT (*psrc) ;
01787       COLAMD_ASSERT (r >= 0 && r < n_row) ;
01788       /* restore first column index */
01789       *psrc = Row [r].shared2.first_column ;
01790       COLAMD_ASSERT (ROW_IS_ALIVE (r)) ;
01791 
01792       /* move and compact the row */
01793       COLAMD_ASSERT (pdest <= psrc) ;
01794       Row [r].start = (Index) (pdest - &A [0]) ;
01795       length = Row [r].length ;
01796       for (j = 0 ; j < length ; j++)
01797       {
01798         c = *psrc++ ;
01799         if (COL_IS_ALIVE (c))
01800         {
01801           *pdest++ = c ;
01802         }
01803       }
01804       Row [r].length = (Index) (pdest - &A [Row [r].start]) ;
01805 
01806     }
01807   }
01808   /* ensure we found all the rows */
01809   COLAMD_ASSERT (debug_rows == 0) ;
01810 
01811   /* === Return the new value of pfree ==================================== */
01812 
01813   return ((Index) (pdest - &A [0])) ;
01814 }
01815 
01816 
01817 /* ========================================================================== */
01818 /* === clear_mark =========================================================== */
01819 /* ========================================================================== */
01820 
01821 /*
01822   Clears the Row [].shared2.mark array, and returns the new tag_mark.
01823   Return value is the new tag_mark.  Not user-callable.
01824 */
01825 template <typename Index>
01826 static inline  Index clear_mark  /* return the new value for tag_mark */
01827   (
01828       /* === Parameters ======================================================= */
01829 
01830     Index n_row,    /* number of rows in A */
01831     Colamd_Row<Index> Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
01832     )
01833 {
01834   /* === Local variables ================================================== */
01835 
01836   Index r ;
01837 
01838   for (r = 0 ; r < n_row ; r++)
01839   {
01840     if (ROW_IS_ALIVE (r))
01841     {
01842       Row [r].shared2.mark = 0 ;
01843     }
01844   }
01845   return (1) ;
01846 }
01847 
01848 
01849 } // namespace internal 
01850 #endif


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:58:07