Commit 3c199079 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Reverse-Cuthill-McKee reordering in MPTILUC preconditioner

Added Reverse-Cuthill-McKee reordering into INNER_MPTILUC preconditioner
as default reordering when Metis_NodeND is not provided. On the tests
RCM shows 10% more fill-in in inverse factors then with Metis_NodeND,
but the RCM reordering procedure is 4 times faster then Metis_NodeND and
overall computational time decreases with RCM.

Fix type mismatch warnings in autodiff.cpp
parent 9c538219
......@@ -255,8 +255,8 @@ namespace INMOST
Storage::integer_array indarr = jt->IntegerArray(it->indices);
for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
{
if( *qt < first_num ) Pre.insert(*qt);
else if( *qt >= last_num ) Post.insert(*qt);
if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
}
}
}
......@@ -271,8 +271,8 @@ namespace INMOST
indarr.resize(jt->RealArray(it->d.t).size());
for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
{
if( *qt < first_num ) Pre.insert(*qt);
else if( *qt >= last_num ) Post.insert(*qt);
if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
}
}
}
......@@ -289,8 +289,8 @@ namespace INMOST
Storage::integer_array indarr = jt->IntegerArray(it->indices);
for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
{
if( *qt < first_num ) Pre.insert(*qt);
else if( *qt >= last_num ) Post.insert(*qt);
if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
}
}
}
......@@ -304,8 +304,8 @@ namespace INMOST
Storage::integer_array indarr = jt->IntegerArray(it->indices);
for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
{
if( *qt < first_num ) Pre.insert(*qt);
else if( *qt >= last_num ) Post.insert(*qt);
if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
}
}
}
......
......@@ -4,6 +4,7 @@
#if defined(USE_SOLVER)
#include "solver_mtiluc2.hpp"
#include <sstream>
#include <deque>
//#define REPORT_ILU
//#undef REPORT_ILU
//#define REPORT_ILU_PROGRESS
......@@ -20,7 +21,7 @@ using namespace INMOST;
#endif
//#define REORDER_RCM
#define REORDER_RCM
//#define REORDER_NNZ
#if defined(USE_SOLVER_METIS)
#define REORDER_METIS_ND
......@@ -77,6 +78,20 @@ __INLINE static bool compare(INMOST_DATA_REAL_TYPE * a, INMOST_DATA_REAL_TYPE *
(*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 BinaryHeap
{
......@@ -507,6 +522,9 @@ public:
double tfactor = 0.0, trescale = 0.0, treorder = 0.0, treassamble = 0.0, ttotal, tt;
#if defined(REORDER_METIS_ND)
double tmetisgraph = 0, tmetisnd = 0;
#endif
#if defined(REORDER_RCM)
double trcmgraph = 0, trcmorder = 0;
#endif
double tlfactor, tlrescale, tlreorder, tlreassamble;
ttotal = Timer();
......@@ -936,6 +954,7 @@ public:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////// END MAXIMUM TRANSVERSE REORDERING /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
tt = Timer();
#if defined(REORDER_DDPQ)
#if defined(REPORT_ILU)
......@@ -1152,7 +1171,6 @@ public:
cend = wend;
i = wend;
treorder += Timer() - tt;
#elif defined(REORDER_NNZ)
//for (k = cbeg; k < cend; ++k) V[k] = DR[localQ[k]];
//for (k = cbeg; k < cend; ++k) DR[k] = V[k];
......@@ -1202,31 +1220,153 @@ public:
}
cend = i;
#elif defined(REORDER_RCM)
//compute order of each entry
std::fill(Ulist.begin() + wbeg - mobeg, Ulist.begin() + wend - mobeg, 0);
std::fill(Llist.begin() + wbeg - mobeg, Llist.begin() + wend - mobeg, 0);
for (k = wbeg; k < wend; ++k)
{
for (INMOST_DATA_ENUM_TYPE it = A_Address[k].first; it != A_Address[k].last; ++it)
//create a symmetric graph of the matrix A + A^T
std::vector<INMOST_DATA_ENUM_TYPE> xadj(wend-wbeg+1), adjncy;
adjncy.reserve(nzA*2);
for (k = cbeg; k < cend; ++k) if( localQ[k] == ENUMUNDEF ) printf("%s:%d No column permutation for row %d. Matrix is structurally singular\n",__FILE__,__LINE__,k);
for (k = cbeg; k < cend; ++k) localP[k] = k;
for (k = cbeg; k < cend; ++k) V[localQ[k]] = DR[k]; //if you expirience a bug here then the matrix is structurally singular
for (k = cbeg; k < cend; ++k) DR[k] = V[k];
for(k = wbeg; k < wend; ++k)
{
i = A_Entries[it].first;
//nonzeros
Ulist[k]++; //row nnz
Llist[i]++; //col nnz
for (INMOST_DATA_ENUM_TYPE jt = A_Address[k].first; jt != A_Address[k].last; ++jt)
{
A_Entries[jt].first = localQ[A_Entries[jt].first];
//A_Entries[jt].second *= DL[k]*DR[A_Entries[jt].first];
}
std::sort(A_Entries.begin()+A_Address[k].first,A_Entries.begin()+A_Address[k].last);
}
inversePQ(wbeg,wend,localP,localQ, invP,invQ);
applyPQ(wbeg, wend, localP, localQ, invP, invQ);
trcmgraph = Timer();
for (k = wend; k > wbeg; --k)
{
//vwgt[k-1] = 1;
if (A_Address[k-1].Size() > 0)
{
i = A_Entries[A_Address[k-1].first].first;
Blist[k-1] = Bbeg[i];
Bbeg[i] = k-1;
}
Ulist[k-1] = A_Address[k-1].first;
}
xadj[0] = 0;
for(i = wbeg; i < wend; ++i)
{
for (INMOST_DATA_ENUM_TYPE jt = A_Address[i].first; jt != A_Address[i].last; ++jt)
{
if( A_Entries[jt].first != i )
adjncy.push_back(static_cast<INMOST_DATA_ENUM_TYPE>(A_Entries[jt].first));
}
Li = Bbeg[i];
while (Li != EOL)
{
if( Li != i )
adjncy.push_back(static_cast<INMOST_DATA_ENUM_TYPE>(Li));
Li = Blist[Li];
}
Li = Bbeg[i];
while (Li != EOL)
{
Ui = Blist[Li];
Ulist[Li]++;
if (A_Address[Li].last - Ulist[Li] > 0)
{
k = A_Entries[Ulist[Li]].first;
Blist[Li] = Bbeg[k];
Bbeg[k] = Li;
}
Li = Ui;
}
std::sort(adjncy.begin()+xadj[i-wbeg],adjncy.end());
adjncy.resize(std::unique(adjncy.begin()+xadj[i-wbeg],adjncy.end())-adjncy.begin());
xadj[i-wbeg+1] = static_cast<INMOST_DATA_ENUM_TYPE>(adjncy.size());
}
std::fill(Bbeg.begin(),Bbeg.end(),EOL);
trcmgraph = Timer()-trcmgraph;
trcmorder = Timer();
std::fill(Ulist.begin() + wbeg - mobeg, Ulist.begin() + wend - mobeg, ENUMUNDEF);
//find node with the lowest order
INMOST_DATA_ENUM_TYPE start = wbeg;
INMOST_DATA_ENUM_TYPE index = wbeg;
INMOST_DATA_ENUM_TYPE cur = ENUMUNDEF;
std::deque<INMOST_DATA_ENUM_TYPE> q;
std::vector<INMOST_DATA_ENUM_TYPE> conns;
//for(k = wbeg+1; k < wend; ++k)
// if( RCM_Comparator(wbeg,xadj)(k,start) )
// start = k;
//Ulist[start] = index++;
do
{
cur = ENUMUNDEF;
for(k = wbeg; k < wend && cur == ENUMUNDEF; ++k)
{
if( Ulist[k] == ENUMUNDEF )
cur = k;
}
assert(cur != ENUMUNDEF);
for(k = cur+1; k < wend; ++k) if( Ulist[k] == ENUMUNDEF )
{
if( RCM_Comparator(wbeg,xadj)(k,cur) )
cur = k;
}
q.push_back(cur);
Ulist[cur] = index++;
while(!q.empty())
{
cur = q.front();
q.pop_front();
for (INMOST_DATA_ENUM_TYPE it = xadj[cur-wbeg]; it < xadj[cur-wbeg+1]; ++it)
if( Ulist[adjncy[it]] == ENUMUNDEF )
conns.push_back(adjncy[it]);
std::sort(conns.begin(),conns.end(),RCM_Comparator(wbeg,xadj));
for (k = 0; k < static_cast<INMOST_DATA_ENUM_TYPE>(conns.size()); ++k)
{
Ulist[conns[k]] = index++;
q.push_back(conns[k]);
}
conns.clear();
}
}
while( index < wend );
for(k = wbeg; k < wend; ++k)
Ulist[k] = wend-(Ulist[k]-wbeg)-1;
for(k = wbeg; k < wend; ++k)
if( Ulist[k] + Llist[Blist[k]] < Ulist[start] + Llist[Blist[start]] )
start = k;
localP[start] = localQ[Blist[start]] = index++;
{
localP[k] = Ulist[k];//wend - Ulist[k] + 1;
localQ[k] = Ulist[k];//wend - Ulist[k] + 1;
//localP[Ulist[k]] = k;
//localQ[Ulist[k]] = k;
}
cend = wend;
i = wend;
trcmorder = Timer() - trcmorder;
}
#else
cend = wend;
i = wbeg;
#endif
tt = Timer() - tt;
treorder += tt;
if (cbeg == cend && cbeg != wend)
{
//std::cout << __FILE__ << ":" << __LINE__ << " singular matrix, factored " << mobeg << ".." << cend << " out of " << mobeg << ".." << moend << std::endl;
......@@ -1299,7 +1439,21 @@ public:
nzA += A_Address[k].Size();
}
//DumpMatrix(B_Address,B_Entries,cbeg,cend,"mat_reorder.mtx");
/*
{
int rank = 0;
#if defined(USE_MPI)
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
#endif
std::stringstream str;
str << "mat_reorder" << rank << ".mtx";
DumpMatrix(B_Address,B_Entries,cbeg,cend,str.str());
}
#if defined(USE_MPI)
MPI_Barrier(MPI_COMM_WORLD);
#endif
exit(0);
*/
//Setup column addressing for B,F, in descending order to keep list ordered
for (k = cend; k > cbeg; --k)
{
......@@ -1622,9 +1776,9 @@ public:
#if defined(ILUC2)
U2_Address[cbeg].last = static_cast<INMOST_DATA_ENUM_TYPE>(LU2_Entries.size());
#endif
///////////////////////////////////////////////////////////////////////////////////
// form first line of L and second-order L //
///////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////
// form first line of L and second-order L //
///////////////////////////////////////////////////////////////////////////////////
L_Address[cbeg].first = static_cast<INMOST_DATA_ENUM_TYPE>(LU_Entries.size());
#if defined(ILUC2)
L2_Address[cbeg].first = static_cast<INMOST_DATA_ENUM_TYPE>(LU2_Entries.size());
......@@ -3032,10 +3186,14 @@ swap_algorithm:
printf("reorder %f (%6.2f%%)\n", treorder, 100.0*treorder/ttotal);
#if defined(REORDER_METIS_ND)
printf("metis graph %f nd %f\n", tmetisgraph, tmetisnd);
#endif
#if defined(REORDER_RCM)
printf("rcm graph %f reorder %f\n", trcmgraph, trcmorder);
#endif
printf("reassamble %f (%6.2f%%)\n", treassamble, 100.0*treassamble / ttotal);
printf("rescale %f (%6.2f%%)\n", trescale, 100.0*trescale / ttotal);
printf("factor %f (%6.2f%%)\n", tfactor, 100.0*tfactor / ttotal);
printf("nnz A %d LU %d LU2 %d\n",nzA,nzLU,nzLU2);
#endif
init = true;
/*
......
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