Commit 8954a322 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Experimental greedy hypergraph partitioner in inner solver

parent 91921dc1
......@@ -106,10 +106,116 @@ namespace INMOST
{
while( size ) keys[PopHeap()] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
}
bool Contains(INMOST_DATA_ENUM_TYPE i)
bool Contains(INMOST_DATA_ENUM_TYPE i) const
{
return index[i] != ENUMUNDEF;
}
INMOST_DATA_ENUM_TYPE Size() const {return size;}
bool Empty() const {return size == 0;}
};
template<typename KeyType, typename Comparator>
class BinaryHeapCustom
{
INMOST_DATA_ENUM_TYPE size_max, size;
array<INMOST_DATA_ENUM_TYPE> heap;
array<INMOST_DATA_ENUM_TYPE> index;
array<KeyType> keys;
void swap(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j)
{
INMOST_DATA_ENUM_TYPE t = heap[i];
heap[i] = heap[j];
heap[j] = t;
index[heap[i]] = i;
index[heap[j]] = j;
}
void bubbleUp(INMOST_DATA_ENUM_TYPE k)
{
while(k > 1 && !Comparator()(keys[heap[k/2]],keys[heap[k]]))
{
swap(k, k/2);
k = k/2;
}
}
void bubbleDown(INMOST_DATA_ENUM_TYPE k)
{
INMOST_DATA_ENUM_TYPE j;
while(2*k <= size)
{
j = 2*k;
if(j < size && !Comparator()(keys[heap[j]],keys[heap[j+1]]) )
j++;
if(Comparator()(keys[heap[k]],keys[heap[j]]) )
break;
swap(k, j);
k = j;
}
}
public:
BinaryHeapCustom(INMOST_DATA_ENUM_TYPE len)
{
size_max = len;
size = 0;
keys.resize(size_max);
heap.resize(size_max+1);
index.resize(size_max+1);
for(INMOST_DATA_ENUM_TYPE i = 0; i <= size_max; i++)
index[i] = ENUMUNDEF;
}
void PushHeap(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
size++;
index[i] = size;
heap[size] = i;
keys[i] = key;
bubbleUp(size);
}
INMOST_DATA_ENUM_TYPE PopHeap()
{
if( size == 0 ) return ENUMUNDEF;
INMOST_DATA_ENUM_TYPE min = heap[1];
swap(1, size--);
bubbleDown(1);
index[min] = ENUMUNDEF;
heap[size+1] = ENUMUNDEF;
return min;
}
INMOST_DATA_ENUM_TYPE PeekHeap()
{
if( size == 0 ) return ENUMUNDEF;
return heap[1];
}
void DecreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
keys[i] = key;
bubbleUp(index[i]);
}
void IncreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
keys[i] = key;
bubbleDown(index[i]);
}
void ChangeKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
keys[i] = key;
if( !Comparator()(keys[i],key) )
bubbleUp(index[i]);
else
bubbleDown(index[i]);
}
KeyType & GetKey(INMOST_DATA_ENUM_TYPE i) {return keys[i];}
const KeyType & GetKey(INMOST_DATA_ENUM_TYPE i) const {return keys[i];}
void Clear()
{
while( size ) PopHeap();
}
bool Contains(INMOST_DATA_ENUM_TYPE i) const
{
return index[i] != ENUMUNDEF;
}
INMOST_DATA_ENUM_TYPE Size() const {return size;}
bool Empty() const {return size == 0;}
};
//this structure is used to provide sorting routing for Reverse-Cuthill-McKee reordering
......@@ -124,7 +230,36 @@ namespace INMOST
RCM_Comparator & operator = (RCM_Comparator const & b) {wbeg = b.wbeg; xadj = b.xadj; return *this;}
bool operator () (INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j)
{return xadj[i+1-wbeg] -xadj[i-wbeg] < xadj[j+1-wbeg] - xadj[j-wbeg];}
bool operator () (const std::pair<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> & i, const std::pair<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> & j)
{
INMOST_DATA_ENUM_TYPE ci = xadj[i.first+1-wbeg] - xadj[i.first-wbeg];
INMOST_DATA_ENUM_TYPE cj = xadj[j.first+1-wbeg] - xadj[j.first-wbeg];
return (ci < cj) || (ci == cj && i.second < i.first);
}
};
struct WRCM_Comparator
{
INMOST_DATA_ENUM_TYPE wbeg;
std::vector<INMOST_DATA_REAL_TYPE> & wadj;
public:
WRCM_Comparator(INMOST_DATA_ENUM_TYPE _wbeg, std::vector<INMOST_DATA_REAL_TYPE> & _wadj)
: wbeg(_wbeg), wadj(_wadj) {}
WRCM_Comparator(const WRCM_Comparator & b) : wbeg(b.wbeg), wadj(b.wadj) {}
WRCM_Comparator & operator = (WRCM_Comparator const & b) {wbeg = b.wbeg; wadj = b.wadj; return *this;}
bool operator () (INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j)
{return wadj[i-wbeg] < wadj[j-wbeg];}
bool operator () (const std::pair<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> & i, const std::pair<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> & j)
{
INMOST_DATA_REAL_TYPE ci = wadj[i.first-wbeg];
INMOST_DATA_REAL_TYPE cj = wadj[j.first-wbeg];
return (ci < cj) || (ci == cj && i.second < i.first);
}
};
}
#endif //INMOST_UTILS_H
......@@ -6,15 +6,78 @@
#include "inmost_solver.h"
#include "../solver_prototypes.hpp"
/*
class Range
{
INMOST_DATA_ENUM_TYPE first, last;
public:
Range(const Range & b):first(b.first), last(b.last) {}
Range(INMOST_DATA_ENUM_TYPE first, INMOST_DATA_ENUM_TYPE last) : first(first), last(last) {}
Range & operator =(Range const & b) {first = b.first; last = b.last; return *this;}
~Range() {}
INMOST_DATA_ENUM_TYPE Size() { return last - first; }
INMOST_DATA_ENUM_TYPE First() {return first;}
INMOST_DATA_ENUM_TYPE Last() {return last;}
};
class CSR_Row
{
Range & range;
std::vector<Sparse::Row::entry> & entries;
public:
}
class CSR_Matrix
{
interval<INMOST_DATA_ENUM_TYPE, Range> Address;
std::vector<Sparse::Row::entry> Entries;
public:
};
*/
class MLMTILUC_preconditioner : public Method
{
typedef struct Block_t
{
INMOST_DATA_ENUM_TYPE row_start, row_end;
INMOST_DATA_ENUM_TYPE col_start, col_end;
bool separator;
Block_t(INMOST_DATA_ENUM_TYPE rows,
INMOST_DATA_ENUM_TYPE rowe,
INMOST_DATA_ENUM_TYPE cols,
INMOST_DATA_ENUM_TYPE cole,
bool separator)
: row_start(rows), row_end(rowe),
col_start(cols), col_end(cole),
separator(separator) {}
Block_t(INMOST_DATA_ENUM_TYPE rows,
INMOST_DATA_ENUM_TYPE rowe,
INMOST_DATA_ENUM_TYPE cols,
INMOST_DATA_ENUM_TYPE cole)
: row_start(rows), row_end(rowe),
col_start(cols), col_end(cole),
separator(false) {}
INMOST_DATA_ENUM_TYPE RowSize() const {return row_end - row_start;}
INMOST_DATA_ENUM_TYPE ColSize() const {return col_end - col_start;}
// Block_t(const Block_t & b)
// : row_start(b.row_start), row_end(b.row_end),
// col_start(b.col_start), col_end(b.col_end) {}
// Block_t & operator =(Block_t const & b)
// {
// row_start = b.row_start;
// row_end = b.row_end;
// col_start = b.col_start;
// col_end = b.col_end;
// return *this;
// }
} Block;
typedef struct Interval_t
{
INMOST_DATA_ENUM_TYPE first, last;
INMOST_DATA_ENUM_TYPE Size() { return last - first; }
INMOST_DATA_ENUM_TYPE Size() const { return last - first; }
} Interval;
typedef dynarray<INMOST_DATA_ENUM_TYPE,256> levels_t;
std::vector<Sparse::Row::entry> LU_Entries, B_Entries;
std::vector<Sparse::Row::entry> L_Entries, U_Entries, B_Entries;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> LU_Diag;
interval<INMOST_DATA_ENUM_TYPE, Interval> U_Address, L_Address, B_Address;
std::vector<interval<INMOST_DATA_ENUM_TYPE, Interval> *> F_Address, E_Address;
......@@ -59,6 +122,139 @@ class MLMTILUC_preconditioner : public Method
interval<INMOST_DATA_ENUM_TYPE, bool> & donePQ,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ);
INMOST_DATA_REAL_TYPE AddListOrdered(INMOST_DATA_ENUM_TYPE cbeg,
Interval & Address,
std::vector<Sparse::Row::entry> & Entries,
INMOST_DATA_REAL_TYPE coef,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues,
INMOST_DATA_REAL_TYPE droptol);
INMOST_DATA_REAL_TYPE AddListUnordered(INMOST_DATA_ENUM_TYPE & Sbeg,
Interval & Address,
std::vector<Sparse::Row::entry> & Entries,
INMOST_DATA_REAL_TYPE coef,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues,
INMOST_DATA_REAL_TYPE droptol);
void OrderList(INMOST_DATA_ENUM_TYPE & Sbeg,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
std::vector<INMOST_DATA_ENUM_TYPE> & indices);
void ScaleList(INMOST_DATA_REAL_TYPE coef,
INMOST_DATA_ENUM_TYPE Sbeg,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues);
INMOST_DATA_REAL_TYPE Estimator1(INMOST_DATA_ENUM_TYPE k,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & Est,
INMOST_DATA_REAL_TYPE & mu_update);
INMOST_DATA_REAL_TYPE Estimator2(INMOST_DATA_ENUM_TYPE k,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & Est,
INMOST_DATA_REAL_TYPE & mu_update);
void EstimatorUpdate(INMOST_DATA_ENUM_TYPE k,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & Est,
INMOST_DATA_REAL_TYPE & mu_update);
void DiagonalUpdate(INMOST_DATA_ENUM_TYPE k,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & Diag,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndecesL,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValuesL,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeceU,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValuesU);
void ClearList(INMOST_DATA_ENUM_TYPE beg,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & LineIndeces,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & LineValues);
void PrepareGraph(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
const interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
const std::vector<Sparse::Row::entry> & Entries,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G);
// tG should be preallocated by number of columns, that may be wider then wbeg:wend
void PrepareGraphTranspose(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
const interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & tG);
// pG should be preallocated by wbeg:wend, tind computed with PrepareGraphTranspose
void PrepareGraphProduct(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
const interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G,
const interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & tG,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & pG);
// moves wbeg:wend subset of rows from tG_in to tG_out and filters entries outside wbeg:wend
// tG_in is not preserved to save memory
void FilterGraphProduct(const Block & b,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & pG_in,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & pG_out);
void FilterGraphTranspose(const Block & b,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invQ,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & tG_in,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & tG_out);
void FilterGraph(const Block & b,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & invP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G_in,
interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G_out);
// finds permutation that separates matrix into blocks by running GreedyDissection recursively
void NestedDissection(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
const interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
const std::vector<Sparse::Row::entry> & Entries,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
std::vector<Block> & blocks,
INMOST_DATA_ENUM_TYPE max_size);
void KwayDissection(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
const interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
const std::vector<Sparse::Row::entry> & Entries,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
std::vector<Block> & blocks,
int parts);
// finds permutation that separates matrix into blocks
void GreedyDissection(const Block & b,
const interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G,
const interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & tG,
const interval< INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & pG,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
std::vector<Block> & blocks, int kway_parts);
// int wgt_sep, int wgt_blk, int kway_parts);
//compute interval of column indices
void ColumnInterval(INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
const interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
const std::vector<Sparse::Row::entry> & Entries,
INMOST_DATA_ENUM_TYPE & cbeg,
INMOST_DATA_ENUM_TYPE & cend);
// column permutations
// has to call functions to permute scaling and E,F blocks
void ReorderMatrixQ(interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
std::vector<Sparse::Row::entry> & Entries,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ);
// rows and columns permutations
// has to call functions to permute scaling and E,F blocks
void ReorderMatrixPQ(interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
std::vector<Sparse::Row::entry> & Entries,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP);
void PrepareTranspose(INMOST_DATA_ENUM_TYPE cbeg,
INMOST_DATA_ENUM_TYPE cend,
interval<INMOST_DATA_ENUM_TYPE, Interval> & Address,
std::vector<Sparse::Row::entry> & Entries,
interval<INMOST_DATA_ENUM_TYPE, std::vector< std::pair<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> > > & Indices );
void DumpGraph(std::string name, interval<INMOST_DATA_ENUM_TYPE, std::vector<INMOST_DATA_ENUM_TYPE> > & G);
int Thread();
int Threads();
public:
INMOST_DATA_ENUM_TYPE & EnumParameter(std::string name);
INMOST_DATA_REAL_TYPE & RealParameter(std::string name);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment