Commit d98ecb00 authored by Kirill Terekhov's avatar Kirill Terekhov

Colors in DrawMatrix and MPTILUC2

parent 2f21c004
...@@ -17,8 +17,9 @@ Sparse::Matrix * m = NULL; ...@@ -17,8 +17,9 @@ Sparse::Matrix * m = NULL;
int zoom = 200; int zoom = 200;
int block_size = 0; int block_size = 0;
int draw_color = 0;
//~ Storage::real min = 1e20,max = -1e20; double amin = 1e20,amax = -1e20;
void printtext(const char * fmt, ...) void printtext(const char * fmt, ...)
{ {
...@@ -379,6 +380,18 @@ void keyboard(unsigned char key, int x, int y) ...@@ -379,6 +380,18 @@ void keyboard(unsigned char key, int x, int y)
} }
} }
if( key == 'd' )
{
draw_color = (draw_color+1)%5;
switch(draw_color)
{
case 0: std::cout << "color diagonal" << std::endl; break;
case 1: std::cout << "color by min:max" << std::endl; break;
case 2: std::cout << "color by min:max in log scale" << std::endl; break;
case 3: std::cout << "color by min:max per row" << std::endl; break;
case 4: std::cout << "color by min:max per row in log scale" << std::endl; break;
}
}
if (key == 'c') ord->clear(); if (key == 'c') ord->clear();
if (key == 'r') if (key == 'r')
{ {
...@@ -415,53 +428,128 @@ void draw() ...@@ -415,53 +428,128 @@ void draw()
glVertex2i(m->Size(), m->Size() - ord->GetMatrixPartSize()); glVertex2i(m->Size(), m->Size() - ord->GetMatrixPartSize());
glEnd(); glEnd();
*/ */
glColor3f(0,0,0);
glBegin(GL_QUADS);
for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
if( jt->first != it - m->Begin() && jt->second)
//DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(jt->first));//, sqrt((jt->second-min)/(max-min)));
DrawEntry(jt->first, m->Size() - (it - m->Begin()));//, sqrt((jt->second-min)/(max-min)));
glEnd();
/* /*
glColor3f(0.0, 1.0, 0); glColor3f(0.0, 1.0, 0);
glBegin(GL_QUADS); glBegin(GL_QUADS);
for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it) for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{
int ind = it - m->Begin();
if (fabs((*it)[ind]) > row_sum[it-m->Begin()]) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind));
DrawEntry((it - m->Begin()), m->Size() - ind);
}
glEnd();
*/
if( draw_color == 0 )
{ {
int ind = it - m->Begin();
if (fabs((*it)[ind]) > row_sum[it-m->Begin()]) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind)); glColor3f(0,0,0);
DrawEntry((it - m->Begin()), m->Size() - ind); glBegin(GL_QUADS);
for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
if( jt->first != it - m->Begin() && jt->second)
{
//double r = jt->second
//DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(jt->first));//, sqrt((jt->second-min)/(max-min)));
DrawEntry(jt->first, m->Size() - (it - m->Begin()));//, sqrt((jt->second-min)/(max-min)));
}
glEnd();
glBegin(GL_QUADS);
for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{
int ind = it - m->Begin();
double t = fabs((*it)[ind]) / row_sum[ind];
//if (fabs((*it)[ind]) < row_sum[ind]) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind));
glColor3f(0.0, t, 1.0-t);
DrawEntry(ind, m->Size() - ind);
}
glEnd();
//zoom += 1;
glColor3f(1.0, 0, 0);
glBegin(GL_QUADS);
for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{
int ind = it - m->Begin();
if (fabs((*it)[ind]) < 1e-13) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind));
DrawEntry((it - m->Begin()), m->Size() - ind);
}
glEnd();
//zoom -= 1;
} }
glEnd(); else if( draw_color == 1 )
*/
glBegin(GL_QUADS);
for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{ {
glBegin(GL_QUADS);
int ind = it - m->Begin(); for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
double t = fabs((*it)[ind]) / row_sum[ind]; for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
//if (fabs((*it)[ind]) < row_sum[ind]) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind)); {
glColor3f(0.0, t, 1.0-t); double r = (fabs(jt->second) - amin)/(amax-amin);
DrawEntry(ind, m->Size() - ind); glColor3f(1-r,1-r,1-r);
DrawEntry(jt->first, m->Size() - (it - m->Begin()));
}
glEnd();
} }
glEnd(); else if( draw_color == 2 )
//zoom += 1;
glColor3f(1.0, 0, 0);
glBegin(GL_QUADS);
for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{ {
int ind = it - m->Begin(); glBegin(GL_QUADS);
if (fabs((*it)[ind]) < 1e-13) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind)); for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
DrawEntry((it - m->Begin()), m->Size() - ind); for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
{
double r = log((fabs(jt->second) - amin)/(amax-amin)*63 + 1)/log(64);
glColor3f(1-r,1-r,1-r);
DrawEntry(jt->first, m->Size() - (it - m->Begin()));
}
glEnd();
}
else if( draw_color == 3 )
{
glBegin(GL_QUADS);
for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{
double rmax = -1.0e20, rmin = 1.0e20;
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
{
if( fabs(jt->second) < rmin ) rmin = fabs(jt->second);
if( fabs(jt->second) > rmax ) rmax = fabs(jt->second);
}
if( rmin > rmax ) continue;
//if( rmax - rmin < 1.0e-13 ) continue;
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
{
double r = (fabs(jt->second) - rmin)/(rmax-rmin);
glColor3f(1-r,1-r,1-r);
DrawEntry(jt->first, m->Size() - (it - m->Begin()));
}
}
glEnd();
}
else if( draw_color == 4 )
{
glBegin(GL_QUADS);
for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
{
double rmax = -1.0e20, rmin = 1.0e20;
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
{
if( fabs(jt->second) < rmin ) rmin = fabs(jt->second);
if( fabs(jt->second) > rmax ) rmax = fabs(jt->second);
}
if( rmin > rmax ) continue;
//if( rmax - rmin < 1.0e-13 ) continue;
for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
{
double r = log((fabs(jt->second) - rmin)/(rmax-rmin)*63 + 1)/log(64);
glColor3f(1-r,1-r,1-r);
DrawEntry(jt->first, m->Size() - (it - m->Begin()));
}
}
glEnd();
} }
glEnd();
//zoom -= 1;
if (CommonInput != NULL) if (CommonInput != NULL)
{ {
...@@ -575,23 +663,23 @@ int main(int argc, char ** argv) ...@@ -575,23 +663,23 @@ int main(int argc, char ** argv)
std::cout << "Nonzeros: " << nnz << std::endl; std::cout << "Nonzeros: " << nnz << std::endl;
zoom = m->Size() / 1000; zoom = 1;//m->Size() / 1000;
//~ for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it) for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
//~ for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt) for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
//~ { {
//~ if(jt->second < min ) if(fabs(jt->second) < amin )
//~ { {
//~ min = jt->second; amin = fabs(jt->second);
//~ std::cout << (it-m->Begin()) << " " << jt->first << " " << jt->second << std::endl; //std::cout << (it-m->Begin()) << " " << jt->first << " " << jt->second << std::endl;
//~ } }
//~ if(jt->second > max ) if(fabs(jt->second) > amax )
//~ { {
//~ max = jt->second; amax = fabs(jt->second);
//~ std::cout << (it-m->Begin()) << " " << jt->first << " " << jt->second << std::endl; //std::cout << (it-m->Begin()) << " " << jt->first << " " << jt->second << std::endl;
//~ } }
//~ } }
//~ std::cout << "max: " << max << " min: " << min << std::endl; std::cout << "absolute max: " << amax << " min: " << amin << std::endl;
glutInit(&argc,argv); glutInit(&argc,argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
glutInitWindowSize(width, height); glutInitWindowSize(width, height);
......
...@@ -613,6 +613,8 @@ static bool allow_pivot = true; ...@@ -613,6 +613,8 @@ static bool allow_pivot = true;
{ {
//ddPQ reordering on current Schur complement //ddPQ reordering on current Schur complement
INMOST_DATA_ENUM_TYPE cbeg = wbeg, cend = wend; //next size of factored B block INMOST_DATA_ENUM_TYPE cbeg = wbeg, cend = wend; //next size of factored B block
//DumpMatrix(A_Address,A_Entries,cbeg,cend,"mat"+to_string(level_size.size())+".mtx");
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
/// MAXIMUM TRANSVERSE REORDERING /// /// MAXIMUM TRANSVERSE REORDERING ///
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
...@@ -1672,6 +1674,7 @@ static bool allow_pivot = true; ...@@ -1672,6 +1674,7 @@ static bool allow_pivot = true;
DR0.Save("DR0_"+to_string(level_size.size())+".mtx"); DR0.Save("DR0_"+to_string(level_size.size())+".mtx");
DumpMatrix(B_Address,B_Entries,cbeg,cend,"mat_equil"+to_string(level_size.size())+".mtx"); DumpMatrix(B_Address,B_Entries,cbeg,cend,"mat_equil"+to_string(level_size.size())+".mtx");
*/ */
//DumpMatrix(B_Address,B_Entries,cbeg,cend,"rescale"+to_string(level_size.size())+".mtx");
//exit(-1); //exit(-1);
//rescale EF //rescale EF
if( !level_size.empty() ) if( !level_size.empty() )
...@@ -1825,17 +1828,21 @@ static bool allow_pivot = true; ...@@ -1825,17 +1828,21 @@ static bool allow_pivot = true;
//assert(B_Entries[B_Address[k].first].first == k); //assert(B_Entries[B_Address[k].first].first == k);
LineIndecesU[cbeg] = k; LineIndecesU[cbeg] = k;
LineIndecesU[k] = EOL; LineIndecesU[k] = EOL;
if (B_Entries[B_Address[k].first].first == k) LineValuesU[k] = 0.0;
LineValuesU[k] = B_Entries[B_Address[k].first].second; if( B_Address[k].Size() )
else
LineValuesU[k] = 0.0;
Ui = k;
for (INMOST_DATA_ENUM_TYPE it = B_Address[k].first + (B_Entries[B_Address[k].first].first == k ? 1 : 0); it < B_Address[k].last; ++it)
{ {
LineValuesU[B_Entries[it].first] = B_Entries[it].second; if (B_Entries[B_Address[k].first].first == k)
Ui = LineIndecesU[Ui] = B_Entries[it].first; LineValuesU[k] = B_Entries[B_Address[k].first].second;
else
LineValuesU[k] = 0.0;
Ui = k;
for (INMOST_DATA_ENUM_TYPE it = B_Address[k].first + (B_Entries[B_Address[k].first].first == k ? 1 : 0); it < B_Address[k].last; ++it)
{
LineValuesU[B_Entries[it].first] = B_Entries[it].second;
Ui = LineIndecesU[Ui] = B_Entries[it].first;
}
LineIndecesU[Ui] = EOL;
} }
LineIndecesU[Ui] = EOL;
Sbeg = k; Sbeg = k;
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// U-part elimination with L // // U-part elimination with L //
...@@ -2027,11 +2034,14 @@ static bool allow_pivot = true; ...@@ -2027,11 +2034,14 @@ static bool allow_pivot = true;
indicesU.push_back(Ui); indicesU.push_back(Ui);
Ui = LineIndecesU[Ui]; Ui = LineIndecesU[Ui];
} }
std::sort(indicesU.begin(),indicesU.end()); if( !indicesU.empty() )
//Sbeg = indices[0]; {
for(size_t qt = 1; qt < indicesU.size(); ++qt) std::sort(indicesU.begin(),indicesU.end());
LineIndecesU[indicesU[qt-1]] = indicesU[qt]; //Sbeg = indices[0];
LineIndecesU[indicesU.back()] = EOL; for(size_t qt = 1; qt < indicesU.size(); ++qt)
LineIndecesU[indicesU[qt-1]] = indicesU[qt];
LineIndecesU[indicesU.back()] = EOL;
}
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// Retrieve diagonal value // // Retrieve diagonal value //
...@@ -2093,21 +2103,25 @@ static bool allow_pivot = true; ...@@ -2093,21 +2103,25 @@ static bool allow_pivot = true;
//insert diagonal value first //insert diagonal value first
LineIndecesL[cbeg] = k; LineIndecesL[cbeg] = k;
LineIndecesL[k] = EOL; LineIndecesL[k] = EOL;
if (B_Entries[B_Address[k].first].first == k) LineValuesL[k] = 0.0;
LineValuesL[k] = B_Entries[B_Address[k].first].second; if( B_Address[k].Size() )
else
LineValuesL[k] = 0.0;
//start from diagonal
Ui = k;
Li = Bbeg[k];
while (Li != EOL)
{ {
assert(B_Entries[B_Address[Li].first].first == k); if (B_Entries[B_Address[k].first].first == k)
LineValuesL[Li] = B_Entries[B_Address[Li].first].second; LineValuesL[k] = B_Entries[B_Address[k].first].second;
Ui = LineIndecesL[Ui] = Li; else
Li = Blist[Li]; LineValuesL[k] = 0.0;
//start from diagonal
Ui = k;
Li = Bbeg[k];
while (Li != EOL)
{
assert(B_Entries[B_Address[Li].first].first == k);
LineValuesL[Li] = B_Entries[B_Address[Li].first].second;
Ui = LineIndecesL[Ui] = Li;
Li = Blist[Li];
}
LineIndecesL[Ui] = EOL;
} }
LineIndecesL[Ui] = EOL;
Sbeg = k; Sbeg = k;
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// L-part elimination with U // // L-part elimination with U //
...@@ -2301,11 +2315,14 @@ static bool allow_pivot = true; ...@@ -2301,11 +2315,14 @@ static bool allow_pivot = true;
indicesL.push_back(Ui); indicesL.push_back(Ui);
Ui = LineIndecesL[Ui]; Ui = LineIndecesL[Ui];
} }
std::sort(indicesL.begin(),indicesL.end()); if( !indicesL.empty() )
//Sbeg = indices[0]; {
for(size_t qt = 1; qt < indicesL.size(); ++qt) std::sort(indicesL.begin(),indicesL.end());
LineIndecesL[indicesL[qt-1]] = indicesL[qt]; //Sbeg = indices[0];
LineIndecesL[indicesL.back()] = EOL; for(size_t qt = 1; qt < indicesL.size(); ++qt)
LineIndecesL[indicesL[qt-1]] = indicesL[qt];
LineIndecesL[indicesL.back()] = EOL;
}
//check that diagonal value is the same(must be!) //check that diagonal value is the same(must be!)
//assert(fabs(LineValues[k] / udiag - 1.0) < 1.0e-10); //assert(fabs(LineValues[k] / udiag - 1.0) < 1.0e-10);
...@@ -3687,7 +3704,7 @@ static bool allow_pivot = true; ...@@ -3687,7 +3704,7 @@ static bool allow_pivot = true;
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// sort contents of linked list // // sort contents of linked list //
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
/*
indicesS.clear(); indicesS.clear();
Ui = Sbeg; Ui = Sbeg;
while(Ui != EOL) while(Ui != EOL)
...@@ -3695,34 +3712,45 @@ static bool allow_pivot = true; ...@@ -3695,34 +3712,45 @@ static bool allow_pivot = true;
indicesS.push_back(Ui); indicesS.push_back(Ui);
Ui = LineIndecesL[Ui]; Ui = LineIndecesL[Ui];
} }
std::sort(indicesS.begin(),indicesS.end()); if( !indicesS.empty() )
*/ {
//Sbeg = indices[0]; std::sort(indicesS.begin(),indicesS.end());
//for(size_t qt = 1; qt < indices.size(); ++qt)
// LineIndecesL[indices[qt-1]] = indices[qt]; Sbeg = indicesS[0];
//LineIndecesL[indices.back()] = EOL; for(size_t qt = 1; qt < indicesS.size(); ++qt)
LineIndecesL[indicesS[qt-1]] = indicesS[qt];
LineIndecesL[indicesS.back()] = EOL;
}
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// prepare row norm for dropping // // prepare row norm for dropping //
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
#if defined(SCHUR_DROPPING_S) #if defined(SCHUR_DROPPING_S)
INMOST_DATA_REAL_TYPE Snorm = 0, Snum = 0, tauS = std::min(tau*tau,tau2)/std::max(NuU_max,NuL_max); INMOST_DATA_REAL_TYPE Snorm = 0, Snum = 0, tauS, Smax = 0, Smin = 1.0e+54;
Ui = Sbeg; Ui = Sbeg;
while (Ui != EOL) while (Ui != EOL)
//for(size_t qt = 0; qt < indicesS.size(); ++qt) //for(size_t qt = 0; qt < indicesS.size(); ++qt)
{ {
//Ui = indicesS[qt]; //Ui = indicesS[qt];
u = LineValuesL[Ui]; u = fabs(LineValuesL[Ui]);
if( u > Smax ) Smax = u;
if( u < Smin ) Smin = u;
Snorm += u*u; Snorm += u*u;
Snum++; Snum++;
Ui = LineIndecesL[Ui]; Ui = LineIndecesL[Ui];
} }
if( Snum ) Snorm = sqrt(Snorm/Snum); if( Snum ) Snorm = sqrt(Snorm/Snum);
Snorm = std::min(1.0,Snorm); Snorm = std::min(1.0,Snorm);
//tauS = std::min(tau*tau,tau2)/std::max(NuU_max,NuL_max)*Snorm;
tauS = Smin - 1.0e-7 + (Smax-Smin)*std::min(tau*tau,tau2);
//exp(log((fabs(u)-Smin)/(Smax-Smin)*127+1)) > exp(0.5*log(128))
//fabs(u) > exp(0.5)/127*(Smax-Smin) + Smax
#endif #endif
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// put calculated row to Schur complement // // put calculated row to Schur complement //
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
INMOST_DATA_REAL_TYPE Sdrop = 0;
Ui = Sbeg; Ui = Sbeg;
S_Address[k].first = static_cast<INMOST_DATA_ENUM_TYPE>(S_Entries.size()); S_Address[k].first = static_cast<INMOST_DATA_ENUM_TYPE>(S_Entries.size());
while (Ui != EOL) while (Ui != EOL)
...@@ -3731,19 +3759,34 @@ static bool allow_pivot = true; ...@@ -3731,19 +3759,34 @@ static bool allow_pivot = true;
//Ui = indicesS[qt]; //Ui = indicesS[qt];
u = LineValuesL[Ui]; u = LineValuesL[Ui];
#if defined(SCHUR_DROPPING_S) #if defined(SCHUR_DROPPING_S)
//if( fabs(u)*NuU*NuL > tau * Snorm ) //if( fabs(Smax-Smin) < 1.0e-5 || log((fabs(u)-Smin)/(Smax-Smin)*127+1)/log(128) > 0.5 )
//if( fabs(u)*NuU*NuL*NuU_acc*NuL_acc*NuD*NuD_acc > tau * Snorm ) if( fabs(u) > tauS )
if( fabs(u) > tauS * Snorm )
//if( fabs(u) > tau*tau * Snorm )
#else #else
if( 1+u != 1 ) if( 1+u != 1 )
#endif #endif
S_Entries.push_back(Sparse::Row::make_entry(Ui,u)); S_Entries.push_back(Sparse::Row::make_entry(Ui,u));
else else
{
Sdrop += u;
ndrops_s++; ndrops_s++;
}
Ui = LineIndecesL[Ui]; Ui = LineIndecesL[Ui];
} }
S_Address[k].last = static_cast<INMOST_DATA_ENUM_TYPE>(S_Entries.size()); S_Address[k].last = static_cast<INMOST_DATA_ENUM_TYPE>(S_Entries.size());
///////////////////////////////////////////////////////////////////////////////////
// Substract dropped elements from maximal element to keep sum of row //
// (TODO: may collect dropping per row/column and preserve row-column sum) //
///////////////////////////////////////////////////////////////////////////////////
if( Sdrop )
{
INMOST_DATA_ENUM_TYPE rmax = S_Address[k].first, r;
INMOST_DATA_REAL_TYPE Sadd = Sdrop / (double)S_Address[k].Size();
for(r = S_Address[k].first; r < S_Address[k].last; ++r)
S_Entries[r].second += Sadd;
// if( fabs(S_Entries[r].second) > fabs(S_Entries[rmax].second) )
// rmax = r;
//S_Entries[rmax].second += Sdrop;
}
//assert(std::is_sorted(S_Entries.begin()+S_Addres[k].first,S_Entries.end())); //assert(std::is_sorted(S_Entries.begin()+S_Addres[k].first,S_Entries.end()));
/////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////
// clean up temporary linked list // // clean up temporary linked list //
......
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