// hnumcolamd.h
#pragma once
#ifndef __HNUMCOLAMD_H__
#define __HNUMCOLAMD_H__
#include "hnumdef.h"
#include "hwinexception.h"
namespace harlinn
{
/*
* This is a C++ implementation of the "column approximate minimum degree ordering" algorthm.
*
* This implementation is based on the original c implementation by
* Stefan I. Larimore and Timothy A. Davis
*
* colamd: An approximate minimum degree column ordering algorithm.
*
* Purpose:
*
* Colamd computes a permutation Q such that the Cholesky factorization of
* (AQ)'(AQ) has less fill-in and requires fewer floating point operations
* than A'A. This also provides a good ordering for sparse partial
* pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
* factorization, and P is computed during numerical factorization via
* conventional partial pivoting with row interchanges. Colamd is the
* column ordering method used in SuperLU, part of the ScaLAPACK library.
* It is also available as user-contributed software for Matlab 5.2,
* available from MathWorks, Inc. (http://www.mathworks.com). This
* routine can be used in place of COLMMD in Matlab. By default, the \
* and / operators in Matlab perform a column ordering (using COLMMD)
* prior to LU factorization using sparse partial pivoting, in the
* built-in Matlab LU(A) routine.
*
* Authors:
*
* The authors of the code itself are Stefan I. Larimore and Timothy A.
* Davis (davis@cise.ufl.edu), University of Florida. The algorithm was
* developed in collaboration with John Gilbert, Xerox PARC, and Esmond
* Ng, Oak Ridge National Laboratory.
*
* Date:
*
* August 3, 1998. Version 1.0.
*
* Acknowledgements:
*
* This work was supported by the National Science Foundation, under
* grants DMS-9504974 and DMS-9803599.
*
* Notice:
*
* Copyright (c) 1998 by the University of Florida. All Rights Reserved.
*
* THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
* EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
*
* Permission is hereby granted to use or copy this program for any
* purpose, provided the above notices are retained on all copies.
* User documentation of any code that uses this code must cite the
* Authors, the Copyright, and "Used by permission." If this code is
* accessible from within Matlab, then typing "help colamd" or "colamd"
* (with no arguments) must cite the Authors. Permission to modify the
* code and to distribute modified code is granted, provided the above
* notices are retained, and a notice that the code was modified is
* included with the above copyright notice. You must also retain the
* Availability information below, of the original version.
*
* This software is provided free of charge.
*
* Availability:
*
* This file is located at
*
* http://www.cise.ufl.edu/~davis/colamd/colamd.c
*
* The colamd.h file is required, located in the same directory.
* The colamdmex.c file provides a Matlab interface for colamd.
* The symamdmex.c file provides a Matlab interface for symamd, which is
* a symmetric ordering based on this code, colamd.c. All codes are
* purely ANSI C compliant (they use no Unix-specific routines, include
* files, etc.).
*
*/
namespace numerics
{
using namespace harlinn::windows;
class ColumnApproximateMinimumDegreeOrdering
{
static const int STATS = 20;
static const int DENSE_ROW = 0;
static const int DENSE_COL = 1;
static const int DEFRAG_COUNT = 2;
static const int JUMBLED_COLS = 3;
static const int EMPTY = -1;
// rowData & column status
static const int ALIVE = 0;
static const int DEAD = UINT32_MAX;
// Column status
static const int DEAD_PRINCIPAL = UINT32_MAX;
static const int DEAD_NON_PRINCIPAL = UINT32_MAX - 1;
static const int FLAG = 0x80000000;
static const int MASK = 0x7FFFFFFF;
inline static void Debug(const char*,...)
{
}
struct ColumnData
{
// index for workArea of first row in this column, or DEAD
// if column is dead
int start;
// number of rows in this column
int length;
union
{
// number of original columns represented by this column, if the column is alive
int thickness;
// parent in parent tree super-column structure, if the column is dead
int parent;
};
union
{
// the score used to maintain heap, if col is alive
int score;
// pivot ordering of this column, if col is dead
int order;
};
union
{
// head of a hash bucket, if col is at the head of a degree list
int headhash;
// hash value, if col is not in a degree list
int hash;
// previous column in degree list, if col is in a degree list (but not at the head of a degree list)
int prev;
};
union
{
// next column, if col is in a degree list
int degree_next;
// next column, if col is in a hash list
int hash_next;
};
bool IsDead() const
{
return start >= ColumnApproximateMinimumDegreeOrdering::DEAD_NON_PRINCIPAL;
}
bool IsDeadPrincipal() const
{
return start == ColumnApproximateMinimumDegreeOrdering::DEAD_PRINCIPAL;
}
bool IsAlive() const
{
return start < ColumnApproximateMinimumDegreeOrdering::DEAD_NON_PRINCIPAL;
}
void KillPrincipal() { start = ColumnApproximateMinimumDegreeOrdering::DEAD_PRINCIPAL; }
void KillNonPrincipal() { start = ColumnApproximateMinimumDegreeOrdering::DEAD_NON_PRINCIPAL; }
};
struct RowData
{
// index for workArea of first col in this row
int start;
// number of principal columns in this row
int length;
union
{
// number of principal & non-principal columns in row
int degree;
// used as a row pointer in init_rows_cols ()
int p;
};
union
{
// for computing set differences and marking dead rows
int mark;
// first column in row (used in garbage collection)
int first_column;
};
bool IsDead() const
{
return mark == ColumnApproximateMinimumDegreeOrdering::DEAD;
}
bool IsAlive() const
{
return mark < ColumnApproximateMinimumDegreeOrdering::DEAD;
}
void Kill() { mark = ColumnApproximateMinimumDegreeOrdering::DEAD; }
};
int numberOfRows;
int numberOfColumns;
int workAreaLength;
int* workArea;
ColumnData* columnData;
RowData* rowData;
double denseRow; // 0.5
double denseCol; // 0.5
public:
ColumnApproximateMinimumDegreeOrdering()
: numberOfRows(0),
numberOfColumns(0),
workAreaLength(0),
workArea(nullptr),
columnData(nullptr),
rowData(nullptr),
denseRow(0.5),
denseCol(0.5)
{}
private:
// InitializeRowsAndColumns
//
// Takes the column form of the matrix in workArea and creates the row form of the
// matrix. Also, row and column attributes are stored in the columnData and rowData
// structs. If the columns are un-sorted or contain duplicate row indexes,
// this routine will also sort and remove duplicate row indexes from the
// column form of the matrix.
//
// int columnIndexes [] array of the column indexes for the workArea, of size numberOfColumns+1
//
// returns true if the columns are jumbled
bool InitializeRowsAndColumns(int columnIndexes[])
{
int *cp; // a column pointer
int *cp_end; // a pointer to the end of a column
int *rp; // a row pointer
int *rp_end; // a pointer to the end of a row
int last_start; // start index of previous column in workArea
int start; // start index of column in workArea
bool jumbled_columns; // indicates if columns are jumbled
// Initialize the columnData
last_start = 0 ;
for (int i = 0 ; i < numberOfColumns ; i++)
{
start = columnIndexes[i];
// verify that the column indexes are ordered
if (start < last_start)
{
throw ArgumentException(Exception::Format("Column pointers must not be decreasing, last %d current %d",last_start,start));
}
columnData[i].start = start;
columnData[i].length = columnIndexes[i+1] - start;
columnData[i].thickness = 1;
columnData[i].score = 0;
columnData[i].prev = EMPTY;
columnData[i].degree_next = EMPTY;
last_start = start;
}
// verify the column index for the last column
if (columnIndexes[numberOfColumns] < last_start)
{
throw ArgumentException(Exception::Format("Column pointers must not be decreasing, last %d current %d",columnIndexes[numberOfColumns],last_start));
}
jumbled_columns = false;
for (int i = 0 ; i < numberOfRows ; i++)
{
rowData[i].length = 0;
rowData[i].mark = -1;
}
// Scan columns, compute row degrees, and check row indexes
for (int i = 0 ; i < numberOfColumns ; i++)
{
int previousRowIndex = -1 ;
cp = &workArea [columnIndexes [i]];
cp_end = &workArea [columnIndexes [i+1]];
while (cp < cp_end)
{
int rowIndex = *cp;
// make sure row indexes within range
if (rowIndex < 0 || rowIndex >= numberOfRows)
{
throw ArgumentException(Exception::Format("row index out of range, column %d row %d previous row %d",i, rowIndex, previousRowIndex));
}
else if (rowIndex <= previousRowIndex )
{
// row indexes are not sorted or repeated, thus cols are jumbled
jumbled_columns = true;
}
// prevent repeated row from being counted
if (rowData [rowIndex].mark != i)
{
rowData [rowIndex].length++ ;
rowData [rowIndex].mark = i;
previousRowIndex = rowIndex;
}
else
{
// this is a repeated entry in the column, it will be removed
columnData[i].length--;
}
cp++;
}
}
// Compute row pointers
// row form of the matrix starts directly after the column form of matrix in workArea
rowData [0].start = columnIndexes [numberOfColumns];
rowData [0].p = rowData [0].start ;
rowData [0].mark = -1 ;
for (int i = 1 ; i < numberOfRows ; i++)
{
rowData [i].start = rowData [i-1].start + rowData [i-1].length;
rowData [i].p = rowData [i].start;
rowData [i].mark = -1;
}
// Create row form
if (jumbled_columns)
{
// if cols jumbled, watch for repeated row indexes
for (int i = 0 ; i < numberOfColumns ; i++)
{
cp = &workArea [columnIndexes[i]];
cp_end = &workArea [columnIndexes[i+1]];
while (cp < cp_end)
{
int row = *cp++ ;
if (rowData [row].mark != i)
{
workArea [(rowData [row].p)++] = i;
rowData [row].mark = i;
}
}
}
}
else
{
// if cols not jumbled, we don't need the mark (this is faster)
for (int i = 0 ; i < numberOfColumns ; i++)
{
cp = &workArea [columnIndexes[i]];
cp_end = &workArea [columnIndexes[i+1]];
while (cp < cp_end)
{
workArea [(rowData [*cp++].p)++] = i;
}
}
}
// Clear the row marks and set row degrees
for (int i = 0 ; i < numberOfRows ; i++)
{
rowData [i].mark = 0 ;
rowData [i].degree = rowData [i].length;
}
// See if we need to re-create columns
if (jumbled_columns)
{
// Compute col pointers
// col form of the matrix starts at workArea [0].
// Note, we may have a gap between the col form and the row
// form if there were duplicate entries, if so, it will be
// removed upon the first garbage collection
columnData [0].start = 0 ;
columnIndexes[0] = columnData [0].start ;
for (int i = 1 ; i < numberOfColumns ; i++)
{
// note that the lengths here are for pruned columns, i.e.
// no duplicate row indexes will exist for these columns
columnData [i].start = columnData [i-1].start + columnData [i-1].length ;
columnIndexes[i] = columnData [i].start ;
}
// Re-create col form
for (int i = 0 ; i < numberOfRows; i++)
{
rp = &workArea [rowData [i].start];
rp_end = rp + rowData [i].length ;
while (rp < rp_end)
{
workArea [(columnIndexes[*rp++])++] = i;
}
}
return true;
}
else
{
// no columns jumbled (this is faster)
return false;
}
}
// InitializeScoring
//
//
// Kills dense or empty columns and rows, calculates an initial score for
// each column, and places all columns in the degree lists.
//
// Parameters:
//
// int head [], of size numberOfColumns+1
// int *p_n_row2, number of non-dense, non-empty rows
// int *p_n_col2, number of non-dense, non-empty columns
// int *p_max_deg maximum row degree
void InitializeScoring(int head [],int *p_n_row2,int *p_n_col2,int *p_max_deg)
{
int c; // a column index
int r, row; // a row index
int *cp; // a column pointer
int *cp_end; // a pointer to the end of a column
int *new_cp; // new column pointer
int col_length;// length of pruned column
int score; // current column score
int n_col2; // number of non-dense, non-empty columns
int n_row2; // number of non-dense, non-empty rows
int dense_row_count;// remove rows with more entries than this
int dense_col_count;// remove cols with more entries than this
int min_score; // smallest column score
int max_deg; // maximum row degree
int next_col; // Used to add to degree list.
#ifndef NDEBUG
int debug_count ; // debug only.
#endif
dense_row_count = std::max( 0, std::min( int( denseRow * numberOfColumns ), numberOfColumns ) );
dense_col_count = std::max( 0, std::min( int( denseCol * numberOfRows ), numberOfRows ) );
max_deg = 0;
n_col2 = numberOfColumns;
n_row2 = numberOfRows;
// Kill empty columns
// Put the empty columns at the end in their natural, so that LU
// factorization can proceed as far as possible.
for (c = numberOfColumns-1 ; c >= 0 ; c--)
{
if (columnData [c].length == 0)
{
// this is a empty column, kill and order it last
n_col2--;
columnData[c].order = n_col2;
columnData[c].KillPrincipal();
}
}
// Kill dense columns
// Put the dense columns at the end, in their natural order
for (c = numberOfColumns-1 ; c >= 0 ; c--)
{
// skip any dead columns
if (columnData[c].IsDead())
{
continue ;
}
int deg = columnData [c].length ;
if (deg > dense_col_count)
{
// this is a dense column, kill and order it last
columnData [c].order = --n_col2 ;
// decrement the row degrees
cp = &workArea [columnData [c].start] ;
cp_end = cp + columnData [c].length ;
while (cp < cp_end)
{
rowData [*cp++].degree--;
}
columnData[c].KillPrincipal();
}
}
Debug("Dense and null columns killed: %d\n", numberOfColumns - n_col2);
// Kill dense and empty rows
for (r = 0 ; r < numberOfRows ; r++)
{
int deg = rowData [r].degree ;
assert (deg >= 0 && deg <= numberOfColumns) ;
if (deg > dense_row_count || deg == 0)
{
// kill a dense or empty row
rowData [r].Kill();
--n_row2 ;
}
else
{
// keep track of max degree of remaining rows
max_deg = std::max( max_deg, deg );
}
}
Debug("Dense and null rows killed: %d\n", numberOfRows - n_row2);
// Compute initial column scores
// At this point the row degrees are accurate. They reflect the number
// of "live" (non-dense) columns in each row. No empty rows exist.
// Some "live" columns may contain only dead rows, however. These are
// pruned in the code below.
// now find the initial matlab score for each column
for (c = numberOfColumns-1 ; c >= 0 ; c--)
{
// skip dead column
if (columnData[c].IsDead())
{
continue ;
}
score = 0;
cp = &workArea [columnData [c].start];
new_cp = cp;
cp_end = cp + columnData [c].length;
while (cp < cp_end)
{
// get a row
row = *cp++ ;
// skip if dead
if (rowData[row].IsDead())
{
continue;
}
// compact the column
*new_cp++ = row ;
// add row's external degree
score += rowData [row].degree - 1 ;
// guard against integer overflow
score = std::min( score, numberOfColumns );
}
// determine pruned column length
col_length = (int) (new_cp - &workArea [columnData [c].start]) ;
if (col_length == 0)
{
// a newly-made null column (all rows in this col are "dense"
// and have already been killed)
Debug("Newly null killed: %d\n", c);
columnData [c].order = --n_col2 ;
columnData[c].KillPrincipal();
}
else
{
// set column length and set score
assert (score >= 0) ;
assert (score <= numberOfColumns) ;
columnData [c].length = col_length ;
columnData [c].score = score ;
}
}
Debug("Dense, null, and newly-null columns killed: %d\n",numberOfColumns-n_col2);
// At this point, all empty rows and columns are dead. All live columns
// are "clean" (containing no dead rows) and simplicial (no supercolumns
// yet). Rows may contain dead columns, but all live rows contain at
// least one live column.
// Initialize degree lists
#ifndef NDEBUG
debug_count = 0 ;
#endif
// clear the hash buckets
for (c = 0 ; c <= numberOfColumns ; c++)
{
head [c] = EMPTY;
}
min_score = numberOfColumns ;
// place in reverse order, so low column indexes are at the front
// of the lists. This will encourage natural tie-breaking
for (c = numberOfColumns-1 ; c >= 0 ; c--)
{
// only add principal columns to degree lists
if (columnData[c].IsAlive())
{
Debug("place %d score %d minscore %d ncol %d\n", c, columnData [c].score, min_score, numberOfColumns);
// Add columns score to DList
score = columnData [c].score ;
assert (min_score >= 0) ;
assert (min_score <= numberOfColumns) ;
assert (score >= 0) ;
assert (score <= numberOfColumns) ;
assert (head [score] >= EMPTY) ;
// now add this column to dList at proper score location
next_col = head [score] ;
columnData [c].prev = EMPTY ;
columnData [c].degree_next = next_col ;
// if there already was a column with the same score, set its
// previous pointer to this new column
if (next_col != EMPTY)
{
columnData [next_col].prev = c ;
}
head [score] = c ;
// see if this score is less than current min
min_score = std::min( min_score, score );
#ifndef NDEBUG
debug_count++ ;
#endif
}
}
// Return number of remaining columns, and max row degree
*p_n_col2 = n_col2 ;
*p_n_row2 = n_row2 ;
*p_max_deg = max_deg ;
}
// DetermineOrder
//
// Order the principal columns of the supercolumn form of the matrix
// (no supercolumns on input). Uses a minimum approximate column minimum
// degree ordering method. Not user-callable.
//
// Parameters
//
// int head [], of size numberOfColumns+1
// int n_col2, Remaining columns to order
// int max_deg, Maximum row degree
// int pfree index of first free slot (2*nnz on entry)
//
// returns the number of garbage collections
int DetermineOrder(int head [],int n_col2,int max_deg,int pfree)
{
int k; // current pivot ordering step
int pivot_col; // current pivot column
int *cp; // a column pointer
int *rp; // a row pointer
int pivot_row; // current pivot row
int *new_cp; // modified column pointer
int *new_rp; // modified row pointer
int pivot_row_start; // pointer to start of pivot row
int pivot_row_degree; // # of columns in pivot row
int pivot_row_length; // # of supercolumns in pivot row
int pivot_col_score; // score of pivot column
int needed_memory; // free space needed for pivot row
int *cp_end; // pointer to the end of a column
int *rp_end; // pointer to the end of a row
int row; // a row index
int col; // a column index
int max_score; // maximum possible score
int cur_score; // score of current column
int hash; // hash value for supernode detection
int head_column; // head of hash bucket
int first_col; // first column in hash bucket
int tag_mark; // marker value for mark array
int row_mark; // rowData [row].mark
int set_difference; // set difference size of row with pivot row
int min_score; // smallest column score
int col_thickness; // "thickness" (# of columns in a supercol)
int max_mark; // maximum value of tag_mark
int pivot_col_thickness; // number of columns represented by pivot col
int prev_col; // Used by Dlist operations.
int next_col; // Used by Dlist operations.
int ngarbage; // number of garbage collections performed
#ifndef NDEBUG
int debug_d ; // debug loop counter
int debug_step = 0 ; // debug loop counter
#endif
// Initialization and clear mark
max_mark = INT_MAX - numberOfColumns ; // INT_MAX defined in <limits.h>
tag_mark = ClearMark();
min_score = 0 ;
ngarbage = 0 ;
Debug("Ordering.. n_col2=%d\n", n_col2);
// Order the columns
for (k = 0 ; k < n_col2 ; )// 'k' is incremented below
{
// Select pivot column, and order it
// make sure degree list isn't empty
assert (min_score >= 0) ;
assert (min_score <= numberOfColumns) ;
assert (head [min_score] >= EMPTY) ;
#ifndef NDEBUG
for (debug_d = 0 ; debug_d < min_score ; debug_d++)
{
assert (head [debug_d] == EMPTY) ;
}
#endif
// get pivot column from head of minimum degree list
while (head [min_score] == EMPTY && min_score < numberOfColumns)
{
min_score++ ;
}
pivot_col = head [min_score] ;
assert (pivot_col >= 0 && pivot_col <= numberOfColumns) ;
next_col = columnData [pivot_col].degree_next ;
head [min_score] = next_col ;
if (next_col != EMPTY)
{
columnData [next_col].prev = EMPTY ;
}
assert (columnData[pivot_col].IsAlive()) ;
Debug("Pivot col: %d\n", pivot_col);
// remember score for defrag check
pivot_col_score = columnData [pivot_col].score ;
// the pivot column is the kth column in the pivot order
columnData [pivot_col].order = k ;
// increment order count by column thickness
pivot_col_thickness = columnData [pivot_col].thickness ;
k += pivot_col_thickness ;
assert (pivot_col_thickness > 0) ;
// Garbage_collection, if necessary
needed_memory = std::min( pivot_col_score, numberOfColumns - k );
if (pfree + needed_memory >= workAreaLength)
{
pfree = CompactWorkArea(&workArea [pfree]);
ngarbage++ ;
// after garbage collection we will have enough
assert (pfree + needed_memory < workAreaLength) ;
// garbage collection has wiped out the rowData[].mark array
tag_mark = ClearMark() ;
}
// Compute pivot row pattern
// get starting location for this new merged row
pivot_row_start = pfree ;
// initialize new row counts to zero
pivot_row_degree = 0 ;
// tag pivot column as having been visited so it isn't included
// in merged pivot row
columnData [pivot_col].thickness = -pivot_col_thickness ;
// pivot row is the union of all rows in the pivot column pattern
cp = &workArea [columnData [pivot_col].start] ;
cp_end = cp + columnData [pivot_col].length ;
while (cp < cp_end)
{
// get a row
row = *cp++ ;
Debug("Pivot col pattern %d %d\n", rowData[row].IsAlive(), row);
// skip if row is dead
if (rowData[row].IsDead())
{
continue ;
}
rp = &workArea [rowData [row].start] ;
rp_end = rp + rowData [row].length ;
while (rp < rp_end)
{
// get a column
col = *rp++ ;
// add the column, if alive and untagged
col_thickness = columnData [col].thickness ;
if (col_thickness > 0 && columnData[col].IsAlive() )
{
// tag column in pivot row
columnData [col].thickness = -col_thickness ;
assert (pfree < workAreaLength) ;
// place column in pivot row
workArea [pfree++] = col ;
pivot_row_degree += col_thickness ;
}
}
}
// clear tag on pivot column
columnData [pivot_col].thickness = pivot_col_thickness ;
max_deg = std::max( max_deg, pivot_row_degree );
// Kill all rows used to construct pivot row
// also kill pivot row, temporarily
cp = &workArea [columnData [pivot_col].start] ;
cp_end = cp + columnData [pivot_col].length ;
while (cp < cp_end)
{
// may be killing an already dead row
row = *cp++ ;
Debug("Kill row in pivot col: %d\n", row);
rowData [row].Kill();
}
// Select a row index to use as the new pivot row
pivot_row_length = pfree - pivot_row_start ;
if (pivot_row_length > 0)
{
// pick the "pivot" row arbitrarily (first row in col)
pivot_row = workArea [columnData [pivot_col].start] ;
Debug("Pivotal row is %d\n", pivot_row);
}
else
{
// there is no pivot row, since it is of zero length
pivot_row = EMPTY ;
assert (pivot_row_length == 0) ;
}
assert (columnData [pivot_col].length > 0 || pivot_row_length == 0) ;
// Approximate degree computation
// Here begins the computation of the approximate degree. The column
// score is the sum of the pivot row "length", plus the size of the
// set differences of each row in the column minus the pattern of the
// pivot row itself. The column ("thickness") itself is also
// excluded from the column score (we thus use an approximate
// external degree).
// The time taken by the following code (compute set differences, and
// add them up) is proportional to the size of the data structure
// being scanned - that is, the sum of the sizes of each column in
// the pivot row. Thus, the amortized time to compute a column score
// is proportional to the size of that column (where size, in this
// context, is the column "length", or the number of row indexes
// in that column). The number of row indexes in a column is
// monotonically non-decreasing, from the length of the original
// column on input to colamd.
// Compute set differences
Debug("** Computing set differences phase. **\n");
// pivot row is currently dead - it will be revived later.
Debug("Pivot row: ");
// for each column in pivot row
rp = &workArea [pivot_row_start] ;
rp_end = rp + pivot_row_length ;
while (rp < rp_end)
{
col = *rp++ ;
assert (columnData[col].IsAlive() && col != pivot_col) ;
Debug("columnData: %d\n", col);
// clear tags used to construct pivot row pattern
col_thickness = -columnData [col].thickness ;
assert (col_thickness > 0) ;
columnData [col].thickness = col_thickness ;
// Remove column from degree list
cur_score = columnData [col].score ;
prev_col = columnData [col].prev ;
next_col = columnData [col].degree_next ;
assert (cur_score >= 0) ;
assert (cur_score <= numberOfColumns) ;
assert (cur_score >= EMPTY) ;
if (prev_col == EMPTY)
{
head [cur_score] = next_col ;
}
else
{
columnData [prev_col].degree_next = next_col ;
}
if (next_col != EMPTY)
{
columnData [next_col].prev = prev_col ;
}
// Scan the column
cp = &workArea [columnData [col].start] ;
cp_end = cp + columnData [col].length ;
while (cp < cp_end)
{
/* get a row */
row = *cp++ ;
/* skip if dead */
if (rowData [row].IsDead())
{
continue ;
}
row_mark = rowData [row].mark ;
assert (row != pivot_row) ;
set_difference = row_mark - tag_mark ;
// check if the row has been seen yet
if (set_difference < 0)
{
assert (rowData [row].degree <= max_deg) ;
set_difference = rowData [row].degree ;
}
// subtract column thickness from this row's set difference
set_difference -= col_thickness ;
assert (set_difference >= 0) ;
// absorb this row if the set difference becomes zero
if (set_difference == 0)
{
Debug("aggressive absorption. rowData: %d\n", row);
rowData [row].Kill();
}
else
{
// save the new mark
rowData [row].mark = set_difference + tag_mark ;
}
}
}
// Add up set differences for each column
Debug("** Adding set differences phase. **\n");
// for each column in pivot row
rp = &workArea [pivot_row_start] ;
rp_end = rp + pivot_row_length ;
while (rp < rp_end)
{
// get a column
col = *rp++ ;
assert (columnData[col].IsAlive() && col != pivot_col) ;
hash = 0 ;
cur_score = 0 ;
cp = &workArea [columnData [col].start] ;
// compact the column
new_cp = cp ;
cp_end = cp + columnData [col].length ;
Debug("Adding set diffs for columnData: %d.\n", col);
while (cp < cp_end)
{
// get a row
row = *cp++ ;
assert(row >= 0 && row < numberOfRows) ;
// skip if dead
if (rowData [row].IsDead())
{
continue ;
}
row_mark = rowData [row].mark ;
assert (row_mark > tag_mark) ;
// compact the column
*new_cp++ = row ;
// compute hash function
hash += row ;
// add set difference
cur_score += row_mark - tag_mark ;
// integer overflow...
cur_score = std::min( cur_score, numberOfColumns );
}
// recompute the column's length
columnData [col].length = (int) (new_cp - &workArea [columnData [col].start]) ;
// Further mass elimination
if (columnData [col].length == 0)
{
Debug("further mass elimination. columnData: %d\n", col);
// nothing left but the pivot row in this column
columnData[col].KillPrincipal();
pivot_row_degree -= columnData [col].thickness;
assert (pivot_row_degree >= 0) ;
// order it
columnData [col].order = k ;
// increment order count by column thickness
k += columnData [col].thickness ;
}
else
{
// Prepare for supercolumn detection
Debug("Preparing supercol detection for columnData: %d.\n", col);
// save score so far
columnData [col].score = cur_score ;
// add column to hash table, for supercolumn detection
hash %= numberOfColumns + 1 ;
Debug(" Hash = %d, numberOfColumns = %d.\n", hash, numberOfColumns);
assert (hash <= numberOfColumns) ;
head_column = head [hash] ;
if (head_column > EMPTY)
{
// degree list "hash" is non-empty, use prev of
// first column in degree list as head of hash bucket
first_col = columnData [head_column].headhash ;
columnData [head_column].headhash = col ;
}
else
{
// degree list "hash" is empty, use head as hash bucket
first_col = - (head_column + 2) ;
head [hash] = - (col + 2) ;
}
columnData [col].hash_next = first_col ;
// save hash function in columnData [col].hash
columnData [col].hash = (int) hash ;
assert (columnData[col].IsAlive()) ;
}
}
// The approximate external column degree is now computed.
// Supercolumn detection
Debug("** Supercolumn detection phase. **\n");
detect_super_cols (head, pivot_row_start, pivot_row_length);
// Kill the pivotal column
columnData[pivot_col].KillPrincipal();
// Clear mark
tag_mark += (max_deg + 1) ;
if (tag_mark >= max_mark)
{
Debug("clearing tag_mark\n");
tag_mark = ClearMark() ;
}
// Finalize the new pivot row, and column scores
Debug("** Finalize scores phase. **\n");
// for each column in pivot row
rp = &workArea [pivot_row_start] ;
// compact the pivot row
new_rp = rp ;
rp_end = rp + pivot_row_length ;
while (rp < rp_end)
{
col = *rp++ ;
// skip dead columns
if (columnData[col].IsDead())
{
continue;
}
*new_rp++ = col ;
// add new pivot row to column
workArea [columnData [col].start + (columnData [col].length++)] = pivot_row;
// retrieve score so far and add on pivot row's degree.
// (we wait until here for this in case the pivot
// row's degree was reduced due to mass elimination).
cur_score = columnData [col].score + pivot_row_degree;
// calculate the max possible score as the number of
// external columns minus the 'k' value minus the
// columns thickness
max_score = numberOfColumns - k - columnData [col].thickness;
// make the score the external degree of the union-of-rows
cur_score -= columnData [col].thickness;
// make sure score is less or equal than the max score
cur_score = std::min( cur_score, max_score );
assert (cur_score >= 0);
// store updated score
columnData [col].score = cur_score;
// Place column back in degree list
assert (min_score >= 0);
assert (min_score <= numberOfColumns);
assert (cur_score >= 0);
assert (cur_score <= numberOfColumns);
assert (head [cur_score] >= EMPTY);
next_col = head [cur_score];
columnData [col].degree_next = next_col;
columnData [col].prev = EMPTY;
if (next_col != EMPTY)
{
columnData [next_col].prev = col;
}
head [cur_score] = col;
// see if this score is less than current min
min_score = std::min( min_score, cur_score );
}
// Resurrect the new pivot row
if (pivot_row_degree > 0)
{
// update pivot row length to reflect any cols that were killed
// during super-col detection and mass elimination
rowData [pivot_row].start = pivot_row_start ;
rowData [pivot_row].length = (int) (new_rp - &workArea[pivot_row_start]) ;
rowData [pivot_row].degree = pivot_row_degree ;
rowData [pivot_row].mark = 0 ;
// pivot row is no longer dead
}
}
// All principal columns have now been ordered
return ngarbage;
}
// order_children
//
// The find_ordering routine has ordered all of the principal columns (the
// representatives of the supercolumns). The non-principal columns have not
// yet been ordered. This routine orders those columns by walking up the
// parent tree (a column is a child of the column which absorbed it). The
// final permutation vector is then placed in p [0 ... numberOfColumns-1], with p [0]
// being the first column, and p [numberOfColumns-1] being the last. It doesn't look
// like it at first glance, but be assured that this routine takes time linear
// in the number of columns. Although not immediately obvious, the time
// taken by this routine is O (numberOfColumns), that is, linear in the number of
// columns.
//
// Parameters
//
// int numberOfColumns, number of columns of workArea
// ColumnData columnData [], of size numberOfColumns+1
// int p [] p [0 ... numberOfColumns-1] is the column permutation
void order_children(int p [])
{
int i; // loop counter for all columns
int c; // column index
int parent; // index of column's parent
int order; // column's order
// Order each non-principal column
for (i = 0 ; i < numberOfColumns ; i++)
{
// find an un-ordered non-principal column
assert (columnData[i].IsDead()) ;
if (!columnData[i].IsDeadPrincipal() && columnData [i].order == EMPTY)
{
parent = i ;
// once found, find its principal parent
do
{
parent = columnData [parent].parent ;
} while (!columnData[parent].IsDeadPrincipal()) ;
// now, order all un-ordered non-principal columns along path
// to this parent. collapse tree at the same time
c = i ;
// get order of parent
order = columnData [parent].order ;
do
{
assert (columnData [c].order == EMPTY) ;
// order this column
columnData [c].order = order++ ;
// collaps tree
columnData [c].parent = parent ;
// get immediate parent of this column
c = columnData [c].parent ;
// continue until we hit an ordered column. There are
// guarranteed not to be anymore unordered columns
// above an ordered column
} while (columnData [c].order == EMPTY) ;
// re-order the super_col parent to largest order for this group
columnData [parent].order = order ;
}
}
// Generate the permutation
for (c = 0 ; c < numberOfColumns ; c++)
{
p [columnData [c].order] = c ;
}
}
// detect_super_cols
//
// Detects supercolumns by finding matches between columns in the hash buckets.
// Check amongst columns in the set workArea [row_start ... row_start + row_length-1].
// The columns under consideration are currently *not* in the degree lists,
// and have already been placed in the hash buckets.
//
// The hash bucket for columns whose hash function is equal to h is stored
// as follows:
//
// if head [h] is >= 0, then head [h] contains a degree list, so:
//
// head [h] is the first column in degree bucket h.
// columnData [head [h]].headhash gives the first column in hash bucket h.
//
// otherwise, the degree list is empty, and:
//
// -(head [h] + 2) is the first column in hash bucket h.
//
// For a column c in a hash bucket, columnData [c].prev is NOT a "previous
// column" pointer. columnData [c].hash is used instead as the hash number
// for that column. The value of columnData [c].hash_next is the next column
// in the same hash bucket.
//
// Assuming no, or "few" hash collisions, the time taken by this routine is
// linear in the sum of the sizes (lengths) of each column whose score has
// just been computed in the approximate degree computation.
//
// Parameters
//
// int head [] head of degree lists and hash buckets
// int row_start pointer to set of columns to check
// int row_length number of columns to check
//
void detect_super_cols(int head [],int row_start,int row_length)
{
int hash; // hash # for a column
int *rp; // pointer to a row
int c; // a column index
int super_c; // column index of the column to absorb into
int *cp1; // column pointer for column super_c
int *cp2; // column pointer for column c
int length; // length of column super_c
int prev_c; // column preceding c in hash bucket
int i; // loop counter
int *rp_end; // pointer to the end of the row
int col; // a column index in the row to check
int head_column;// first column in hash bucket or degree list
int first_col; // first column in hash bucket
// Consider each column in the row
rp = &workArea [row_start];
rp_end = rp + row_length;
while (rp < rp_end)
{
col = *rp++;
if (columnData[col].IsDead())
{
continue;
}
// get hash number for this column
hash = columnData[col].hash;
assert (hash <= numberOfColumns);
// Get the first column in this hash bucket
head_column = head [hash];
if (head_column > EMPTY)
{
first_col = columnData [head_column].headhash;
}
else
{
first_col = - (head_column + 2);
}
// Consider each column in the hash bucket
for (super_c = first_col ; super_c != EMPTY; super_c = columnData [super_c].hash_next)
{
assert (columnData[super_c].IsAlive() ) ;
assert (columnData [super_c].hash == hash) ;
length = columnData [super_c].length ;
// prev_c is the column preceding column c in the hash bucket
prev_c = super_c ;
// Compare super_c with all columns after it
for (c = columnData [super_c].hash_next ;c != EMPTY ; c = columnData [c].hash_next)
{
assert (c != super_c);
assert (columnData[c].IsAlive());
assert (columnData [c].hash == hash);
// not identical if lengths or scores are different
if (columnData [c].length != length ||columnData [c].score != columnData [super_c].score)
{
prev_c = c;
continue;
}
// compare the two columns
cp1 = &workArea [columnData [super_c].start];
cp2 = &workArea [columnData [c].start];
for (i = 0 ; i < length ; i++)
{
// the columns are "clean" (no dead rows)
assert (rowData[(*cp1)].IsAlive() ) ;
assert (rowData[(*cp2)].IsAlive() ) ;
// row indexes will same order for both supercols,
// no gather scatter nessasary
if (*cp1++ != *cp2++)
{
break ;
}
}
// the two columns are different if the for-loop "broke"
if (i != length)
{
prev_c = c ;
continue ;
}
// Got it! two columns are identical
assert (columnData [c].score == columnData [super_c].score) ;
columnData [super_c].thickness += columnData [c].thickness ;
columnData [c].parent = super_c ;
columnData[c].KillNonPrincipal();
// order c later, in order_children()
columnData [c].order = EMPTY ;
// remove c from hash bucket
columnData [prev_c].hash_next = columnData [c].hash_next ;
}
}
// Empty this hash bucket
if (head_column > EMPTY)
{
// corresponding degree list "hash" is not empty
columnData [head_column].headhash = EMPTY ;
}
else
{
// corresponding degree list "hash" is empty
head [hash] = EMPTY ;
}
}
}
// CompactWorkArea
//
// Defragments and compacts columns and rows in the workspace workArea. Used when
// all avaliable memory has been used while performing row merging. Returns
// the index of the first free position in workArea, after garbage collection. The
// time taken by this routine is linear is the size of the array workArea, which is
// itself linear in the number of nonzeros in the input matrix.
//
// Parameters
//
// int *pfree &workArea [0] ... pfree is in use
//
// returns the new value of pfree
//
int CompactWorkArea(int *pfree)
{
int *psrc; // source pointer
int *pdest; // destination pointer
int j; // counter
int r; // a row index
int c; // a column index
int length; // length of a row or column
#ifndef NDEBUG
int debug_rows ;
Debug("Defrag..\n");
for (psrc = &workArea[0] ; psrc < pfree ; psrc++) assert (*psrc >= 0) ;
debug_rows = 0 ;
#endif
// Defragment the columns
pdest = &workArea[0];
for (c = 0 ; c < numberOfColumns ; c++)
{
if (columnData[c].IsAlive())
{
psrc = &workArea[columnData[c].start];
// move and compact the column
assert (pdest <= psrc) ;
columnData[c].start = (int) (pdest - &workArea[0]);
length = columnData[c].length;
for (j = 0 ; j < length ; j++)
{
r = *psrc++ ;
if (rowData[r].IsAlive() )
{
*pdest++ = r ;
}
}
columnData[c].length = (int)(pdest - &workArea [columnData[c].start]);
}
}
// Prepare to defragment the rows
for (r = 0 ; r < numberOfRows ; r++)
{
if (rowData [r].IsAlive())
{
if (rowData [r].length == 0)
{
// this row is of zero length. cannot compact it, so kill it
Debug("Defrag row kill\n");
rowData [r].Kill();
}
else
{
// save first column index in rowData [r].first_column
psrc = &workArea [rowData [r].start];
rowData [r].first_column = *psrc;
assert (rowData[r].IsAlive());
// flag the start of the row
*psrc |= FLAG;
#ifndef NDEBUG
debug_rows++;
#endif
}
}
}
// Defragment the rows
psrc = pdest ;
while (psrc < pfree)
{
// find a negative number ... the start of a row
if (*psrc & FLAG)
{
// get back the row index by removing the flag bit
r = *psrc & MASK;
assert (r >= 0 && r < numberOfRows) ;
// restore first column index
*psrc = rowData [r].first_column ;
assert (rowData[r].IsAlive()) ;
// move and compact the row
assert (pdest <= psrc) ;
rowData [r].start = (int) (pdest - &workArea [0]) ;
length = rowData [r].length ;
for (j = 0 ; j < length ; j++)
{
c = *psrc++ ;
if ( columnData[c].IsAlive() )
{
*pdest++ = c ;
}
}
rowData [r].length = (int) (pdest - &workArea [rowData [r].start]) ;
#ifndef NDEBUG
debug_rows-- ;
#endif
}
psrc++;
}
// ensure we found all the rows
assert (debug_rows == 0) ;
// Return the new value of pfree
return ((int) (pdest - &workArea [0])) ;
}
// ClearMark
//
// Clears the rowData [].mark array, and returns the new tag_mark.
// Return value is the new tag_mark
//
// Parameters
//
// int numberOfRows number of rows in workArea
// RowData rowData [] rowData [0 ... numberOfRows-1].mark is set to zero
//
// return the new value for tag_mark
int ClearMark()
{
int r;
Debug("Clear mark\n");
for (r = 0 ; r < numberOfRows ; r++)
{
if (rowData[r].IsAlive())
{
rowData [r].mark = 0 ;
}
}
return 1;
}
public:
static size_t Recommended(size_t nnz, size_t numberOfRows, size_t numberOfColumns)
{
size_t minimum; // bare minimum requirements
size_t recommended; // recommended value of Alen
minimum =
2 * (nnz) // for workArea
+ (((numberOfColumns) + 1) * sizeof (ColumnData) / sizeof (int)) // for columnData
+ (((numberOfRows) + 1) * sizeof (RowData) / sizeof (int)) // for rowData
+ numberOfColumns // minimum elbow room to guarrantee success
+ STATS ; // for output statistics
/* recommended is equal to the minumum plus enough memory to keep the */
/* number garbage collections low */
recommended = minimum + nnz/5;
return recommended;
}
// Parameters
//
// int theNumberOfRows number of rows in theMatrixData
// int theNumberOfColumns number of columns in theMatrixData
// int theMatrixDataLength length of theMatrixData
// int theMatrixData[] row indexes of theMatrixData
// int p [] column indexes in theMatrixData
//
// returns true if successful
bool Execute(int theNumberOfRows,int theNumberOfColumns,int theMatrixDataLength,int theMatrixData [],int p [])
{
int nnz; // nonzeros in theMatrixData
int Row_size; // size of rowData [], in integers
int Col_size; // size of columnData [], in integers
int elbow_room; // remaining free space
int n_col2; // number of non-dense, non-empty columns
int n_row2; // number of non-dense, non-empty rows
int ngarbage; // number of garbage collections performed
int max_deg; // maximum row degree
int init_result; // return code from initialization
#if 0
#ifndef NDEBUG
debug_colamd = 0; // no debug printing
// get "D" environment variable, which gives the debug printing level
if (getenv ("D")) debug_colamd = atoi (getenv ("D")) ;
Debug("debug version, D = %d (THIS WILL BE SLOOOOW!)\n", debug_colamd);
#endif
#endif
// Check the input arguments
if (!theNumberOfRows || !theNumberOfColumns || !theMatrixDataLength || !theMatrixData || !p)
{
// numberOfRows and numberOfColumns must be non-negative, workArea and p must be present
Debug("error! %d %d %d\n", theNumberOfRows, theNumberOfColumns, theMatrixDataLength);
return (FALSE) ;
}
this->workArea = theMatrixData;
numberOfRows = theNumberOfRows;
numberOfColumns = theNumberOfColumns;
workAreaLength = theMatrixDataLength;
nnz = p [numberOfColumns] ;
if (nnz < 0 || p [0] != 0)
{
// nnz must be non-negative, and p [0] must be zero
Debug("colamd error! %d %d\n", nnz, p [0]);
return false;
}
// Allocate the rowData and columnData arrays from array workArea
Col_size = (numberOfColumns + 1) * sizeof (ColumnData) / sizeof (int);
Row_size = (numberOfRows + 1) * sizeof (RowData) / sizeof (int);
elbow_room = theMatrixDataLength - (2*nnz + Col_size + Row_size);
if (elbow_room < numberOfColumns + STATS)
{
// not enough space in array workArea to perform the ordering
Debug("error! elbow_room %d, %d\n", elbow_room,numberOfColumns);
return (FALSE) ;
}
workAreaLength = 2*nnz + elbow_room;
columnData = (ColumnData *) &workArea[workAreaLength];
rowData = (RowData *) &workArea[workAreaLength + Col_size];
// Construct the row and column data structures
init_result = InitializeRowsAndColumns(p);
if (init_result == -1)
{
// input matrix is invalid
Debug("error! matrix invalid\n");
return (FALSE) ;
}
// Initialize scores, kill dense rows/columns
InitializeScoring(p, &n_row2, &n_col2, &max_deg) ;
// Order the supercolumns
ngarbage = DetermineOrder(p,n_col2, max_deg, 2*nnz);
// Order the non-principal columns
order_children (p);
// Return statistics in theMatrixData
for (int i = 0 ; i < STATS ; i++)
{
theMatrixData [i] = 0 ;
}
theMatrixData [DENSE_ROW] = numberOfRows - n_row2 ;
theMatrixData [DENSE_COL] = numberOfColumns - n_col2 ;
theMatrixData [DEFRAG_COUNT] = ngarbage ;
theMatrixData [JUMBLED_COLS] = init_result ;
return true;
}
};
};
};
#endif //__HNUMCOLAMD_H__
Submit an article or tip