Skip to content

GitLab

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

MPT permutation for structurally singular matrix

parent 8f29e471
Hide whitespace changes
Inline Side-by-side
Source/Solver/solver.cpp
View file @ cc1eb6b3
......@@ -46,8 +46,8 @@
//#define USE_OMP
#define KSOLVER BCGSL_solver
//#define KSOLVER BCGS_solver
//#define KSOLVER BCGSL_solver
#define KSOLVER BCGS_solver
//#define ACCELERATED_CONDEST
//#define PRINT_CONDEST
......
Source/Solver/solver_mtilu2.cpp
View file @ cc1eb6b3
......@@ -515,44 +515,65 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
//for(k = mobeg; k != moend; ++k)
//if( Perm[k] != k ) std::cout << "Column " << k << " go to " << Perm[k] << " inverse " << IPerm[k] << std::endl;
{ //fill gaps in perm
std::fill(IPerm.begin(), IPerm.end(), ENUMUNDEF);
for(k = mobeg; k < moend; ++k)
{
if( Perm[k] != ENUMUNDEF )
IPerm[Perm[k]] = 0;
}
std::vector<INMOST_DATA_ENUM_TYPE> gaps;
for(k = mobeg; k < moend; ++k)
if( IPerm[k] == ENUMUNDEF )
gaps.push_back(k);
for(k = mobeg; k < moend; ++k)
if( Perm[k] == ENUMUNDEF )
{
Perm[k] = gaps.back();
gaps.pop_back();
}
}
cnt = 0;
#if defined(RESCALE_MAXIMUM_TRANSVERSAL)
trescale = Timer();
for (k = mobeg; k < moend; ++k)
{
bool flip_sign = false;
for(INMOST_DATA_ENUM_TYPE qt = 0; qt < A[k].Size(); ++qt)
{
INMOST_DATA_ENUM_TYPE i = A[k].GetIndex(qt), j = Perm[i];
INMOST_DATA_REAL_TYPE l = exp(V[k]), u = exp(U[i])/Cmax[i];
DL[k] = l;
DR[i] = u;
B_Entries[cnt++] = Sparse::Row::make_entry(j, l*u*A[k].GetValue(qt));
if( j == k && (B_Entries[cnt-1].second) < 0.0 ) flip_sign = true;
}
if( flip_sign )
{
DL[k] *= -1;
for(INMOST_DATA_ENUM_TYPE Li = B_Address[k]; Li < B_Address[k+1]; ++Li) B_Entries[Li].second *= -1;
}
//std::sort(B_Entries.begin()+C_Address[k],B_Entries.begin()+C_Address[k+1]);
}
for (k = mobeg; k < moend; ++k) U[Perm[k]] = DR[k];
for (k = mobeg; k < moend; ++k) DR[k] = U[k];
trescale = Timer() - trescale;
trescale = Timer();
for (k = mobeg; k < moend; ++k)
{
bool flip_sign = false;
for(INMOST_DATA_ENUM_TYPE qt = 0; qt < A[k].Size(); ++qt)
{
INMOST_DATA_ENUM_TYPE i = A[k].GetIndex(qt), j = Perm[i];
INMOST_DATA_REAL_TYPE l = exp(V[k]), u = exp(U[i])/Cmax[i];
DL[k] = l;
DR[i] = u;
B_Entries[cnt++] = Sparse::Row::make_entry(j, l*u*A[k].GetValue(qt));
if( j == k && (B_Entries[cnt-1].second) < 0.0 ) flip_sign = true;
}
if( flip_sign )
{
DL[k] *= -1;
for(INMOST_DATA_ENUM_TYPE Li = B_Address[k]; Li < B_Address[k+1]; ++Li) B_Entries[Li].second *= -1;
}
//std::sort(B_Entries.begin()+C_Address[k],B_Entries.begin()+C_Address[k+1]);
}
for (k = mobeg; k < moend; ++k) U[Perm[k]] = DR[k];
for (k = mobeg; k < moend; ++k) DR[k] = U[k];
trescale = Timer() - trescale;
#else
for (k = mobeg; k < moend; ++k)
{
for(INMOST_DATA_ENUM_TYPE qt = 0; qt < A[k].Size(); ++qt)
{
if( A[k].GetIndex(qt) >= mobeg && A[k].GetIndex(qt) < moend )
B_Entries[cnt++] = Sparse::Row::make_entry(Perm[A[k].GetIndex(qt)], A[k].GetValue(qt));
}
//std::sort(B_Entries.begin()+C_Address[k],B_Entries.begin()+C_Address[k+1]);
}
for (k = mobeg; k < moend; ++k)
{
for(INMOST_DATA_ENUM_TYPE qt = 0; qt < A[k].Size(); ++qt)
{
if( A[k].GetIndex(qt) >= mobeg && A[k].GetIndex(qt) < moend )
B_Entries[cnt++] = Sparse::Row::make_entry(Perm[A[k].GetIndex(qt)], A[k].GetValue(qt));
}
//std::sort(B_Entries.begin()+C_Address[k],B_Entries.begin()+C_Address[k+1]);
}
#endif
//for (k = mobeg; k < moend; ++k) IPerm[Perm[k]] = k; //inversePQ
......@@ -563,21 +584,21 @@ void MTILU2_preconditioner::DumpMatrix(interval<INMOST_DATA_ENUM_TYPE, INMOST_DA
}
#else
{
treorder = Timer();
INMOST_DATA_ENUM_TYPE cnt = 0;
B_Address[mobeg] = 0;
for(k = mobeg; k < moend; ++k)
{
Perm[k] = k;
for (Sparse::Row::iterator it = A[k].Begin(); it != A[k].End(); ++it)
{
if( it->first >= mobeg && it->first < moend )
B_Entries[cnt++] = Sparse::Row::make_entry(it->first,it->second);
}
B_Address[k+1] = cnt;
}
//DumpMatrix(B_Address,B_Entries,mobeg,moend,"NOMC64.mtx");
treorder = Timer() - treorder;
treorder = Timer();
INMOST_DATA_ENUM_TYPE cnt = 0;
B_Address[mobeg] = 0;
for(k = mobeg; k < moend; ++k)
{
Perm[k] = k;
for (Sparse::Row::iterator it = A[k].Begin(); it != A[k].End(); ++it)
{
if( it->first >= mobeg && it->first < moend )
B_Entries[cnt++] = Sparse::Row::make_entry(it->first,it->second);
}
B_Address[k+1] = cnt;
}
//DumpMatrix(B_Address,B_Entries,mobeg,moend,"NOMC64.mtx");
treorder = Timer() - treorder;
}
#endif
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
......
Source/Solver/solver_mtiluc2.cpp
View file @ cc1eb6b3
......@@ -905,6 +905,26 @@ public:
//DumpMatrix(B_Address,B_Entries,wbeg,wend,"MC64.mtx");
std::fill(localP.begin() + (wbeg - mobeg), localP.begin() + (wend - mobeg), ENUMUNDEF);
{ //check that there are no gaps in Perm
for(k = cbeg; k < cend; ++k)
{
if( Perm[k] != ENUMUNDEF )
localP[Perm[k]] = 0;
}
std::vector<INMOST_DATA_ENUM_TYPE> gaps;
for(k = cbeg; k < cend; ++k)
if( localP[k] == ENUMUNDEF )
gaps.push_back(k);
for(k = cbeg; k < cend; ++k)
if( Perm[k] == ENUMUNDEF )
{
Perm[k] = gaps.back();
gaps.pop_back();
}
std::fill(localP.begin() + (wbeg - mobeg), localP.begin() + (wend - mobeg), ENUMUNDEF);
}
//exit(-1);
......@@ -1973,7 +1993,7 @@ swap_algorithm:
LineValuesU[k] = B_Entries[B_Address[k].first].second;
else
{
std::cout << __LINE__ << " No diagonal value! " << k << " " << B_Entries[B_Address[k].first].first << std::endl;
//std::cout << __LINE__ << " No diagonal value! " << k << " " << B_Entries[B_Address[k].first].first << std::endl;
LineValuesU[k] = 0.0;
}
......@@ -2179,7 +2199,7 @@ swap_algorithm:
LineValuesL[k] = B_Entries[B_Address[k].first].second;
else
{
std::cout << __LINE__ << " No diagonal value! " << k << std::endl;
//std::cout << __LINE__ << " No diagonal value! " << k << std::endl;
LineValuesL[k] = 0.0;
}
//start from diagonal
......
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