mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2024-11-14 09:45:36 +00:00
146 lines
4.5 KiB
C++
146 lines
4.5 KiB
C++
/*
|
|
* GslHybridRootFinder.cpp, part of LatAnalyze 3
|
|
*
|
|
* Copyright (C) 2013 - 2020 Antonin Portelli
|
|
*
|
|
* LatAnalyze 3 is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* LatAnalyze 3 is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include <LatAnalyze/Numerical/GslHybridRootFinder.hpp>
|
|
#include <LatAnalyze/includes.hpp>
|
|
#include <gsl/gsl_vector.h>
|
|
#include <gsl/gsl_multiroots.h>
|
|
|
|
using namespace std;
|
|
using namespace Latan;
|
|
|
|
/******************************************************************************
|
|
* GslHybridRootFinder implementation *
|
|
******************************************************************************/
|
|
// constructor /////////////////////////////////////////////////////////////////
|
|
GslHybridRootFinder::GslHybridRootFinder(const Index dim)
|
|
: RootFinder(dim)
|
|
{}
|
|
|
|
// output //////////////////////////////////////////////////////////////////////
|
|
void GslHybridRootFinder::printState(void)
|
|
{
|
|
if (solver_)
|
|
{
|
|
cout << "x=";
|
|
for (size_t i = 0; i < solver_->x->size; ++i)
|
|
{
|
|
cout << " " << scientific << gsl_vector_get(solver_->x, i);
|
|
}
|
|
cout << endl;
|
|
cout << "f=";
|
|
for (size_t i = 0; i < solver_->f->size; ++i)
|
|
{
|
|
cout << " " << scientific << gsl_vector_get(solver_->f, i);
|
|
}
|
|
cout << endl;
|
|
}
|
|
}
|
|
|
|
// solver //////////////////////////////////////////////////////////////////////
|
|
const DVec &
|
|
GslHybridRootFinder::operator()(const vector<DoubleFunction *> &func)
|
|
{
|
|
DVec &res = getState();
|
|
Verbosity verbosity = getVerbosity();
|
|
int status;
|
|
unsigned int iter = 0;
|
|
const size_t nFunc = func.size();
|
|
Index nArg;
|
|
gsl_vector *x;
|
|
gsl_multiroot_function fStruct;
|
|
|
|
int (*fWrap)(const gsl_vector *, void *, gsl_vector *) =
|
|
[](const gsl_vector *var, void *vFunc, gsl_vector *f)->int
|
|
{
|
|
vector<DoubleFunction *> &fPt =
|
|
*static_cast<vector<DoubleFunction *> *>(vFunc);
|
|
|
|
for (unsigned int i = 0; i < fPt.size(); ++i)
|
|
{
|
|
gsl_vector_set(f, i, (*fPt[i])(var->data));
|
|
}
|
|
|
|
return GSL_SUCCESS;
|
|
};
|
|
|
|
nArg = func[0]->getNArg();
|
|
for (auto f: func)
|
|
{
|
|
if (f->getNArg() != nArg)
|
|
{
|
|
LATAN_ERROR(Size,
|
|
"equations do not have the same number of unknown");
|
|
}
|
|
}
|
|
if (nArg != static_cast<Index>(nFunc))
|
|
{
|
|
LATAN_ERROR(Size, "equation and unknown number mismatch");
|
|
}
|
|
if (res.size() != nArg)
|
|
{
|
|
resize(nArg);
|
|
}
|
|
solver_ = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrids, nFunc);
|
|
x = gsl_vector_alloc(nFunc);
|
|
FOR_VEC(res, i)
|
|
{
|
|
gsl_vector_set(x, static_cast<size_t>(i), res(i));
|
|
}
|
|
fStruct.n = nFunc;
|
|
fStruct.params = reinterpret_cast<void *>(
|
|
const_cast<vector<DoubleFunction *> *>(&func));
|
|
fStruct.f = fWrap;
|
|
gsl_multiroot_fsolver_set(solver_, &fStruct, x);
|
|
do
|
|
{
|
|
iter++;
|
|
status = gsl_multiroot_fsolver_iterate(solver_);
|
|
if (verbosity >= Verbosity::Debug)
|
|
{
|
|
cout << "--- iteration " << iter << endl;
|
|
printState();
|
|
}
|
|
if (status)
|
|
{
|
|
break;
|
|
}
|
|
status = gsl_multiroot_test_residual(solver_->f, getPrecision());
|
|
} while ((status == GSL_CONTINUE) and (iter < getMaxIteration()));
|
|
if (verbosity >= Verbosity::Debug)
|
|
{
|
|
cout << "--- done" << endl;
|
|
cout << "end status: " << gsl_strerror(status) << endl;
|
|
}
|
|
if (status)
|
|
{
|
|
LATAN_WARNING("GSL hybrid root finder ended with status '" +
|
|
strFrom(gsl_strerror(status)) + "'");
|
|
}
|
|
FOR_VEC(res, i)
|
|
{
|
|
res(i) = gsl_vector_get(solver_->x, static_cast<size_t>(i));
|
|
}
|
|
gsl_vector_free(x);
|
|
gsl_multiroot_fsolver_free(solver_);
|
|
solver_ = nullptr;
|
|
|
|
return res;
|
|
}
|