Skip to content

Commit

Permalink
FIX adding custom constrains matrix into lazy evaluation process of r…
Browse files Browse the repository at this point in the history 
…egularization settings (fixes #659)
  • Loading branch information
@carsten-forty2
carsten-forty2 committed Feb 28, 2024
1 parent 320f5fc commit d72b828
Show file tree
Hide file tree
Showing 17 changed files with 188 additions and 92 deletions.
52 changes: 37 additions & 15 deletions core/src/bert/dcfemmodelling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -556,14 +556,14 @@ void DCMultiElectrodeModelling::init_(){
JIsRMatrix_ = true;
JIsCMatrix_ = false;

solver_ = nullptr;
buildCompleteElectrodeModel_ = false;
dipoleCurrentPattern_ = false;

primDataMap_ = new DataMap();

byPassFile_ = "bypass.map";


Index nThreads = getEnvironment("BERTTHREADS", 0, verbose_);
nThreads = getEnvironment("BERT_NUM_THREADS", 0, verbose_);
if (nThreads > 0) setThreadCount(nThreads);
Expand Down Expand Up @@ -1838,14 +1838,18 @@ void DCMultiElectrodeModelling::calculateK_(const std::vector < ElectrodeShape *
//** END assemble matrix

//** START solving

LinSolver solver(debug);
//solver.setSolverType(LDL);
// std::cout << "solver: " << solver.solverName() << std::endl;

solver.setMatrix(S_, 1);
if (debug) std::cout << "Factorizing (" << solver.solverName() << ") system matrix ... ";

SolverWrapper * solver;
bool solverNeedsDelete = false;
if (solver_ != nullptr){
solver = solver_;
solver->setMatrix(S_);
} else {
solver = new LinSolver(debug);
solverNeedsDelete = true;
dynamic_cast< LinSolver * >(solver)->setMatrix(S_, 1);
if (debug) std::cout << "Factorize (" << solver->name() << ") matrix ... " << swatch.duration() << std::endl;
}

// MEMINFO

Vector < ValueType > sol(S_.cols());
Expand All @@ -1860,7 +1864,7 @@ void DCMultiElectrodeModelling::calculateK_(const std::vector < ElectrodeShape *
if (eB[i]) eB[i]->assembleRHS(rTmp, -1.0, oldMatSize);
Vector < ValueType >rhs(rTmp);

solver.solve(rhs, sol);
solver->solve(rhs, sol);

/*if (i==4){
__MS(*mesh_)
Expand Down Expand Up @@ -1893,7 +1897,7 @@ void DCMultiElectrodeModelling::calculateK_(const std::vector < ElectrodeShape *
// S_.save("S.matrix");
// save(rhs[i], "rhs.vec");
// save(sol, "sol.vec");
std::cout << " Ooops: Warning!!!! Solver: " << solver.solverName()
std::cout << " Ooops: Warning!!!! Solver: " << solver->name()
<< " fails with rms(A *x -b)/rms(b) > tol: "
<< norml2(S_ * sol - rhs)<< std::endl;

Expand All @@ -1918,6 +1922,9 @@ void DCMultiElectrodeModelling::calculateK_(const std::vector < ElectrodeShape *
// MEMINFO
// we dont need reserve the memory
S_.clean();
if (solverNeedsDelete == true){
delete solver;
}
}


Expand Down Expand Up @@ -2189,9 +2196,19 @@ void DCSRMultiElectrodeModelling::calculateK(const std::vector < ElectrodeShape
//mesh_->setCellAttributes(tmpRho);

// MEMINFO
LinSolver solver(debug);
solver.setMatrix(S_, 1);
if (debug) std::cout << "Factorize (" << solver.solverName() << ") matrix ... " << swatch.duration() << std::endl;
SolverWrapper *solver;
bool solverNeedsDelete = false;

if (solver_ != nullptr){
solver = solver_;
solver->setMatrix(S_);
} else {
solver = new LinSolver(debug);
solverNeedsDelete = true;
dynamic_cast< LinSolver * >(solver)->setMatrix(S_, 1);
if (debug) std::cout << "Factorize (" << solver->name() << ") matrix ... " << swatch.duration() << std::endl;
}


// MEMINFO

Expand Down Expand Up @@ -2265,12 +2282,17 @@ void DCSRMultiElectrodeModelling::calculateK(const std::vector < ElectrodeShape
rhs[calibrationSourceIdx_[j]] = 0.0;
}
solutionK[i + (kIdx * nCurrentPattern)] *= 0.0;
solver.solve(rhs, solutionK[i + (kIdx * nCurrentPattern)]);
solver->solve(rhs, solutionK[i + (kIdx * nCurrentPattern)]);

solutionK[i + (kIdx * nCurrentPattern)] += prim;

if (singleVerbose) singleVerbose = false;
}

if (solverNeedsDelete == true){
delete solver;
}

}

} // namespace GIMLI{
5 changes: 5 additions & 0 deletions core/src/bert/dcfemmodelling.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,6 +262,9 @@ class DLLEXPORT DCMultiElectrodeModelling : public GIMLI::ModellingBase {
/*! Return true if singular value estimation is switched on.*/
bool isSetSingValue() const { return setSingValue_;}

/*! Set a custom solver if you don't want the default Choldmod or UMFPACK. */
void setSolver(SolverWrapper *solver){ solver_ = solver; }

private:
void init_();

Expand Down Expand Up @@ -334,6 +337,8 @@ class DLLEXPORT DCMultiElectrodeModelling : public GIMLI::ModellingBase {
RVector vContactImpedance_;

DataMap * primDataMap_;

SolverWrapper *solver_;
};

class DLLEXPORT DCSRMultiElectrodeModelling : public DCMultiElectrodeModelling {
Expand Down
70 changes: 55 additions & 15 deletions core/src/cholmodWrapper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,17 +45,43 @@ namespace GIMLI{

CHOLMODWrapper::CHOLMODWrapper(RSparseMatrix & S, bool verbose, int stype,
bool forceUmfpack)
: SolverWrapper(S, verbose), forceUmfpack_(forceUmfpack){
init_(S, stype);
: SolverWrapper(verbose), stype_(stype), forceUmfpack_(forceUmfpack){

Numeric_ = nullptr;
NumericD_= nullptr;
c_ = NULL;
A_ = NULL;
L_ = NULL;
AxV_ = NULL;
AzV_ = NULL;
Ap_ = NULL;
Ai_ = NULL;
ApR_ = NULL;
AiR_ = NULL;
setMatrix(S);
}

CHOLMODWrapper::CHOLMODWrapper(CSparseMatrix & S, bool verbose, int stype,
bool forceUmfpack)
: SolverWrapper(S, verbose), forceUmfpack_(forceUmfpack){
init_(S, stype);
: SolverWrapper(verbose), stype_(stype), forceUmfpack_(forceUmfpack){
Numeric_ = nullptr;
NumericD_= nullptr;
c_ = NULL;
A_ = NULL;
L_ = NULL;
AxV_ = NULL;
AzV_ = NULL;
Ap_ = NULL;
Ai_ = NULL;
ApR_ = NULL;
AiR_ = NULL;
setMatrix(S);
}

CHOLMODWrapper::~CHOLMODWrapper(){
this->free();
}

void CHOLMODWrapper::free(){
#if USE_CHOLMOD
if (L_) cholmod_free_factor((cholmod_factor**)(&L_), (cholmod_common*)c_);
// ** We did not allocate the matrix so we dont need to free it
Expand All @@ -64,27 +90,48 @@ CHOLMODWrapper::~CHOLMODWrapper(){
cholmod_finish((cholmod_common*)c_);

if (A_) delete (cholmod_sparse*)A_;
A_ = nullptr;
if (c_) delete (cholmod_common*)c_;
c_ = nullptr;

#if USE_UMFPACK
if (Numeric_) umfpack_zi_free_numeric (&Numeric_);
Numeric_ = nullptr;
if (NumericD_) umfpack_di_free_numeric (&NumericD_);
NumericD_ = nullptr;
#endif
// try allways
if (AxV_) delete AxV_;
AxV_ = nullptr;
if (AzV_) delete AzV_;
AzV_ = nullptr;

if (ApR_) delete [] ApR_;
ApR_ = nullptr;

if (AiR_) delete [] AiR_;
AiR_ = nullptr;

#else
std::cerr << WHERE_AM_I << " cholmod not installed" << std::endl;
#endif
}

void CHOLMODWrapper::setMatrix(RSparseMatrix & S){
this->free();
init_(S, stype_);
}
void CHOLMODWrapper::setMatrix(CSparseMatrix & S){
this->free();
init_(S, stype_);
}

template < class ValueType >
void CHOLMODWrapper::init_(SparseMatrix < ValueType > & S, int stype){

dim_ = S.size();
nVals_ = S.nVals();

useUmfpack_ = false;
Numeric_ = 0;
NumericD_ = 0;
Expand Down Expand Up @@ -330,7 +377,7 @@ int CHOLMODWrapper::factorise_(){
}

template < class ValueType >
int CHOLMODWrapper::solveCHOL_(const Vector < ValueType > & rhs,
void CHOLMODWrapper::solveCHOL_(const Vector < ValueType > & rhs,
Vector < ValueType > & solution){
ASSERT_VEC_SIZE(rhs, this->dim_)
ASSERT_VEC_SIZE(solution, this->dim_)
Expand Down Expand Up @@ -368,15 +415,13 @@ template < class ValueType >
cholmod_free_dense(&x, (cholmod_common*)c_);
cholmod_free_dense(&r, (cholmod_common*)c_);
cholmod_free_dense(&b, (cholmod_common*)c_);
return 1;
#else
std::cerr << WHERE_AM_I << " cholmod not installed" << std::endl;
#endif
}
return 0;
}

int CHOLMODWrapper::solve(const RVector & rhs, RVector & solution){
void CHOLMODWrapper::solve(const RVector & rhs, RVector & solution){
ASSERT_VEC_SIZE(rhs, this->dim_)
ASSERT_VEC_SIZE(solution, this->dim_)
if (!dummy_){
Expand Down Expand Up @@ -406,17 +451,14 @@ int CHOLMODWrapper::solve(const RVector & rhs, RVector & solution){
// umfpack_di_free_numeric (&Numeric) ;
//


return 1;
#endif
} else {
return solveCHOL_(rhs, solution);
}
}
return 0;
}

int CHOLMODWrapper::solve(const CVector & rhs, CVector & solution){
void CHOLMODWrapper::solve(const CVector & rhs, CVector & solution){
ASSERT_VEC_SIZE(rhs, this->dim_)
ASSERT_VEC_SIZE(solution, this->dim_)
if (!dummy_){
Expand All @@ -442,13 +484,11 @@ int CHOLMODWrapper::solve(const CVector & rhs, CVector & solution){
// set booth to NULL or they will be wrongly deleted ..

solution = toComplex(xx, xz);
return 1;
#endif
} else {
return solveCHOL_(rhs, solution);
}
}
return 0;
}


Expand Down
21 changes: 12 additions & 9 deletions core/src/cholmodWrapper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,23 +26,26 @@ namespace GIMLI{

class DLLEXPORT CHOLMODWrapper : public SolverWrapper {
public:
CHOLMODWrapper(RSparseMatrix & S, bool verbose=false, int stype=-2,
bool forceUmfpack=false);

CHOLMODWrapper(CSparseMatrix & S, bool verbose=false, int stype=-2,
bool forceUmfpack=false);
CHOLMODWrapper(RSparseMatrix & S, bool verbose=false, int stype=-2, bool forceUmfpack=false);
CHOLMODWrapper(CSparseMatrix & S, bool verbose=false, int stype=-2, bool forceUmfpack=false);

virtual ~CHOLMODWrapper();

static bool valid();

virtual int solve(const RVector & rhs, RVector & solution);
virtual void setMatrix(RSparseMatrix & S);

virtual void setMatrix(CSparseMatrix & S);

virtual int solve(const CVector & rhs, CVector & solution);
virtual void solve(const RVector & rhs, RVector & solution);

virtual void solve(const CVector & rhs, CVector & solution);

protected:
void init();

void free();

int initializeMatrix_(RSparseMatrix & S);

int initializeMatrix_(CSparseMatrix & S);
Expand All @@ -54,10 +57,10 @@ class DLLEXPORT CHOLMODWrapper : public SolverWrapper {
int initMatrixChol_(SparseMatrix < ValueType > & S, int xType);

template < class ValueType >
int solveCHOL_(const Vector < ValueType > & rhs, Vector < ValueType > & solution);
void solveCHOL_(const Vector < ValueType > & rhs, Vector < ValueType > & solution);

template < class ValueType >
int solveUmf_(const Vector < ValueType > & rhs, Vector < ValueType > & solution);
void solveUmf_(const Vector < ValueType > & rhs, Vector < ValueType > & solution);

int factorise_();

Expand Down
1 change: 1 addition & 0 deletions core/src/gimli.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,6 +304,7 @@ class RegionManager;
class Shape;
class Stopwatch;
class Pos;
class SolverWrapper;

typedef Pos RVector3;
typedef std::complex < double > Complex;
Expand Down
Loading

0 comments on commit d72b828

Please sign in to comment.