Permalink
Browse files

Improvements

Pseudoinverse for polynomial matrix in bicgstab(l) to improve stability
(as in petsc);
Added RowMerger class for sum of sparse rows.
  • Loading branch information...
kirill-terekhov committed Mar 18, 2015
1 parent 6177cfc commit 82230d19a385dcbb7c7515ff40b23cc664aa5dcb
Showing with 840 additions and 273 deletions.
  1. +1 −1 autodiff.cpp
  2. +3 −3 examples/MatSolve/main.cpp
  3. +5 −5 inmost_autodiff.h
  4. +87 −1 inmost_solver.h
  5. +120 −0 solver.cpp
  6. +403 −78 solver_bcgsl.hpp
  7. +7 −3 solver_ddpqiluc2.cpp
  8. +4 −2 solver_ilu2.hpp
  9. +204 −178 tests.txt
  10. +2 −2 tests/solver_test001/main.cpp
  11. +4 −0 tests/solver_test002/CMakeLists.txt
View
@@ -19,7 +19,7 @@ namespace INMOST
void Automatizator::DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Solver::Row & entries, INMOST_DATA_REAL_TYPE multval, void * user_data)
{
INMOST_DATA_ENUM_TYPE voffset = var.values_offset(element), doffset = var.derivatives_offset(element);
INMOST_DATA_ENUM_TYPE k = var.data.size()-1;
INMOST_DATA_ENUM_TYPE k = static_cast<INMOST_DATA_ENUM_TYPE>(var.data.size()-1);
Storage e = Storage(m,var.current_stencil[element].first);
INMOST_DATA_REAL_TYPE lval, rval, ret;
var.values[doffset+k] = multval;
@@ -94,7 +94,7 @@ int main(int argc, char ** argv)
t = Timer();
success = s.Solve(b,x); // Solve the linear system with the previously computted preconditioner
BARRIER
if( !rank ) std::cout << "solver: " << Timer() - t << std::endl;
if( !rank ) std::cout << "solver: " << Timer() - t << "\t\t\t" << std::endl;
iters = s.Iterations(); // Get the number of iterations performed
resid = s.Residual(); // Get the final residual achieved
reason = s.GetReason();
@@ -125,9 +125,9 @@ int main(int argc, char ** argv)
double temp[2] = {aresid,bresid}, recv[2] = {aresid,bresid};
#if defined(USE_MPI)
MPI_Reduce(temp,recv,2,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
if( info.GetRank() == 0 ) std::cout << "||Ax-b|| " << sqrt(recv[0]) << " ||b|| " << sqrt(recv[1]) << " ||Ax-b||/||b|| " << sqrt(recv[0]/recv[1]) << std::endl;
if( info.GetRank() == 0 ) std::cout << "||Ax-b|| " << sqrt(recv[0]) << " ||b|| " << sqrt(recv[1]) << " ||Ax-b||/||b|| " << sqrt(recv[0]/(recv[1]+1.0e-100)) << std::endl;
#endif
realresid = sqrt(recv[0]/recv[1]);
realresid = sqrt(recv[0]/(recv[1]+1.0e-100));
//realresid = sqrt(realresid);
info.RestoreVector(x);
View
@@ -156,7 +156,7 @@ namespace INMOST
__INLINE INMOST_DATA_REAL_TYPE GetStaticValue(const Storage & e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->RealArray(GetStaticValueTag(ind))[comp]; }
__INLINE bool isStaticValid(const Storage & e, INMOST_DATA_ENUM_TYPE ind) { MarkerType mask = GetStaticMask(ind); return mask == 0 || e->GetMarker(mask); }
#if defined(NEW_VERSION)
INMOST_DATA_REAL_TYPE Evaluate(expr & var, const const Storage & e, void * user_data);
INMOST_DATA_REAL_TYPE Evaluate(expr & var, const Storage & e, void * user_data);
INMOST_DATA_REAL_TYPE Derivative(expr & var, const Storage & e, Solver::Row & out, Storage::real multiply, void * user_data);
#else
INMOST_DATA_REAL_TYPE Evaluate(const expr & var, const Storage & e, void * user_data);
@@ -547,17 +547,17 @@ namespace INMOST
expr(const expr * l, const expr * r) :data() { data.push_back(expr_data(l,r)); }
expr(INMOST_DATA_REAL_TYPE val) : data() { data.push_back(expr_data(val)); }
expr(INMOST_DATA_ENUM_TYPE op, INMOST_DATA_ENUM_TYPE comp) : data() { data.push_back(expr_data(op, comp, ENUMUNDEF)); }
expr(INMOST_DATA_ENUM_TYPE op, const expr & operand) : data(operand.data) { relink_data(); data.push_back(expr_data(op, data.size() - 1, ENUMUNDEF)); }
expr(const expr & operand, const expr & multiplyer) : data(operand.data) { relink_data(); data.push_back(expr_data(AD_VAL, data.size() - 1, multiplyer)); }
expr(INMOST_DATA_ENUM_TYPE op, const expr & operand) : data(operand.data) { relink_data(); data.push_back(expr_data(op, static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), ENUMUNDEF)); }
expr(const expr & operand, const expr & multiplyer) : data(operand.data) { relink_data(); data.push_back(expr_data(AD_VAL, static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), multiplyer)); }
expr(const expr & cond, const expr & if_true, const expr & if_false) :data(cond.data)
{
relink_data();
data.push_back(expr_data(data.size() - 1, expr_data(&if_true, &if_false)));
data.push_back(expr_data(static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1), expr_data(&if_true, &if_false)));
}
expr(INMOST_DATA_ENUM_TYPE op, const expr & l, const expr & r) :data(l.data)
{
relink_data();
INMOST_DATA_ENUM_TYPE lp = data.size() - 1;
INMOST_DATA_ENUM_TYPE lp = static_cast<INMOST_DATA_ENUM_TYPE>(data.size() - 1);
INMOST_DATA_ENUM_TYPE rp = merge_data(r.data);
data.push_back(expr_data(op, lp, rp));
}
View
@@ -225,12 +225,15 @@ namespace INMOST
ret.second = val;
return ret;
}
private:
typedef dynarray<entry,16> Entries; //replace later with more memory-efficient chunk_array, with first chunk in stack
//typedef array<entry> Entries;
//typedef std::vector<entry> Entries;
//typedef sparse_data<INMOST_DATA_ENUM_TYPE,INMOST_DATA_REAL_TYPE> Entries;
//typedef Entries::pair entry; //for sparse_data
public:
typedef Entries::iterator iterator;
typedef Entries::const_iterator const_iterator;
typedef Entries::reverse_iterator reverse_iterator;
@@ -348,6 +351,9 @@ namespace INMOST
/// @param size New size of the row.
void Resize(INMOST_DATA_ENUM_TYPE size) {data.resize(size);}
};
/// Class to store the distributed sparse matrix by compressed rows.
/// The format used to store sparse matrix is analogous to Compressed Row Storage format (CRS).
@@ -428,6 +434,86 @@ namespace INMOST
std::string GetName() {return name;}
//~ friend class Solver;
};
/// This class may be used to sum multiple sparse rows.
/// \warning
/// In parallel column indices of the matrix may span wider then
/// local row indices, to prevent any problem you are currently
/// advised to set total size of the matrix as interval of the
/// RowMerger. In future this may change, see todo 2 below.
/// \todo
/// 1. Add iterators over entries.
/// 2. Implement multiple intervals for distributed computation,
/// then in parallel the user may specify additional range of indexes
/// for elements that lay on the borders between each pair of processors.
class RowMerger
{
public:
static const INMOST_DATA_ENUM_TYPE EOL = ENUMUNDEF-1; ///< End of linked list.
static const INMOST_DATA_ENUM_TYPE UNDEF = ENUMUNDEF; ///< Value not defined in linked list.
private:
bool Sorted; ///< Contents of linked list should be sorted.
INMOST_DATA_ENUM_TYPE Nonzeros; ///< Number of nonzero in linked list.
interval< INMOST_DATA_ENUM_TYPE, Row::entry > LinkedList; ///< Storage for linked list.
public:
/// Default constructor without size specfied.
RowMerger();
/// Constructor with size specified.
/// @param interval_begin First index in linked list.
/// @param interval_end Last index in linked list.
/// @param Sorted Result should be sorted.
RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
/// Constructor that gets sizes from matrix
/// @param A Matrix to get sizes from.
/// @param Sorted Result should be sorted.
RowMerger(Matrix & A, bool Sorted = true);
/// Destructor.
~RowMerger();
/// Resize linked list for new interval.
/// \warning
/// All contents of linked list will be lost after resize.
/// This behavior may be changed in future.
/// @param interval_begin First index in linked list.
/// @param interval_end Last index in linked list.
/// @param Sorted Result should be sorted.
void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
/// Resize linked list for new matrix.
/// \warning
/// All contents of linked list will be lost after resize.
/// This behavior may be changed in future.
/// @param A Matrix to get sizes from.
/// @param Sorted Result should be sorted.
void Resize(Matrix & A, bool Sorted = true);
/// Clear linked list.
void Clear();
/// Add a row with a coefficient into empty linked list.
/// This routing should be a bit faster then Solver::RowMerger::AddRow
/// for empty linked list. It may result in an unexpected behavior
/// for non-empty linked list, asserts will fire in debug mode.
/// @param coef Coefficient to multiply row values.
/// @param r A row to be added.
/// @param PreSortRow Sort values of the row before adding. Will be activated only for sorted linked lists.
void PushRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow = false);
/// Add a row with a coefficient into non-empty linked list.
/// Use Solver::RowMerger::PushRow for empty linked list.
/// @param coef Coefficient to multiply row values.
/// @param r A row to be added.
/// @param PreSortRow Sort values of the row before adding. Will be activated only for sorted linked lists.
void AddRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow = false);
/// Multiply all entries of linked list by a coefficient.
/// @param coef A coefficient for multiplication.
void Multiply(INMOST_DATA_REAL_TYPE coef);
/// Place entries from linked list into row.
/// \warning
/// All contents of the row will be overwritten.
/// If you want contents of the row to be added
/// use AddRow with this row in advance.
/// @param r A row to be filled.
void RetriveRow(Row & r);
//INMOST_DATA_REAL_TYPE ScalarProd(RowMerger & other);
/// Get current number of nonzeros from linked list.
INMOST_DATA_ENUM_TYPE Size() {return Nonzeros;}
};
private:
static bool is_initialized, is_finalized;
View
@@ -94,6 +94,126 @@ namespace INMOST
#endif
}
Solver::RowMerger::RowMerger() : Sorted(true), Nonzeros(0) {}
Solver::RowMerger::RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted)
: Sorted(Sorted), Nonzeros(0), LinkedList(interval_begin,interval_end+1,Row::make_entry(UNDEF,0.0))
{
LinkedList.begin()->first = EOL;
}
void Solver::RowMerger::Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool _Sorted)
{
LinkedList.set_interval_beg(interval_begin);
LinkedList.set_interval_end(interval_end+1);
std::fill(LinkedList.begin(),LinkedList.end(),Row::make_entry(UNDEF,0.0));
LinkedList.begin()->first = EOL;
Nonzeros = 0;
Sorted = _Sorted;
}
void Solver::RowMerger::Resize(Matrix & A, bool _Sorted)
{
INMOST_DATA_ENUM_TYPE mbeg, mend;
A.GetInterval(mbeg,mend);
Resize(mbeg,mend,_Sorted);
}
Solver::RowMerger::RowMerger(Matrix & A, bool Sorted) : Sorted(Sorted), Nonzeros(0)
{
INMOST_DATA_ENUM_TYPE mbeg, mend;
A.GetInterval(mbeg,mend);
LinkedList.set_interval_beg(mbeg);
LinkedList.set_interval_end(mend+1);
std::fill(LinkedList.begin(),LinkedList.end(),Row::make_entry(UNDEF,0.0));
LinkedList.begin()->first = EOL;
}
Solver::RowMerger::~RowMerger() {}
void Solver::RowMerger::Clear()
{
INMOST_DATA_ENUM_TYPE i = LinkedList.begin()->first, j;
LinkedList.begin()->first = EOL;
while( i != EOL )
{
j = LinkedList[i].first;
LinkedList[i].first = UNDEF;
i = j;
}
Nonzeros = 0;
}
void Solver::RowMerger::PushRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow)
{
if( Sorted && PreSortRow ) std::sort(r.Begin(),r.End());
assert(Nonzeros == 0); //Linked list should be empty
assert(LinkedList.begin()->first == EOL); //again check that list is empty
INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg();
Row::iterator it = r.Begin(), jt;
while( it != r.End() )
{
LinkedList[index].first = it->first+1;
LinkedList[it->first+1].first = EOL;
LinkedList[it->first+1].second = it->second*coef;
index = it->first+1;
++Nonzeros;
jt = it;
++it;
assert(!Sorted || it == r.End() || jt->first < it->first);
}
}
void Solver::RowMerger::AddRow(INMOST_DATA_REAL_TYPE coef, Row & r, bool PreSortRow)
{
if( Sorted && PreSortRow ) std::sort(r.Begin(),r.End());
INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg(), next;
Row::iterator it = r.Begin(), jt;
while( it != r.End() )
{
if( LinkedList[it->first+1].first != UNDEF )
LinkedList[it->first+1].second += coef*it->second;
else if( Sorted )
{
next = index;
while(next < it->first+1)
{
index = next;
next = LinkedList[index].first;
}
assert(index < it->first+1);
assert(it->first+1 < next);
LinkedList[index].first = it->first+1;
LinkedList[it->first+1].first = next;
LinkedList[it->first+1].second = coef*it->second;
++Nonzeros;
}
else
{
LinkedList[it->first+1].first = LinkedList[index].first;
LinkedList[it->first+1].second = coef*it->second;
LinkedList[index].first = it->first+1;
++Nonzeros;
}
jt = it;
++it;
assert(!Sorted || it == r.End() || jt->first < it->first);
}
}
void Solver::RowMerger::RetriveRow(Row & r)
{
r.Resize(Nonzeros);
INMOST_DATA_ENUM_TYPE i = LinkedList.begin()->first, k = 0;
while( i != EOL )
{
r.GetIndex(k) = i-1;
r.GetValue(k) = LinkedList[i].second;
i = LinkedList[i].first;
++k;
}
}
void Solver::OrderInfo::PrepareMatrix(Matrix & m, INMOST_DATA_ENUM_TYPE overlap)
{
have_matrix = true;
Oops, something went wrong.

0 comments on commit 82230d1

Please sign in to comment.