newton_correct.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//newton_correct.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.
//
//newton_correct.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 newton_correct.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) 2015 - 2021 by Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license, 
// as well as COPYING.  Bertini2 is provided with permitted 
// additional terms in the b2/licenses/ directory.
 
// individual authors of this file include:
// silviana amethyst, university of wisconsin eau claire
 

 
 
#ifndef BERTINI_NEWTON_CORRECT_HPP
#define BERTINI_NEWTON_CORRECT_HPP
 
 
#include "bertini2/trackers/amp_criteria.hpp"
#include "bertini2/trackers/config.hpp"
#include "bertini2/system.hpp"
 
namespace bertini{
    namespace tracking{
        namespace correct{
 
            

            template <typename ComplexType, typename RealType>
            SuccessCode NewtonLoop(Vec<ComplexType> & next_space,
                                   System const& S,
                                   Vec<ComplexType> const& current_space, // pass by value to get a copy of it
                                   ComplexType const& current_time, 
                                   RealType const& tracking_tolerance,
                                   unsigned min_num_newton_iterations,
                                   unsigned max_num_newton_iterations)
            {
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(max_num_newton_iterations >= min_num_newton_iterations && "max number newton iterations must be at least the min.");
                #endif
 
 
                next_space = current_space;
                for (unsigned ii = 0; ii < max_num_newton_iterations; ++ii)
                {
                    auto f = S.Eval(next_space, current_time);
                    auto J = S.Jacobian(next_space, current_time);
                    auto LU = J.lu();
 
                    if (LUPartialPivotDecompositionSuccessful(LU.matrixLU())!=MatrixSuccessCode::Success)
                        return SuccessCode::MatrixSolveFailure;
 
                    auto delta_z = LU.solve(-f);
                    next_space += delta_z;
 
                    if ( (delta_z.norm() < tracking_tolerance) && (ii >= (min_num_newton_iterations-1)) )
                        return SuccessCode::Success;
                }
 
                return SuccessCode::FailedToConverge;
            }
            
 
 

            template <typename ComplexType, typename RealType>
            SuccessCode NewtonLoop(Vec<ComplexType> & next_space,
                                   System const& S,
                                   Vec<ComplexType> const& current_space, // pass by value to get a copy of it
                                   ComplexType const& current_time, 
                                   RealType const& tracking_tolerance,
                                   unsigned min_num_newton_iterations,
                                   unsigned max_num_newton_iterations,
                                   AdaptiveMultiplePrecisionConfig const& AMP_config)
            {
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(max_num_newton_iterations >= min_num_newton_iterations && "max number newton iterations must be at least the min.");
                #endif
 
 
                next_space = current_space;
                for (unsigned ii = 0; ii < max_num_newton_iterations; ++ii)
                {
                    //TODO: wrap these into a single line.
                    auto f = S.Eval(next_space, current_time);
                    auto J = S.Jacobian(next_space, current_time);
                    auto LU = J.lu();
 
                    if (LUPartialPivotDecompositionSuccessful(LU.matrixLU())!=MatrixSuccessCode::Success)
                        return SuccessCode::MatrixSolveFailure;
 
 
                    auto delta_z = LU.solve(-f);
                    next_space += delta_z;
 
                    if ( (delta_z.norm() < tracking_tolerance) && (ii >= (min_num_newton_iterations-1)) )
                        return SuccessCode::Success;
 
                    auto norm_J_inverse = LU.solve(RandomOfUnits<ComplexType>(S.NumVariables())).norm();
                    if (!amp::CriterionB(J.norm(), norm_J_inverse, max_num_newton_iterations - ii, tracking_tolerance, delta_z.norm(), AMP_config))
                        return SuccessCode::HigherPrecisionNecessary;
 
                    if (!amp::CriterionC(norm_J_inverse, next_space, tracking_tolerance, AMP_config))
                        return SuccessCode::HigherPrecisionNecessary;
                }
 
                return SuccessCode::FailedToConverge;
            }
            
 
 
 
 
 

            template <typename ComplexType, typename RealType>
            SuccessCode NewtonLoop(Vec<ComplexType> & next_space,
                                   RealType & norm_delta_z,
                                   RealType & norm_J,
                                   RealType & norm_J_inverse,
                                   RealType & condition_number_estimate,
                                   System const& S,
                                   Vec<ComplexType> const& current_space, // pass by value to get a copy of it
                                   ComplexType const& current_time, 
                                   RealType const& tracking_tolerance,
                                   unsigned min_num_newton_iterations,
                                   unsigned max_num_newton_iterations,
                                   AdaptiveMultiplePrecisionConfig const& AMP_config)
            {
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(max_num_newton_iterations >= min_num_newton_iterations && "max number newton iterations must be at least the min.");
                #endif
 
 
                next_space = current_space;
                for (unsigned ii = 0; ii < max_num_newton_iterations; ++ii)
                {
                    //TODO: wrap these into a single line.
                    auto f = S.Eval(next_space, current_time);
                    auto J = S.Jacobian(next_space, current_time);
                    auto LU = J.lu();
 
 
                    if (LUPartialPivotDecompositionSuccessful(LU.matrixLU())!=MatrixSuccessCode::Success)
                        return SuccessCode::MatrixSolveFailure;
 
 
                    auto delta_z = LU.solve(-f);
                    next_space += delta_z;
 
 
                    norm_delta_z = delta_z.norm();
                    norm_J = J.norm();
                    norm_J_inverse = LU.solve(RandomOfUnits<ComplexType>(S.NumVariables())).norm();
                    condition_number_estimate = norm_J*norm_J_inverse;
 
 
 
                    if ( (norm_delta_z < tracking_tolerance) && (ii >= (min_num_newton_iterations-1)) )
                        return SuccessCode::Success;
 
                    if (!amp::CriterionB(norm_J, norm_J_inverse, max_num_newton_iterations - ii, tracking_tolerance, norm_delta_z, AMP_config))
                        return SuccessCode::HigherPrecisionNecessary;
 
                    if (!amp::CriterionC(norm_J_inverse, next_space, tracking_tolerance, AMP_config))
                        return SuccessCode::HigherPrecisionNecessary;
                }
 
                return SuccessCode::FailedToConverge;
            }
        }
        
    }
    
}
 
#endif