/**
 * \copyright
 * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 */

#include "NumLib/NewtonRaphson.h"

#include <gtest/gtest.h>

#include <Eigen/LU>
#include <limits>
TEST(NumLibNewtonRaphson, Sqrt3)
{
    static const int N = 1;  // Problem's size.

    using LocalJacobianMatrix = Eigen::Matrix<double, N, N, Eigen::RowMajor>;
    using LocalResidualVector = Eigen::Matrix<double, N, 1>;

    Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(N);
    LocalJacobianMatrix jacobian;

    const int maximum_iterations = 5;
    const double tolerance = 1e-15;

    // Solve f(x) = x^2 - S == 0, where S = 3.
    //
    // Initial value
    double state = 1;

    auto const update_jacobian = [&state](LocalJacobianMatrix& jacobian)
    { jacobian(0, 0) = 2 * state; };

    auto const update_residual = [&state](LocalResidualVector& residual)
    { residual[0] = state * state - 3; };

    auto const update_solution = [&state](LocalResidualVector const& increment)
    { state += increment[0]; };

    auto const newton_solver =
        NumLib::NewtonRaphson<decltype(linear_solver), LocalJacobianMatrix,
                              decltype(update_jacobian), LocalResidualVector,
                              decltype(update_residual),
                              decltype(update_solution)>(
            linear_solver, update_jacobian, update_residual, update_solution,
            {maximum_iterations, tolerance, 0.0});
    auto const success_iterations = newton_solver.solve(jacobian);

    EXPECT_TRUE(static_cast<bool>(success_iterations));
    EXPECT_LE(*success_iterations, maximum_iterations);
    ASSERT_LE(state - std::sqrt(3), std::numeric_limits<double>::epsilon());
}