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