newton_refine.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//bertini2/nag_algorithms/newton_refine.hpp 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.
//
//bertini2/nag_algorithms/newton_refine.hpp 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 bertini2/nag_algorithms/newton_refine.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// silviana amethyst, university of wisconsin eau claire

 
#pragma once
 
#include "bertini2/system/system.hpp"
#include "bertini2/common/config.hpp"
#include "bertini2/linalg/lu_solver.hpp"
 
namespace bertini {
namespace algorithm {

template <typename ComplexT>
struct NewtonRefineResult
{
    SuccessCode code;  
    Vec<ComplexT> point;         
    NumErrorT achieved;          
    unsigned iterations;         
};

template <typename ComplexT>
NewtonRefineResult<ComplexT> NewtonRefine(System const& S,
                                          Vec<ComplexT> const& start,
                                          NumErrorT tolerance,
                                          unsigned max_iterations,
                                          ComplexT time = ComplexT(0))
{
    const auto n_funcs = S.NumTotalFunctions();
    const auto n_vars = S.NumVariables();
    if (n_funcs != n_vars)
    {
        std::stringstream ss;
        ss << "NewtonRefine requires a SQUARE system, but this one has "
           << n_funcs << " total functions (user functions plus patches) over "
           << n_vars << " variables.  Square an overdetermined system by "
           << "randomization (multiply by a generic full-rank matrix) before refining.";
        throw std::runtime_error(ss.str());
    }
    if (start.size() != static_cast<Eigen::Index>(n_vars))
    {
        std::stringstream ss;
        ss << "NewtonRefine: start point has " << start.size()
           << " coordinates but the system has " << n_vars << " variables.";
        throw std::runtime_error(ss.str());
    }
 
    NewtonRefineResult<ComplexT> result{SuccessCode::FailedToConverge,
                                        start, static_cast<NumErrorT>(-1), 0};
 
    Vec<ComplexT> f(n_funcs);
    Mat<ComplexT> J(n_funcs, n_vars);
    Vec<ComplexT> step(n_vars);
    linalg::PartialPivLU<ComplexT> lu;
    lu.ChangeSize(static_cast<unsigned>(n_vars));   // the LU workspace is stateful:
                                                    // without this, Factor/Solve run
                                                    // over a 0-dimensional system
 
    for (unsigned it = 0; it < max_iterations; ++it)
    {
        if (S.HavePathVariable())
            S.SetAndReset(result.point, time);
        else
            S.SetAndReset(result.point);
        S.EvalInPlace(f);
        S.JacobianInPlace(J);
 
        if (lu.Factor(J) != MatrixSuccessCode::Success)
        {
            result.code = SuccessCode::MatrixSolveFailure;
            return result;
        }
        lu.Solve(f, step);          // step = +J^{-1} f = -(Newton step)
        result.point -= step;
        result.iterations = it + 1;
        result.achieved = static_cast<NumErrorT>(
            step.template lpNorm<Eigen::Infinity>());
 
        if (result.achieved <= tolerance)
        {
            result.code = SuccessCode::Success;
            return result;
        }
    }
    return result;
}
 
} // namespace algorithm
} // namespace bertini