Skip to content

GitLab

Commit 84d9eeda authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Updates and fixes in global solvers 

Fix binary heap with key update in maximum product transversal
algorithm in inner solvers.

Correct global rescaling strategy in MLMPTILUC2 inner solver.
parent 5be271af
Expand all Hide whitespace changes
Inline Side-by-side
Source/Headers/inmost_dense.h
View file @ 84d9eeda
......@@ -2422,111 +2422,84 @@ namespace INMOST
class BinaryHeapDense
{
INMOST_DATA_ENUM_TYPE size_max, size;
std::vector<INMOST_DATA_ENUM_TYPE> heap;
std::vector<INMOST_DATA_ENUM_TYPE> index;
INMOST_DATA_REAL_TYPE * keys;
INMOST_DATA_REAL_TYPE * Base;
std::vector<INMOST_DATA_REAL_TYPE *> Array;
std::vector<INMOST_DATA_ENUM_TYPE> Position;
public:
void Clear()
{
while(!Array.empty())
{
Position[Array.back()-Base] = ENUMUNDEF;
Array.pop_back();
}
}
INMOST_DATA_REAL_TYPE * Get(INMOST_DATA_ENUM_TYPE pos) {return Array[pos];}
INMOST_DATA_ENUM_TYPE GetSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(Array.size());}
INMOST_DATA_ENUM_TYPE GetPosition(INMOST_DATA_ENUM_TYPE pos)
void swap(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j)
{
return Position[pos];
INMOST_DATA_ENUM_TYPE t = heap[i];
heap[i] = heap[j];
heap[j] = t;
index[heap[i]] = i;
index[heap[j]] = j;
}
INMOST_DATA_ENUM_TYPE DecreaseKey(INMOST_DATA_ENUM_TYPE pos)
void bubbleUp(INMOST_DATA_ENUM_TYPE k)
{
INMOST_DATA_ENUM_TYPE i = Position[pos];
++i;
while(i > 1)
while(k > 1 && keys[heap[k/2]] > keys[heap[k]])
{
//if((*Array[i-1]) <= (*Array[i/2-1]))
if( compare(Array[i-1],Array[i/2-1]) )
{
Position[(Array[i/2-1]-Base)] = i-1;
Position[(Array[i-1]-Base)] = i/2-1;
std::swap(Array[i/2-1],Array[i-1]);
}
else break;
i = i/2;
swap(k, k/2);
k = k/2;
}
return i;
}
INMOST_DATA_ENUM_TYPE PushHeap(INMOST_DATA_REAL_TYPE * key)
void bubbleDown(INMOST_DATA_ENUM_TYPE k)
{
INMOST_DATA_ENUM_TYPE i = GetSize();
Array.push_back(key);
Position[(key-Base)] = i;
++i;
while(i > 1)
INMOST_DATA_ENUM_TYPE j;
while(2*k <= size)
{
//if((*Array[i-1]) <= (*Array[i/2-1]) )
if( compare(Array[i-1],Array[i/2-1]) )
{
Position[(Array[i-1]-Base)] = i/2-1;
Position[(Array[i/2-1]-Base)] = i-1;
std::swap(Array[i-1],Array[i/2-1]);
}
else break;
i = i/2;
j = 2*k;
if(j < size && keys[heap[j]] > keys[heap[j+1]])
j++;
if(keys[heap[k]] <= keys[heap[j]])
break;
swap(k, j);
k = j;
}
return i;
}
void BalanceHeap(INMOST_DATA_ENUM_TYPE i)
public:
BinaryHeapDense(INMOST_DATA_REAL_TYPE * pkeys, INMOST_DATA_ENUM_TYPE len)
{
INMOST_DATA_ENUM_TYPE Index;
++i;
while(i <= Array.size()/2)
{
if( 2*i+1 > Array.size() )
Index = 2*i;
//else if( (*Array[2*i-1]) <= (*Array[2*i+1-1]) )
else if( compare(Array[2*i-1],Array[2*i+1-1]) )
Index = 2*i;
else
Index = 2*i+1;
//if(!((*Array[i-1]) <= (*Array[Index-1])))
if(!compare(Array[i-1],Array[Index-1]))
{
Position[(Array[i-1]-Base)] = Index-1;
Position[(Array[Index-1]-Base)] = i-1;
std::swap(Array[i-1],Array[Index-1]);
}
else break;
i = Index;
}
size_max = len;
keys = pkeys;
size = 0;
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()
{
INMOST_DATA_ENUM_TYPE Ret = ENUMUNDEF;
if(Array.empty()) return Ret;
Ret = static_cast<INMOST_DATA_ENUM_TYPE>(Array[0]-Base);
Array[0] = Array.back();
Position[Array[0] - Base] = 0;
Array.pop_back();
Position[Ret] = ENUMUNDEF;
BalanceHeap(0);
return Ret;
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;
}
void DecreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
keys[i] = key;
bubbleUp(index[i]);
}
BinaryHeapDense(INMOST_DATA_REAL_TYPE * Base, INMOST_DATA_ENUM_TYPE Size)
: Base(Base)
void Clear()
{
Position.resize(Size,ENUMUNDEF);
Array.reserve(4096);
while( size ) keys[PopHeap()] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
}
~BinaryHeapDense()
bool Contains(INMOST_DATA_ENUM_TYPE i)
{
return index[i] != ENUMUNDEF;
}
};
template<typename Var>
void AbstractMatrix<Var>::MPT(INMOST_DATA_ENUM_TYPE * Perm, INMOST_DATA_REAL_TYPE * SL, INMOST_DATA_REAL_TYPE * SR) const
......@@ -2655,12 +2628,11 @@ namespace INMOST
}
else if( l < Dist[Lit] )
{
Dist[Lit] = l;
Parent[Perm[Lit]] = Li;
if( Heap.GetPosition(Lit) != ENUMUNDEF )
Heap.DecreaseKey(Lit);
if( Heap.Contains(Lit) )
Heap.DecreaseKey(Lit,l);
else
Heap.PushHeap(&Dist[Lit]);
Heap.PushHeap(Lit,l);
}
}
}
......@@ -2711,7 +2683,6 @@ namespace INMOST
Trace = Parent[Trace];
}
for(Ui = 0; Ui < Heap.GetSize(); ++Ui) *Heap.Get(Ui) = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
Heap.Clear();
}
}
......
Source/Misc/utils.cpp
View file @ 84d9eeda
......@@ -4,26 +4,29 @@
namespace INMOST {
template<>
int from_string<int>(std::string str) {
int from_string<int>(std::string str)
{
return atoi(str.c_str());
}
template<>
double from_string<double>(std::string str) {
double from_string<double>(std::string str)
{
return atof(str.c_str());
}
template<>
std::string from_string<std::string>(std::string str) {
std::string from_string<std::string>(std::string str)
{
return str;
}
std::string string_to_lower(const std::string &str) {
std::string string_to_lower(const std::string &str)
{
std::string lower = std::string(str);
for (int i = 0; i < lower.length(); i++) {
for (int i = 0; i < lower.length(); i++)
lower[i] = (char) tolower(lower[i]);
}
return lower;
}
}
\ No newline at end of file
}
Source/Misc/utils.h
View file @ 84d9eeda
......@@ -3,11 +3,14 @@
#include <sstream>
#include <string>
#include "inmost_common.h"
namespace INMOST {
namespace INMOST
{
template<typename T>
T from_string(std::string str) {
T from_string(std::string str)
{
T v;
std::istringstream ss(str);
ss >> v;
......@@ -15,14 +18,112 @@ namespace INMOST {
}
template<typename T>
std::string to_string(T value) {
std::string to_string(T value)
{
std::stringstream ss;
ss << value;
return ss.str();
}
std::string string_to_lower(const std::string &str);
class BinaryHeap
{
INMOST_DATA_ENUM_TYPE size_max, size;
std::vector<INMOST_DATA_ENUM_TYPE> heap;
std::vector<INMOST_DATA_ENUM_TYPE> index;
INMOST_DATA_REAL_TYPE * 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 && 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 && keys[heap[j]] > keys[heap[j+1]])
j++;
if(keys[heap[k]] <= keys[heap[j]])
break;
swap(k, j);
k = j;
}
}
public:
BinaryHeap(INMOST_DATA_REAL_TYPE * pkeys, INMOST_DATA_ENUM_TYPE len)
{
size_max = len;
keys = pkeys;
size = 0;
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;
}
void DecreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
{
keys[i] = key;
bubbleUp(index[i]);
}
void Clear()
{
while( size ) keys[PopHeap()] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
}
bool Contains(INMOST_DATA_ENUM_TYPE i)
{
return index[i] != ENUMUNDEF;
}
};
//this structure is used to provide sorting routing for Reverse-Cuthill-McKee reordering
struct RCM_Comparator
{
INMOST_DATA_ENUM_TYPE wbeg;
std::vector<INMOST_DATA_ENUM_TYPE> & xadj;
public:
RCM_Comparator(INMOST_DATA_ENUM_TYPE _wbeg, std::vector<INMOST_DATA_ENUM_TYPE> & _xadj)
: wbeg(_wbeg), xadj(_xadj) {}
RCM_Comparator(const RCM_Comparator & b) : wbeg(b.wbeg), xadj(b.xadj) {}
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];}
};
}
#endif //INMOST_UTILS_H
Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.hpp
View file @ 84d9eeda
......@@ -14,7 +14,6 @@ class MLMTILUC_preconditioner : public Method
INMOST_DATA_ENUM_TYPE Size() { return last - first; }
} Interval;
typedef dynarray<INMOST_DATA_ENUM_TYPE,256> levels_t;
//Sparse::Vector DL0,DR0;
std::vector<Sparse::Row::entry> LU_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;
......
Source/Solver/solver_inner/solver_mptilu2/solver_mtilu2.cpp
View file @ 84d9eeda
#include "inmost_solver.h"
#if defined(USE_SOLVER)
#include "solver_mtilu2.hpp"
#include "../../../Misc/utils.h"
#define DEFAULT_TAU 0.005
#define DEFAULT_TAU2 0.00001
......@@ -13,124 +14,6 @@
#define REORDER_MAXIMUM_TRANSVERSAL
template<class T> struct make_integer;
template<> struct make_integer<float> {typedef int type;};
template<> struct make_integer<double> {typedef long long type;};
__INLINE static bool compare(INMOST_DATA_REAL_TYPE * a, INMOST_DATA_REAL_TYPE * b)
{
return (*reinterpret_cast< make_integer<INMOST_DATA_REAL_TYPE>::type * >(a)) <=
(*reinterpret_cast< make_integer<INMOST_DATA_REAL_TYPE>::type * >(b));
}
class BinaryHeap
{
INMOST_DATA_REAL_TYPE * Base;
std::vector<INMOST_DATA_REAL_TYPE *> Array;
std::vector<INMOST_DATA_ENUM_TYPE> Position;
public:
void Clear()
{
while(!Array.empty())
{
Position[Array.back()-Base] = ENUMUNDEF;
Array.pop_back();
}
}
INMOST_DATA_REAL_TYPE * Get(INMOST_DATA_ENUM_TYPE pos) {return Array[pos];}
INMOST_DATA_ENUM_TYPE GetSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(Array.size());}
INMOST_DATA_ENUM_TYPE GetPosition(INMOST_DATA_ENUM_TYPE pos)
{
return Position[pos];
}
INMOST_DATA_ENUM_TYPE DecreaseKey(INMOST_DATA_ENUM_TYPE pos)
{
INMOST_DATA_ENUM_TYPE i = Position[pos];
++i;
while(i > 1)
{
//if((*Array[i-1]) <= (*Array[i/2-1]))
if( compare(Array[i-1],Array[i/2-1]) )
{
Position[(Array[i/2-1]-Base)] = i-1;
Position[(Array[i-1]-Base)] = i/2-1;
std::swap(Array[i/2-1],Array[i-1]);
}
else break;
i = i/2;
}
return i;
}
INMOST_DATA_ENUM_TYPE PushHeap(INMOST_DATA_REAL_TYPE * key)
{
INMOST_DATA_ENUM_TYPE i = GetSize();
Array.push_back(key);
Position[(key-Base)] = i;
++i;
while(i > 1)
{
//if((*Array[i-1]) <= (*Array[i/2-1]) )
if( compare(Array[i-1],Array[i/2-1]) )
{
Position[(Array[i-1]-Base)] = i/2-1;
Position[(Array[i/2-1]-Base)] = i-1;
std::swap(Array[i-1],Array[i/2-1]);
}
else break;
i = i/2;
}
return i;
}
void BalanceHeap(INMOST_DATA_ENUM_TYPE i)
{
INMOST_DATA_ENUM_TYPE Index;
++i;
while(i <= Array.size()/2)
{
if( 2*i+1 > Array.size() )
Index = 2*i;
//else if( (*Array[2*i-1]) <= (*Array[2*i+1-1]) )
else if( compare(Array[2*i-1],Array[2*i+1-1]) )
Index = 2*i;
else
Index = 2*i+1;
//if(!((*Array[i-1]) <= (*Array[Index-1])))
if(!compare(Array[i-1],Array[Index-1]))
{
Position[(Array[i-1]-Base)] = Index-1;
Position[(Array[Index-1]-Base)] = i-1;
std::swap(Array[i-1],Array[Index-1]);
}
else break;
i = Index;
}
}
INMOST_DATA_ENUM_TYPE PopHeap()
{
INMOST_DATA_ENUM_TYPE Ret = ENUMUNDEF;
if(Array.empty()) return Ret;
Ret = static_cast<INMOST_DATA_ENUM_TYPE>(Array[0]-Base);
Array[0] = Array.back();
Position[Array[0] - Base] = 0;
Array.pop_back();
Position[Ret] = ENUMUNDEF;
BalanceHeap(0);
return Ret;
}
BinaryHeap(INMOST_DATA_REAL_TYPE * Base, INMOST_DATA_ENUM_TYPE Size)
: Base(Base)
{
Position.resize(Size,ENUMUNDEF);
Array.reserve(4096);
}
~BinaryHeap()
{
}
};
void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & Address,
std::vector<Sparse::Row::entry> & Entries,
......@@ -453,13 +336,12 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
}
else if( l < Dist[Ui] )
{
Dist[Ui] = l;
Parent[Perm[Ui]] = Li;
AugmentPosition[Ui] = Lit;
if( Heap.GetPosition(Ui-mobeg) != ENUMUNDEF )
Heap.DecreaseKey(Ui-mobeg);
if( Heap.Contains(Ui-mobeg) )
Heap.DecreaseKey(Ui-mobeg,l);
else
Heap.PushHeap(&Dist[Ui]);
Heap.PushHeap(Ui-mobeg,l);
}
}
}
......@@ -510,7 +392,6 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
Trace = Parent[Trace];
}
for(Ui = 0; Ui < Heap.GetSize(); ++Ui) *Heap.Get(Ui) = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
Heap.Clear();
}
}
......
Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.cpp
View file @ 84d9eeda
......@@ -6,6 +6,7 @@
#include <sstream>
#include <deque>
#include <iomanip>
#include "../../../Misc/utils.h"
//#define REPORT_ILU
//#undef REPORT_ILU
#define REPORT_ILU_PROGRESS
......@@ -69,137 +70,7 @@ using namespace INMOST;
#include <Mondriaan.h>
#endif
template<class T> struct make_integer;
template<> struct make_integer<float> {typedef int type;};
template<> struct make_integer<double> {typedef long long type;};
__INLINE static bool compare(INMOST_DATA_REAL_TYPE * a, INMOST_DATA_REAL_TYPE * b)
{
return (*reinterpret_cast< make_integer<INMOST_DATA_REAL_TYPE>::type * >(a)) <=
(*reinterpret_cast< make_integer<INMOST_DATA_REAL_TYPE>::type * >(b));
}
//this structure is used to provide sorting routing for Reverse-Cuthill-McKee reordering
struct RCM_Comparator
{
INMOST_DATA_ENUM_TYPE wbeg;
std::vector<INMOST_DATA_ENUM_TYPE> & xadj;
public:
RCM_Comparator(INMOST_DATA_ENUM_TYPE _wbeg, std::vector<INMOST_DATA_ENUM_TYPE> & _xadj)
: wbeg(_wbeg), xadj(_xadj) {}
RCM_Comparator(const RCM_Comparator & b) : wbeg(b.wbeg), xadj(b.xadj) {}
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];}
};
class BinaryHeap1
{
INMOST_DATA_REAL_TYPE * Base;
std::vector<INMOST_DATA_REAL_TYPE *> Array;
std::vector<INMOST_DATA_ENUM_TYPE> Position;
public:
void Clear()
{
while(!Array.empty())
{
Position[Array.back()-Base] = ENUMUNDEF;
Array.pop_back();
}
}
INMOST_DATA_REAL_TYPE * Get(INMOST_DATA_ENUM_TYPE pos) {return Array[pos];}
INMOST_DATA_ENUM_TYPE GetSize() {return static_cast<INMOST_DATA_ENUM_TYPE>(Array.size());}
INMOST_DATA_ENUM_TYPE GetPosition(INMOST_DATA_ENUM_TYPE pos)
{
return Position[pos];
}
INMOST_DATA_ENUM_TYPE DecreaseKey(INMOST_DATA_ENUM_TYPE pos)
{
INMOST_DATA_ENUM_TYPE i = Position[pos];
++i;
while(i > 1)
{
//if((*Array[i-1]) <= (*Array[i/2-1]))