Homec4science
Diffusion SpecMiCP / ReactMiCP (oxydoreduc) Browse

specmicp/src/micpsolver/ncp_function.hpp
7faccfcacb59oxydoreduc

Last Change
Raw File

ncp_function.hpp
View Options

/*-------------------------------------------------------
- Module : micpsolver
- File : ncp_function.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_MICPSOLVE_NCPFUNCTION_HPP
#define SPECMICP_MICPSOLVE_NCPFUNCTION_HPP
#include <cassert>
//! \file ncp_function.hpp implements the ncp function
namespace specmicp {
namespace micpsolver {
//! \brief The Fisher-Burmeister NCP-function
//!
//! @param a a scalar, first argument of the NCP-function
//! @param b a scalar, second argument of the NCP-function
//!
//! References:
//! - \cite Munson2001
//! - \cite Facchinei2003
template <typename ScalarT>
ScalarT fisher_burmeister(ScalarT a, ScalarT b)
{
ScalarT s = std::abs(a) + std::abs(b);
if (s != 0)
{
s = s*std::sqrt((a*a)/(s*s) + (b*b)/(s*s));
}
if ( a + b <= 0)
{
return s - (a + b);
}
else
{
return -2*a*b/(s + (a+b));
}
}
//! \brief The penalized Fisher-Burmeister NCP-function
//!
//! @param t in (0, 1), penalization factor
//! @param a a scalar, first argument of the NCP-function
//! @param b a scalar, second argument of the NCP-function
//!
//! References:
//! - \cite Chen1997a
//! - \cite Chen2000
//! - \cite Munson2001
template <typename ScalarT>
inline ScalarT penalized_fisher_burmeister(ScalarT a, ScalarT b, ScalarT t)
{
assert(t >= 0);
assert(t <= 1);
return t*fisher_burmeister(a,b)-(1-t)*std::max(0.0, a)*std::max(0.0, b);
}
} // end namespace micpsolver
} // end namespace specmicp
#endif // SPECMICP_MICPSOLVE_NCPFUNCTION_HPP
c4science · Help