Skip to content

GitLab

Commit ee6715b1 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

some refactoring in MLMPTILUC2

parent 3231377a
Hide whitespace changes
Inline Side-by-side
Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.cpp
View file @ ee6715b1
......@@ -67,10 +67,11 @@ const INMOST_DATA_ENUM_TYPE UNDEF = ENUMUNDEF, EOL = ENUMUNDEF - 1;
#define NNZ_GROWTH_PARAM 1.05
#if defined(REORDER_METIS_ND)
#if defined(USE_SOLVER_METIS)
#define METIS_EXPORT
#include "metis.h"
#endif
#if defined(REORDER_ZOLTAN_HUND)
#include <zoltan.h>
#endif
......@@ -1353,448 +1354,841 @@ const INMOST_DATA_ENUM_TYPE UNDEF = ENUMUNDEF, EOL = ENUMUNDEF - 1;
}
bool MLMTILUC_preconditioner::isInitialized() { return init; }
bool MLMTILUC_preconditioner::isFinalized() { return !init; }
bool MLMTILUC_preconditioner::Initialize()
void MLMTILUC_preconditioner::MaximalTransversal( INMOST_DATA_ENUM_TYPE wbeg,
INMOST_DATA_ENUM_TYPE wend,
interval<INMOST_DATA_ENUM_TYPE, Interval> & A_Address,
std::vector<Sparse::Row::entry> & A_Entries,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localP,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & localQ,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & U,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & V,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & DL,
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> & DR)
{
if (verbosity > 1) std::cout << __FILE__ << ":" << __LINE__ << " mem " << getCurrentRSS() << " peak " << getPeakRSS() << std::endl;
int nthreads = Threads();
const INMOST_DATA_REAL_TYPE subst = 1.0; (void)subst;
const INMOST_DATA_REAL_TYPE tol_modif = PIVOT_THRESHOLD_VALUE;
bool block_pivot = false;
INMOST_DATA_ENUM_TYPE wbeg, wend; //working interval
INMOST_DATA_ENUM_TYPE mobeg, moend; // total interval
INMOST_DATA_ENUM_TYPE vbeg, vend; // vector interval
INMOST_DATA_INTEGER_TYPE k;
INMOST_DATA_ENUM_TYPE i, j, Li, Ui, curr, next;
INMOST_DATA_REAL_TYPE l,u, max_diag = 0, min_diag = 0;
INMOST_DATA_ENUM_TYPE nzA, nzLU = 0, nzA0;
//Sparse::Vector DL, DR;
Sparse::Vector DL0,DR0;
info->GetOverlapRegion(info->GetRank(), mobeg, moend);
info->GetVectorRegion(vbeg,vend);
//prepare reordering vectors
ddP.set_interval_beg(mobeg);
ddP.set_interval_end(moend);
ddQ.set_interval_beg(mobeg);
ddQ.set_interval_end(moend);
//prepare rescaling vectors
//DL.SetInterval(mobeg, moend);
//DR.SetInterval(mobeg, moend);
DL0.SetInterval(mobeg, moend);
DR0.SetInterval(mobeg, moend);
for(k = mobeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(moend); ++k)
DL0[k] = DR0[k] = 1.0;
//Sparse::Vector & U = DL;
//Sparse::Vector & V = DR;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> Cmax(wbeg,wend,0.0);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & Perm = localQ;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & IPerm = localP;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> ColumnList(wbeg,wend);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> Parent(wbeg,wend);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> AugmentPosition(wbeg,wend,ENUMUNDEF);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> ColumnPosition(wbeg,wend,ENUMUNDEF);
std::vector<INMOST_DATA_REAL_TYPE> C_Entries(A_Entries.size());
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> Dist(wbeg, wend + 1, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
for (k = mobeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(moend); k++)
ddP[k] = ddQ[k] = k;
//interval<INMOST_DATA_ENUM_TYPE, bool> Pivot(mobeg,moend,false);
BinaryHeap Heap(&Dist[wbeg],wend-wbeg);
// arrays U,V,Dist are cleared at the end of schur complement calculation
//std::fill(U.begin() + wbeg - mobeg, U.begin() + wend - mobeg, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
//std::fill(V.begin() + wbeg - mobeg, V.begin() + wend - mobeg, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
//std::fill(Cmax.begin() + wbeg - mobeg, Cmax.begin() + wend - mobeg, 0.0);
//std::fill(Dist.begin() + wbeg - mobeg, Dist.begin() + wend + 1 - mobeg, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
Perm[k] = ENUMUNDEF;
IPerm[k] = ENUMUNDEF;
Parent[k] = ENUMUNDEF;
ColumnList[k] = ENUMUNDEF;
}
if (verbosity > 1) std::cout << __FILE__ << ":" << __LINE__ << " mem " << getCurrentRSS() << " peak " << getPeakRSS() << std::endl;
//supplimentary data structures for ILUC
INMOST_DATA_ENUM_TYPE L_Beg, U_Beg, Sbeg;
U_Address.set_interval_beg(mobeg);
U_Address.set_interval_end(moend);
L_Address.set_interval_beg(mobeg);
L_Address.set_interval_end(moend);
LU_Diag.set_interval_beg(mobeg);
LU_Diag.set_interval_end(moend);
std::vector<Sparse::Row::entry> A_Entries;
//std::vector<INMOST_DATA_ENUM_TYPE> sort_indeces;
interval<INMOST_DATA_ENUM_TYPE, Interval> A_Address(mobeg, moend);
#if defined(NNZ_PIVOT)
INMOST_DATA_ENUM_TYPE nnzL = 0, nnzU = 0, nnzL_max = 0, nnzU_max = 0;
#endif
INMOST_DATA_ENUM_TYPE nnzE = 0, nnzF = 0, nnzA = 0;
INMOST_DATA_REAL_TYPE NuU = 1, NuL = 1, NuD = 1, NuU_max = 1.0, NuL_max = 1.0;
INMOST_DATA_REAL_TYPE NuU_acc = 1, NuL_acc = 1, NuD_acc = 1;
//supplimentary data structures for returning values of dropped elements
//~ INMOST_DATA_REAL_TYPE DropLk, DropUk, SumLk, SumUk;
//~ interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> DropU(mobeg,moend,0.0), DropL(mobeg,moend,0.0);
//~ #if defined(ILUC2)
//~ interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> DropU2(mobeg,moend,0.0), DropL2(mobeg,moend,0.0);
//~ #endif
//data structure for linked list
//interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> U(mobeg,moend,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
//interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> V(mobeg,moend,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
//interval<INMOST_DATA_ENUM_TYPE, std::vector< std::pair<INMOST_DATA_ENUM_TYPE,INMOST_DATA_ENUM_TYPE> > > TransposeU(mobeg,moend), TransposeL(mobeg,moend);
double tfactor = 0.0, trescale = 0.0, treorder = 0.0, ttransversal = 0.0, treassamble = 0.0, ttotal, tt, testimator = 0.0, tschur = 0.0, tlocal;
#if defined(REORDER_METIS_ND)
double tmetisgraph = 0, tmetisnd = 0;
#endif
#if defined(REORDER_RCM) || defined(REORDER_WRCM) || defined(REORDER_BRCM)
double trcmgraph = 0, trcmorder = 0;
#endif
double tlfactor, tlrescale, tlreorder, tlreassamble, tlschur;
ttotal = Timer();
//(*Alink).Save("M.mtx");
//calculate number of nonzeros
nzA = 0;
for (k = mobeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(moend); ++k)
// Initial LOG transformation to dual problem and initial extreme match
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
for (Sparse::Row::iterator r = (*Alink)[k].Begin(); r != (*Alink)[k].End(); ++r)
if (r->first >= mobeg && r->first < moend && fabs(r->second) > 0.0) nzA++;
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
INMOST_DATA_ENUM_TYPE i = A_Entries[it].first;
INMOST_DATA_REAL_TYPE u = C_Entries[it] = fabs(A_Entries[it].second);
if( u > Cmax[i] ) Cmax[i] = u;
//C_Entries.push_back(Sparse::Row::make_entry(i,u));
}
}
nzA0 = nzA;
//sort_indeces.reserve(256);
A_Entries.resize(nzA);
//LU_Entries.reserve(nzA*4);
j = 0;
for (k = mobeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(moend); ++k)
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
A_Address[k].first = j;
for (Sparse::Row::iterator r = (*Alink)[k].Begin(); r != (*Alink)[k].End(); ++r)
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
if (r->first >= mobeg && r->first < moend && fabs(r->second) > 0 )//1.0e-30 )// != 0.0)
A_Entries[j++] = Sparse::Row::make_entry(r->first, r->second);
INMOST_DATA_ENUM_TYPE i = A_Entries[it].first;
if( Cmax[i] == 0 || C_Entries[it] == 0 )
C_Entries[it] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
else
{
C_Entries[it] = log(Cmax[i])-log(C_Entries[it]);
if( C_Entries[it] < U[i] ) U[i] = C_Entries[it];
}
}
A_Address[k].last = j;
//assert(A_Address[k].Size() != 0); //singular matrix
}
if (verbosity > 1) std::cout << __FILE__ << ":" << __LINE__ << " mem " << getCurrentRSS() << " peak " << getPeakRSS() << std::endl;
INMOST_DATA_REAL_TYPE Unorm, Lnorm;
//~ INMOST_DATA_REAL_TYPE Unum, Lnum;
INMOST_DATA_ENUM_TYPE nzLU2 = 0, nzLU2tot = 0, ndrops = 0;
#if defined(ILUC2)
INMOST_DATA_REAL_TYPE tau2 = iluc2_tau;
#else
INMOST_DATA_REAL_TYPE tau2 = tau;
#endif
for(k = mobeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(moend); ++k)
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
U_Address[k].first = U_Address[k].last = ENUMUNDEF;
L_Address[k].first = L_Address[k].last = ENUMUNDEF;
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
INMOST_DATA_REAL_TYPE u = C_Entries[it] - U[A_Entries[it].first];
if( u < V[k] ) V[k] = u;
}
}
if( verbosity > 1 )
std::cout << "nonzeros in A " << nzA << " pivot cond " << pivot_cond << " diag " << pivot_diag << " tau " << tau << " tau2 " << tau2 << std::endl;
int swaps = 0, totswaps = 0;
//set up working interval
wbeg = mobeg;
wend = moend;
while (wbeg < wend) //iterate into levels until all matrix is factored
/// Update cost and match
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
//ddPQ reordering on current Schur complement
INMOST_DATA_ENUM_TYPE cbeg = wbeg, cend = wend; //next size of factored B block
if (verbosity > 1) std::cout << __FILE__ << ":" << __LINE__ << " mem " << getCurrentRSS() << " peak " << getPeakRSS() << std::endl;
//this scope is for reordering and rescaling
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
interval<INMOST_DATA_ENUM_TYPE, bool> donePQ(mobeg, moend, false);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> invP(wbeg, wend), invQ(wbeg, wend);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> localP(wbeg, wend, ENUMUNDEF), localQ(wbeg, wend, ENUMUNDEF);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> U(wbeg,wend,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> V(wbeg,wend,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> DL(wbeg,wend,1.0), DR(wbeg,wend,1.0);
//this scope is for reorderings
INMOST_DATA_REAL_TYPE u = fabs(C_Entries[it] - V[k] - U[A_Entries[it].first]);
if( u < 1.0e-30 && Perm[A_Entries[it].first] == ENUMUNDEF && IPerm[k] == ENUMUNDEF )
{
//DumpMatrix(A_Address,A_Entries,cbeg,cend,"mat"+to_string(level_size.size())+".mtx");
/// MAXIMUM TRANSVERSE REORDERING
if( run_mpt )
Perm[A_Entries[it].first] = k;
IPerm[k] = A_Entries[it].first;
ColumnPosition[k] = it;
}
}
}
/// 1-step augmentation
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
if( IPerm[k] == ENUMUNDEF ) //unmatched row
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last && IPerm[k] == ENUMUNDEF; ++it)
{
INMOST_DATA_REAL_TYPE u = fabs(C_Entries[it] - V[k] - U[A_Entries[it].first]);
if( u <= 1.0e-30 )
{
if( verbosity > 1 )
std::cout << "Reordering with MPT\n";
ttransversal = Timer();
INMOST_DATA_ENUM_TYPE ColumnBegin;
INMOST_DATA_ENUM_TYPE pop_heap_pos;
INMOST_DATA_REAL_TYPE ShortestPath, AugmentPath;
INMOST_DATA_ENUM_TYPE PathEnd, Trace, IPermPrev;
//Sparse::Vector & U = DL;
//Sparse::Vector & V = DR;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> Cmax(wbeg,wend,0.0);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & Perm = localQ;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> & IPerm = localP;
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> ColumnList(wbeg,wend);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> Parent(wbeg,wend);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> AugmentPosition(wbeg,wend,ENUMUNDEF);
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_ENUM_TYPE> ColumnPosition(wbeg,wend,ENUMUNDEF);
std::vector<INMOST_DATA_REAL_TYPE> C_Entries(A_Entries.size());
interval<INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE> Dist(wbeg, wend + 1, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
BinaryHeap Heap(&Dist[wbeg],wend-wbeg);
// arrays U,V,Dist are cleared at the end of schur complement calculation
//std::fill(U.begin() + wbeg - mobeg, U.begin() + wend - mobeg, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
//std::fill(V.begin() + wbeg - mobeg, V.begin() + wend - mobeg, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
//std::fill(Cmax.begin() + wbeg - mobeg, Cmax.begin() + wend - mobeg, 0.0);
//std::fill(Dist.begin() + wbeg - mobeg, Dist.begin() + wend + 1 - mobeg, std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
Perm[k] = ENUMUNDEF;
IPerm[k] = ENUMUNDEF;
Parent[k] = ENUMUNDEF;
ColumnList[k] = ENUMUNDEF;
}
if (verbosity > 1) std::cout << __FILE__ << ":" << __LINE__ << " mem " << getCurrentRSS() << " peak " << getPeakRSS() << std::endl;
// Initial LOG transformation to dual problem and initial extreme match
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
INMOST_DATA_ENUM_TYPE Li = Perm[A_Entries[it].first];
assert(Li != ENUMUNDEF);
// Search other row in C for 0
for (INMOST_DATA_ENUM_TYPE Lit = A_Address[Li].first; Lit < A_Address[Li].last; ++Lit)
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
u = fabs(C_Entries[Lit]- V[Li] - U[A_Entries[Lit].first]);
if( u <= 1.0e-30 && Perm[A_Entries[Lit].first] == ENUMUNDEF )
{
i = A_Entries[it].first;
u = C_Entries[it] = fabs(A_Entries[it].second);
if( u > Cmax[i] ) Cmax[i] = u;
//C_Entries.push_back(Sparse::Row::make_entry(i,u));
}
}
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
i = A_Entries[it].first;
if( Cmax[i] == 0 || C_Entries[it] == 0 )
C_Entries[it] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
else
{
C_Entries[it] = log(Cmax[i])-log(C_Entries[it]);
if( C_Entries[it] < U[i] ) U[i] = C_Entries[it];
}
}
}
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
u = C_Entries[it] - U[A_Entries[it].first];
if( u < V[k] ) V[k] = u;
Perm[A_Entries[it].first] = k;
IPerm[k] = A_Entries[it].first;
ColumnPosition[k] = it;
Perm[A_Entries[Lit].first] = Li;
IPerm[Li] = A_Entries[Lit].first;
ColumnPosition[Li] = Lit;
break;
}
}
/// Update cost and match
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
}
}
}
}
/// Weighted bipartite matching
for(INMOST_DATA_INTEGER_TYPE k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
if( IPerm[k] != ENUMUNDEF )
continue;
INMOST_DATA_ENUM_TYPE Li = k;
INMOST_DATA_ENUM_TYPE ColumnBegin = EOL;
Parent[Li] = ENUMUNDEF;
INMOST_DATA_ENUM_TYPE PathEnd = ENUMUNDEF;
INMOST_DATA_ENUM_TYPE Trace = k;
INMOST_DATA_REAL_TYPE ShortestPath = 0;
INMOST_DATA_REAL_TYPE AugmentPath = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
while(true)
{
for (INMOST_DATA_ENUM_TYPE Lit = A_Address[Li].first; Lit < A_Address[Li].last; ++Lit)
{
INMOST_DATA_ENUM_TYPE Ui = A_Entries[Lit].first;
//if( ColumnList[Ui] == k ) continue;
if( ColumnList[Ui] != ENUMUNDEF ) continue;
INMOST_DATA_REAL_TYPE l = fabs(ShortestPath + C_Entries[Lit] - V[Li] - U[Ui]);
//if( l < 0.0 ) printf("row %d col %d negative l %g Augment %lf Shortest %lf C %lf V %lf U %lf\n",k,Ui,l,AugmentPath,ShortestPath,C_Entries[Lit],V[Li],U[Ui]);
if( l < 0.0 && l > -1.0e-8 ) l = 0;
if( l < 0.0 ) continue;
if( l < AugmentPath )
{
if( Perm[Ui] == ENUMUNDEF )
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last; ++it)
{
u = fabs(C_Entries[it] - V[k] - U[A_Entries[it].first]);
if( u < 1.0e-30 && Perm[A_Entries[it].first] == ENUMUNDEF && IPerm[k] == ENUMUNDEF )
{
Perm[A_Entries[it].first] = k;
IPerm[k] = A_Entries[it].first;
ColumnPosition[k] = it;
}
}
PathEnd = Ui;
Trace = Li;
AugmentPath = l;
AugmentPosition[Ui] = Lit;
}
/// 1-step augmentation
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
else if( l < Dist[Ui] )
{
if( IPerm[k] == ENUMUNDEF ) //unmatched row
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it < A_Address[k].last && IPerm[k] == ENUMUNDEF; ++it)
{
u = fabs(C_Entries[it] - V[k] - U[A_Entries[it].first]);
if( u <= 1.0e-30 )
{
Li = Perm[A_Entries[it].first];
assert(Li != ENUMUNDEF);
// Search other row in C for 0
for (INMOST_DATA_ENUM_TYPE Lit = A_Address[Li].first; Lit < A_Address[Li].last; ++Lit)
{
u = fabs(C_Entries[Lit]- V[Li] - U[A_Entries[Lit].first]);
if( u <= 1.0e-30 && Perm[A_Entries[Lit].first] == ENUMUNDEF )
{
Perm[A_Entries[it].first] = k;
IPerm[k] = A_Entries[it].first;
ColumnPosition[k] = it;
Perm[A_Entries[Lit].first] = Li;
IPerm[Li] = A_Entries[Lit].first;
ColumnPosition[Li] = Lit;
break;
}
}
}
}
}
Parent[Perm[Ui]] = Li;
AugmentPosition[Ui] = Lit;
if( Heap.Contains(Ui-wbeg) )
Heap.DecreaseKey(Ui-wbeg,l);
else
Heap.PushHeap(Ui-wbeg,l);
}
/// Weighted bipartite matching
for(k = wbeg; k < static_cast<INMOST_DATA_INTEGER_TYPE>(wend); ++k)
{
if( IPerm[k] != ENUMUNDEF )
continue;
Li = k;
ColumnBegin = EOL;
Parent[Li] = ENUMUNDEF;
PathEnd = ENUMUNDEF;
Trace = k;
ShortestPath = 0;
AugmentPath = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
while(true)
{
for (INMOST_DATA_ENUM_TYPE Lit = A_Address[Li].first; Lit < A_Address[Li].last; ++Lit)
{
Ui = A_Entries[Lit].first;
//if( ColumnList[Ui] == k ) continue;
if( ColumnList[Ui] != ENUMUNDEF ) continue;
l = fabs(ShortestPath + C_Entries[Lit] - V[Li] - U[Ui]);
//if( l < 0.0 ) printf("row %d col %d negative l %g Augment %lf Shortest %lf C %lf V %lf U %lf\n",k,Ui,l,AugmentPath,ShortestPath,C_Entries[Lit],V[Li],U[Ui]);
if( l < 0.0 && l > -1.0e-8 ) l = 0;
if( l < 0.0 ) continue;
if( l < AugmentPath )
{
if( Perm[Ui] == ENUMUNDEF )
{
PathEnd = Ui;
Trace = Li;
AugmentPath = l;
AugmentPosition[Ui] = Lit;
}
else if( l < Dist[Ui] )
{
Parent[Perm[Ui]] = Li;
AugmentPosition[Ui] = Lit;
if( Heap.Contains(Ui-wbeg) )
Heap.DecreaseKey(Ui-wbeg,l);
else
Heap.PushHeap(Ui-wbeg,l);
}
}
}
}
}
pop_heap_pos = Heap.PopHeap();
if( pop_heap_pos == ENUMUNDEF ) break;
Ui = pop_heap_pos+wbeg;
ShortestPath = Dist[Ui];
INMOST_DATA_ENUM_TYPE pop_heap_pos = Heap.PopHeap();
if( pop_heap_pos == ENUMUNDEF ) break;
INMOST_DATA_ENUM_TYPE Ui = pop_heap_pos+wbeg;
ShortestPath = Dist[Ui];
if( AugmentPath <= ShortestPath )
{
Dist[Ui] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
//Heap.increaseKey(Ui,Dist[Ui]);
break;
}
if( AugmentPath <= ShortestPath )
{
Dist[Ui] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
//Heap.increaseKey(Ui,Dist[Ui]);
break;
}
ColumnList[Ui] = ColumnBegin;
ColumnBegin = Ui;
ColumnList[Ui] = ColumnBegin;
ColumnBegin = Ui;
Li = Perm[Ui];
}
if( PathEnd != ENUMUNDEF )
{
Ui = ColumnBegin;
while(Ui != EOL)
{
U[Ui] += Dist[Ui] - AugmentPath;
if( Perm[Ui] != ENUMUNDEF ) V[Perm[Ui]] = C_Entries[ColumnPosition[Perm[Ui]]] - U[Ui];
Dist[Ui] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
Li = ColumnList[Ui];
ColumnList[Ui] = ENUMUNDEF;
Ui = Li;
}
Li = Perm[Ui];
}
if( PathEnd != ENUMUNDEF )
{
INMOST_DATA_ENUM_TYPE Ui = ColumnBegin;
while(Ui != EOL)
{
U[Ui] += Dist[Ui] - AugmentPath;
if( Perm[Ui] != ENUMUNDEF ) V[Perm[Ui]] = C_Entries[ColumnPosition[Perm[Ui]]] - U[Ui];
Dist[Ui] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
Li = ColumnList[Ui];
ColumnList[Ui] = ENUMUNDEF;
Ui = Li;
}
Ui = PathEnd;
while(Trace != ENUMUNDEF)
{
INMOST_DATA_ENUM_TYPE IPermPrev = IPerm[Trace];
Perm[Ui] = Trace;
IPerm[Trace] = Ui;
ColumnPosition[Trace] = AugmentPosition[Ui];
V[Trace] = C_Entries[ColumnPosition[Trace]] - U[Ui];
Ui = IPermPrev;
Trace = Parent[Trace];