Commit 2e23260c authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Add true residual calculation method

parent 8c1f10fe
......@@ -44,7 +44,7 @@
// in Mesh::Init function change two variables:
// check_shared_mrk - check shared markers.
// check_private_mrk - check private markers.
#define CHECKS_MARKERS
//#define CHECKS_MARKERS
// use additional sets to store elements for parallel
// exchanges, otherwise it will recompute those elements
......
......@@ -261,7 +261,7 @@ namespace INMOST
/// and the preconditioner have been already constructed.
///
/// @see Sparse::SetMatrix
bool Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL);
bool Solve(Sparse::Vector &RHS, Sparse::Vector &SOL);
/// Clear all internal data of the current solver including matrix, preconditioner etc.
bool Clear();
/// Get the solver output parameter
......@@ -342,7 +342,7 @@ namespace INMOST
/// Return the number of iterations performed by the last solution.
/// @see Sparse::Solve
INMOST_DATA_ENUM_TYPE Iterations() const;
/// Return the final residual achieved by the last solution.
/// Return the final precondioned residual achieved by the last solution.
/// @see Sparse::Solve
INMOST_DATA_REAL_TYPE Residual() const;
/// Get the reason of convergence or divergence of the last solution.
......
......@@ -521,6 +521,15 @@ namespace INMOST
const bool & isParallel() const { return is_parallel; }
/// Get the matrix name specified in the main constructor.
std::string GetName() const {return name;}
/// Calculate the real residual.
///
/// @param RHS The right-hand side Vector b.
/// @param SOL The initial guess to the solution on input and the solution Vector x on return.
/// @return ||A*x-b||
///
/// It is assumed that the coefficient matrix A have been set
/// and the preconditioner have been already constructed.
INMOST_DATA_REAL_TYPE Residual(Sparse::Vector &RHS, Sparse::Vector &SOL);
};
#endif //defined(USE_SOLVER)
......
......@@ -220,7 +220,7 @@ namespace INMOST
}
*/
void OperationMinDistance(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
static void OperationMinDistance(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
{
(void)size;
assert(size == 2);
......@@ -237,6 +237,15 @@ namespace INMOST
rdata[1] = idata[1];
}
}
static void OperationMinOwner(const Tag & tag, const Element & element, const INMOST_DATA_BULK_TYPE * data, INMOST_DATA_ENUM_TYPE size)
{
(void)size;
assert(size == 1);
INMOST_DATA_INTEGER_TYPE idata = *(INMOST_DATA_INTEGER_TYPE *)data;
element->Integer(tag) = std::min(element->Integer(tag),idata);
}
void Mesh::EquilibrateGhost(bool only_new)
{
if( GetProcessorsNumber() == 1 ) return;
......@@ -244,7 +253,7 @@ namespace INMOST
#if defined(USE_MPI)
ElementType bridge = FACE;
if( tag_bridge.isValid() ) bridge = ElementType(Integer(GetHandle(),tag_bridge));
TagIntegerArray tag = CreateTag("TEMPORARY_OWNER_DISTANCE_TAG",DATA_INTEGER,CELL,NONE,2);
TagIntegerArray tag = CreateTag("TEMPORARY_OWNER_DISTANCE_TAG",DATA_INTEGER,CELL,CELL,2);
std::queue<HandleType> cells_queue;
ENTER_BLOCK();
//TODO: compute graph distance only in new elements!!!!
......@@ -312,7 +321,6 @@ namespace INMOST
if( !c.Hidden() && c.GetStatus() != Element::Owned )
{
//if ( !only_new || it->nbAdjElements(NODE | EDGE | FACE | CELL, NewMarker()) != 0)
//if( !only_new || c.nbAdjElements(NODE | FACE | CELL, NewMarker()) )
{
int new_owner = tag[c][0];
c.IntegerDF(tag_owner) = new_owner;
......@@ -324,31 +332,46 @@ namespace INMOST
}
}
EXIT_BLOCK();
DeleteTag(tag);
ENTER_BLOCK();
TagInteger new_owner = CreateTag("TEMPORARY_NEW_OWNER",DATA_INTEGER,NODE|EDGE|FACE,NODE|EDGE|FACE,1);
for(ElementType etype = FACE; etype >= NODE; etype = PrevElementType(etype))
{
//ElementType check = etype | NextElementType(etype) | PrevElementType(etype);
//if( etype == NODE ) etype |= CELL;
ElementType check = NODE | FACE | EDGE | CELL;
for(int k = 0; k < LastLocalID(etype); k++) if( isValidElement(etype,k) )
{
Element e = ElementByLocalID(etype,k);
if( !e.Hidden() && e.GetStatus() != Element::Owned )
{
//if( !only_new || e.nbAdjElements(check,NewMarker()) )
//if( !only_new || e.nbAdjElements(NODE | EDGE | FACE | CELL,NewMarker()) )
{
Element::adj_type & hc = HighConn(e.GetHandle());
int new_owner = INT_MAX;
new_owner[e] = INT_MAX;
for(int l = 0; l < hc.size(); ++l)
{
if( !GetMarker(hc[l],HideMarker()) )
new_owner = std::min(new_owner,Integer(hc[l],tag_owner));
new_owner[e] = std::min(new_owner[e],Integer(hc[l],tag_owner));
}
if( new_owner != INT_MAX &&
e.IntegerDF(tag_owner) != new_owner )
}
}
}
ReduceData(new_owner,etype,0,OperationMinOwner);
ExchangeData(new_owner,etype);
for(int k = 0; k < LastLocalID(etype); k++) if( isValidElement(etype,k) )
{
Element e = ElementByLocalID(etype,k);
if( !e.Hidden() && e.GetStatus() != Element::Owned )
{
//if( !only_new || e.nbAdjElements(NODE | EDGE | FACE | CELL,NewMarker()) )
{
if( new_owner[e] != INT_MAX &&
e.IntegerDF(tag_owner) != new_owner[e] )
{
e.IntegerDF(tag_owner) = new_owner;
if (GetProcessorRank() == new_owner)
e.IntegerDF(tag_owner) = new_owner[e];
if (GetProcessorRank() == new_owner[e])
e.SetStatus(Element::Shared);
else
e.SetStatus(Element::Ghost);
......@@ -357,10 +380,8 @@ namespace INMOST
}
}
}
DeleteTag(new_owner);
DeleteTag(tag);
EXIT_BLOCK();
RecomputeParallelStorage(CELL | EDGE | FACE | NODE);
#endif
......
......@@ -3,261 +3,365 @@
#include "../Misc/utils.h"
#if defined(USE_SOLVER)
namespace INMOST {
int *Solver::argc = NULL;
char ***Solver::argv = NULL;
bool Solver::is_initialized = false;
bool Solver::is_finalized = false;
std::vector<SolverParameters> Solver::parameters = std::vector<SolverParameters>();
namespace INMOST
{
int *Solver::argc = NULL;
char ***Solver::argv = NULL;
bool Solver::is_initialized = false;
bool Solver::is_finalized = false;
std::vector<SolverParameters> Solver::parameters = std::vector<SolverParameters>();
void Solver::SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value) {SetParameter(name,to_string(value));}
void Solver::SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value) {SetParameter(name,to_string(value));}
Solver::Solver(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm) {
std::string lowerName = string_to_lower(solverName);
this->solver = SolverMaster::getSolver(solverName);
this->prefix = string_to_lower(prefix);
solver->SetCommunicator(_comm);
//TODO find easiest way
bool parametersFound = false;
if (Solver::parameters.size() > 0) {
for (solver_parameters_iterator_t parameters = Solver::parameters.end() - 1; parameters >= Solver::parameters.begin(); parameters--) {
if ((*parameters).solverName == lowerName && (*parameters).solverPrefix == (this->prefix)) {
solver->Setup(argc, argv, *parameters);
parametersFound = true;
break;
}
}
}
if (!parametersFound) {
SolverParameters emptyParameters(lowerName, this->prefix, "");
solver->Setup(argc, argv, emptyParameters);
}
}
Solver::Solver(const Solver &other) {
this->solver = SolverMaster::getSolver(other.solver->SolverName());
this->solver->Copy(other.solver);
this->prefix = other.prefix;
solver->SetCommunicator(other.solver->GetCommunicator());
//TODO find easiest way
bool parametersFound = false;
if (Solver::parameters.size() > 0) {
for (solver_parameters_iterator_t parameters = Solver::parameters.end() - 1; parameters >= Solver::parameters.begin(); parameters--) {
if ((*parameters).solverName == other.solver->SolverName() && (*parameters).solverPrefix == (this->prefix)) {
solver->Setup(argc, argv, *parameters);
parametersFound = true;
break;
}
}
}
if (!parametersFound) {
SolverParameters emptyParameters(other.solver->SolverName(), this->prefix, "");
solver->Setup(argc, argv, emptyParameters);
}
}
Solver &Solver::operator=(const Solver &other) {
if (this != &other) {
this->solver->SetCommunicator(other.solver->GetCommunicator());
this->prefix = other.prefix;
this->solver->Assign(other.solver);
}
return *this;
}
std::string Solver::SolverName() const {
return solver->SolverName();
}
std::string Solver::SolverPrefix() const {
return prefix;
}
void Solver::Initialize(int *argc, char ***argv, const char *database) {
Solver::argc = argc;
Solver::argv = argv;
Solver::is_initialized = true;
Solver::is_finalized = false;
Solver::Solver(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm)
{
std::string lowerName = string_to_lower(solverName);
this->solver = SolverMaster::getSolver(solverName);
this->prefix = string_to_lower(prefix);
solver->SetCommunicator(_comm);
//TODO find easiest way
bool parametersFound = false;
if (Solver::parameters.size() > 0)
{
for (solver_parameters_iterator_t parameters = Solver::parameters.end() - 1; parameters >= Solver::parameters.begin(); parameters--)
{
if ((*parameters).solverName == lowerName && (*parameters).solverPrefix == (this->prefix))
{
solver->Setup(argc, argv, *parameters);
parametersFound = true;
break;
}
}
}
if (!parametersFound)
{
SolverParameters emptyParameters(lowerName, this->prefix, "");
solver->Setup(argc, argv, emptyParameters);
}
}
Solver::Solver(const Solver &other)
{
this->solver = SolverMaster::getSolver(other.solver->SolverName());
this->solver->Copy(other.solver);
this->prefix = other.prefix;
solver->SetCommunicator(other.solver->GetCommunicator());
//TODO find easiest way
bool parametersFound = false;
if (Solver::parameters.size() > 0)
{
for (solver_parameters_iterator_t parameters = Solver::parameters.end() - 1; parameters >= Solver::parameters.begin(); parameters--)
{
if ((*parameters).solverName == other.solver->SolverName() && (*parameters).solverPrefix == (this->prefix))
{
solver->Setup(argc, argv, *parameters);
parametersFound = true;
break;
}
}
}
if (!parametersFound)
{
SolverParameters emptyParameters(other.solver->SolverName(), this->prefix, "");
solver->Setup(argc, argv, emptyParameters);
}
}
Solver &Solver::operator=(const Solver &other)
{
if (this != &other)
{
this->solver->SetCommunicator(other.solver->GetCommunicator());
this->prefix = other.prefix;
this->solver->Assign(other.solver);
}
return *this;
}
std::string Solver::SolverName() const
{
return solver->SolverName();
}
std::string Solver::SolverPrefix() const
{
return prefix;
}
void Solver::Initialize(int *argc, char ***argv, const char *database)
{
Solver::argc = argc;
Solver::argv = argv;
Solver::is_initialized = true;
Solver::is_finalized = false;
#if defined(USE_MPI)
{
int flag = 0;
int ierr = 0;
MPI_Initialized(&flag);
if (!flag) {
ierr = MPI_Init(argc, argv);
if (ierr != MPI_SUCCESS) {
std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Init" << std::endl;
}
}
}
{
int flag = 0;
int ierr = 0;
MPI_Initialized(&flag);
if (!flag)
{
ierr = MPI_Init(argc, argv);
if (ierr != MPI_SUCCESS)
{
std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Init" << std::endl;
}
}
}
#endif
Solver::parseXMLDatabase(database);
Sparse::CreateRowEntryType();
}
void Solver::Finalize() {
Sparse::DestroyRowEntryType();
Solver::parseXMLDatabase(database);
Sparse::CreateRowEntryType();
}
void Solver::Finalize()
{
Sparse::DestroyRowEntryType();
#if defined(USE_MPI)
{
int flag = 0;
MPI_Finalized(&flag);
if (!flag) {
MPI_Finalize();
}
}
{
int flag = 0;
MPI_Finalized(&flag);
if (!flag)
MPI_Finalize();
}
#endif
Solver::is_finalized = true;
Solver::is_initialized = false;
}
bool Solver::isInitialized() {
return is_initialized;
}
bool Solver::isFinalized() {
return is_finalized;
}
void Solver::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
solver->SetMatrix(A, ModifiedPattern, OldPreconditioner);
}
bool Solver::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
if (!solver->isMatrixSet()) throw MatrixNotSetInSolver;
if (RHS.GetCommunicator() != solver->GetCommunicator() ||
SOL.GetCommunicator() != solver->GetCommunicator())
throw DifferentCommunicatorInSolver;
INMOST_DATA_ENUM_TYPE vbeg, vend;
RHS.GetInterval(vbeg, vend);
if (RHS.Size() != SOL.Size()) {
if (SOL.Size() == 0) {
SOL.SetInterval(vbeg, vend);
for (Sparse::Vector::iterator ri = SOL.Begin(); ri != SOL.End(); ++ri) *ri = 0.0;
} else throw InconsistentSizesInSolver;
}
return solver->Solve(RHS, SOL);
}
bool Solver::Clear() {
return solver->Clear();
}
std::string Solver::GetParameter(std::string name) const {
return solver->GetParameter(name);
}
void Solver::SetParameter(std::string name, std::string value) {
return solver->SetParameter(name, value);
}
INMOST_DATA_ENUM_TYPE Solver::Iterations() const {
return solver->Iterations();
}
INMOST_DATA_REAL_TYPE Solver::Residual() const {
return solver->Residual();
}
const std::string Solver::ReturnReason() const {
return solver->ReturnReason();
}
INMOST_DATA_REAL_TYPE Solver::Condest(INMOST_DATA_REAL_TYPE tol, unsigned int maxits) {
if (!solver->isMatrixSet()) throw MatrixNotSetInSolver;
return solver->Condest(tol, maxits);
}
bool Solver::isSolverAvailable(std::string name) {
return SolverMaster::isSolverAvailable(name);
}
std::vector<std::string> Solver::getAvailableSolvers() {
return SolverMaster::getAvailableSolvers();
}
Solver::~Solver() {
solver->Finalize();
delete solver;
}
void Solver::parseXMLDatabase(const char *xml_database) {
if (xml_database == NULL || xml_database[0] == '\0') return;
std::ifstream input;
input.open(xml_database);
if (input.fail()) {
std::cout << __FILE__ << ":" << __LINE__ << ": XML database file " << xml_database << " not found." << std::endl;
return;
}
XMLReader reader(std::string(xml_database), input);
try {
XMLReader::XMLTree root = reader.ReadXML();
if (root.tag.name != "SolverParameters") {
std::cout << __FILE__ << ":" << __LINE__ << ": Bad XML database file " << xml_database << "!" << std::endl;
return;
}
for (xml_reader_tree_iterator_t solver = root.children.begin(); solver < root.children.end(); solver++) {
std::string solverName = string_to_lower((*solver).tag.name);
std::string internalFile = "";
if ((*solver).tag.attributes.size() != 0) {
for (xml_reader_attrib_iterator_t attr = (*solver).tag.attributes.begin(); attr < (*solver).tag.attributes.end(); attr++) {
if ((*attr).name == "File" || (*attr).name == "file" || (*attr).name == "FILE") {
internalFile = (*attr).value;
}
}
}
if ((*solver).children.size() == 0) {
//Internal solver
parameters.push_back(SolverParameters(solverName, "", internalFile));
} else {
//Inner solver
for (xml_reader_tree_iterator_t prefix = (*solver).children.begin(); prefix < (*solver).children.end(); prefix++) {
internalFile = "";
std::string solverPrefix = string_to_lower((*prefix).tag.name);
if ((*prefix).tag.attributes.size() != 0) {
for (xml_reader_attrib_iterator_t attr = (*prefix).tag.attributes.begin(); attr < (*prefix).tag.attributes.end(); attr++) {
if ((*attr).name == "File" || (*attr).name == "file" || (*attr).name == "FILE") {
internalFile = (*attr).value;
}
}
}
SolverParameters prefix_p = SolverParameters(solverName, solverPrefix, internalFile);
for (xml_reader_tree_iterator_t p = (*prefix).children.begin(); p < (*prefix).children.end(); p++) {
if ((*p).tag.attributes.size() == 1) {
if ((*p).tag.attributes[0].name == "value" || (*p).tag.attributes[0].name == "Value") {
prefix_p.parameters.push_back(std::make_pair((*p).tag.name, (*p).tag.attributes[0].value));
}
}
}
parameters.push_back(prefix_p);
}
}
}
} catch (...) {
std::cout << __FILE__ << ": Error while parsing xml database file" << std::endl;
}
}
Solver::is_finalized = true;
Solver::is_initialized = false;
}
bool Solver::isInitialized()
{
return is_initialized;
}
bool Solver::isFinalized()
{
return is_finalized;
}
void Solver::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner)
{
solver->SetMatrix(A, ModifiedPattern, OldPreconditioner);
}
bool Solver::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL)
{
if (!solver->isMatrixSet()) throw MatrixNotSetInSolver;
if (RHS.GetCommunicator() != solver->GetCommunicator() ||
SOL.GetCommunicator() != solver->GetCommunicator())
throw DifferentCommunicatorInSolver;
INMOST_DATA_ENUM_TYPE vbeg, vend;
RHS.GetInterval(vbeg, vend);
if (RHS.Size() != SOL.Size())
{
if (SOL.Size() == 0)
{
SOL.SetInterval(vbeg, vend);
for (Sparse::Vector::iterator ri = SOL.Begin(); ri != SOL.End(); ++ri) *ri = 0.0;
}
else throw InconsistentSizesInSolver;
}
return solver->Solve(RHS, SOL);
}
bool Solver::Clear()
{
return solver->Clear();
}
std::string Solver::GetParameter(std::string name) const
{
return solver->GetParameter(name);
}
void Solver::SetParameter(std::string name, std::string value)
{
return solver->SetParameter(name, value);
}
INMOST_DATA_ENUM_TYPE Solver::Iterations() const
{
return solver->Iterations();
}
INMOST_DATA_REAL_TYPE Solver::Residual() const
{
return solver->Residual();
}
const std::string Solver::ReturnReason() const
{
return solver->ReturnReason();
}
INMOST_DATA_REAL_TYPE Solver::Condest(INMOST_DATA_REAL_TYPE tol, unsigned int maxits)
{
if (!solver->isMatrixSet()) throw MatrixNotSetInSolver;
return solver->Condest(tol, maxits);
}