diff --git a/src/dfpmsolver/driver_structs.hpp b/src/dfpmsolver/driver_structs.hpp
index 2464af7..104e2a9 100644
--- a/src/dfpmsolver/driver_structs.hpp
+++ b/src/dfpmsolver/driver_structs.hpp
@@ -1,52 +1,52 @@
#ifndef SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
#define SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
#include "common.hpp"
-#include "utils/sparse_solver_structs.hpp"
+#include "utils/sparse_solvers/sparse_solver_structs.hpp"
namespace specmicp {
namespace dfpmsolver {
//! \brief Options of a driver
//!
struct DriverOptions
{
scalar_t residuals_tolerance; //!< Tolerance for the residual
scalar_t step_tolerance; //!< Tolerance for the minimum step length
scalar_t threshold_stationary_point; //!< if ||R||>threshold, the point is classified as stationary
int maximum_iterations; //!< Maximum iterations allowed
scalar_t maximum_step_length; //!< Maximum step length allowed
int max_iterations_at_max_length; //!< Maximum number of iterations at maximum step length
scalar_t coeff_accept_newton_step; //!< Accept Newton step if enough progress is made
- SparseSolver sparse_solver; //!< The sparse solver to use
+ sparse_solvers::SparseSolver sparse_solver; //!< The sparse solver to use
DriverOptions():
residuals_tolerance(5e-5),
step_tolerance(1e-10),
threshold_stationary_point(1e-4),
maximum_iterations(200),
maximum_step_length(1e3),
max_iterations_at_max_length(50),
coeff_accept_newton_step(0.9),
- sparse_solver(SparseSolver::SparseQR)
+ sparse_solver(sparse_solvers::SparseSolver::SparseQR)
{}
};
//! \brief Performance of a driver
struct DriverPerformance
{
int nb_call_residuals;
int nb_call_jacobian;
int nb_iterations;
int nb_consecutive_max_step_taken;
int nb_max_step_taken;
bool maximum_step_taken;
scalar_t current_residual;
scalar_t current_update;
};
} // end namespace dfpmsolver
} // end namespace specmicp
#endif // SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
diff --git a/src/dfpmsolver/parabolic_driver.hpp b/src/dfpmsolver/parabolic_driver.hpp
index 2b5fd80..81a9d2f 100644
--- a/src/dfpmsolver/parabolic_driver.hpp
+++ b/src/dfpmsolver/parabolic_driver.hpp
@@ -1,147 +1,146 @@
#ifndef SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
#define SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
#include "common.hpp"
-#include <Eigen/Sparse>
+#include "utils/sparse_solvers/sparse_solver.hpp"
#include "driver.hpp"
#include "parabolic_structs.hpp"
namespace specmicp {
namespace dfpmsolver {
//! \brief The parabolic driver
//!
//! Parabolic driver for finite element programs
template <class Program>
class ParabolicDriver: public Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>
{
public:
using base=Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>;
using base::program;
using base::get_neq;
using base::get_options;
using base::get_perfs;
using base::scaling;
using base::initialize_scaling;
ParabolicDriver(Program& the_program):
Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>(the_program),
m_velocity(Vector::Zero(the_program.get_tot_ndf()))
{}
//! \brief Solve a timestep of length dt
ParabolicDriverReturnCode solve_timestep(scalar_t dt, Vector& displacement);
//! \brief Restart the current timestep
ParabolicDriverReturnCode restart_timestep(Vector& displacement);
//! \brief Check if the solution has been found
ParabolicDriverReturnCode check_convergence();
//! \brief Compute the residuals
void compute_residuals(const Vector& displacement,
Vector& residual)
{
compute_residuals(displacement, m_velocity, residual);
}
//! \brief Compute the residuals
void compute_residuals(const Vector& displacement,
const Vector& velocity,
Vector& residual)
{
program().compute_residuals(displacement, velocity, residual);
get_perfs().nb_call_residuals += 1;
}
//! \brief Compute the jacobian
void compute_jacobian(Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian
);
//! \brief Return the norm of the current residuals
double norm_residuals() {return m_residuals.norm();}
//! \brief Read/Write reference to the velocity vector
const Vector& get_velocity() const {return m_velocity;}
Vector& velocity() {return m_velocity;}
void set_velocity(Vector& velocity_vector);
//! \brief Initialize the computation
void initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement);
//! \brief Strang Linesearch
//!
//! ref :
//! - Matthies et al. (1979)
//! - JHP course notes
ParabolicLinesearchReturnCode strang_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda
);
scalar_t residuals_0() {return m_norm_0;}
private:
void reset_velocity();
//! \brief Backtracking Linesearch
//!
//! ref :
//! - Algo A6.3.1 : Dennis and Schnabel (1983)
//! - Nocedal & Wrigth (2006)
ParabolicLinesearchReturnCode backtracking_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda_out
);
//! Update the variables in displacement
void update_variable(
const Vector& update,
scalar_t lambda,
Vector& displacement
);
//! \brief Set the predictor
void set_predictor(Vector& displacement);
- //! \brief solve the linear system
- SparseSolverReturnCode solve_linear_system(Eigen::VectorXd& update);
//! \brief perform the linesearch
ParabolicDriverReturnCode linesearch(
Vector &update,
Vector &displacements
);
//! Compute the residuals for the linesearch
double compute_residuals_linesearch(
Vector& update,
scalar_t lambda,
Vector& displacement
);
double compute_residuals_strang_linesearch(
Vector &update,
scalar_t lambda,
Vector &displacement
);
scalar_t update_norm(const Vector& update);
double m_norm_0;
double m_current_dt;
Eigen::VectorXd m_gradient;
Eigen::VectorXd m_velocity;
Eigen::VectorXd m_residuals;
Eigen::VectorXd m_predictor;
Eigen::SparseMatrix<double> m_jacobian;
+ sparse_solvers::SparseSolverPtr<Eigen::SparseMatrix<double>, Vector, Vector> m_solver;
bool m_velocity_is_initialized;
};
} // end namespace dfpmsolver
} // end namespace specmicp
// implementation
// ==============
#include "parabolic_driver.inl"
#endif // SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
diff --git a/src/dfpmsolver/parabolic_driver.inl b/src/dfpmsolver/parabolic_driver.inl
index 6f629c8..0079a30 100644
--- a/src/dfpmsolver/parabolic_driver.inl
+++ b/src/dfpmsolver/parabolic_driver.inl
@@ -1,471 +1,466 @@
#include "parabolic_driver.hpp" // for syntaxic coloration only
#include "utils/log.hpp"
-#include "utils/sparse_solver.hpp"
+//#include "utils/sparse_solver.hpp"
namespace specmicp {
namespace dfpmsolver {
template <class Program>
void ParabolicDriver<Program>::compute_jacobian(Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian
)
{
program().compute_jacobian(displacement, velocity, jacobian, get_options().alpha*m_current_dt);
jacobian = jacobian*scaling().asDiagonal();
get_perfs().nb_call_jacobian += 1;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::solve_timestep(scalar_t dt, Eigen::VectorXd& displacement)
{
initialize_timestep(dt, displacement);
return restart_timestep(displacement);
}
template <class Program>
void ParabolicDriver<Program>::initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement)
{
initialize_scaling();
m_residuals = Eigen::VectorXd::Zero(get_neq());
m_current_dt = dt;
set_predictor(displacement);
reset_velocity();
program().apply_bc(dt, displacement, m_velocity);
compute_residuals(displacement, m_velocity, m_residuals);
m_norm_0 = norm_residuals();
get_perfs().nb_iterations = 0;
}
template <class Program>
void ParabolicDriver<Program>::set_predictor(Vector& displacement)
{
if (get_options().alpha < 1)
m_predictor = displacement + (1-get_options().alpha)*m_current_dt*m_velocity;
else
m_predictor = displacement;
}
template <class Program>
scalar_t ParabolicDriver<Program>::update_norm(const Vector& update)
{
// l-∞ scaled norm
scalar_t norm = 0.0;
for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
{
const index_t id_eq = program().id_equation(dof);
if (id_eq == no_species or update(id_eq) == 0.0) continue;
norm = std::max(norm, std::abs(update(id_eq))/(std::max(std::abs(m_velocity(dof)), 1.0/scaling(id_eq))));
}
return norm;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::restart_timestep(Vector& displacement)
{
ParabolicDriverReturnCode return_code = ParabolicDriverReturnCode::NotConvergedYet;
+ m_solver = sparse_solvers::get_sparse_solver<
+ Eigen::SparseMatrix<double>, Vector, Vector>(get_options().sparse_solver);
Eigen::VectorXd update(get_neq());
update.setZero();
get_perfs().current_update = 0;
+ bool force_recompute_jacobian = true;
while (return_code == ParabolicDriverReturnCode::NotConvergedYet)
{
compute_residuals(displacement, m_velocity, m_residuals);
get_perfs().current_residual = m_residuals.norm()/m_norm_0;
get_perfs().current_update = update_norm(update);
DEBUG << " NB iterations : " << get_perfs().nb_iterations
<< " - res : " << get_perfs().current_residual
<< " - update : " << get_perfs().current_update;
return_code = check_convergence();
if (return_code != ParabolicDriverReturnCode::NotConvergedYet) break;
- if (get_perfs().nb_iterations % get_options().quasi_newton == 0)
+ if (force_recompute_jacobian
+ or get_perfs().nb_iterations % get_options().quasi_newton == 0)
{
- // Todo : avoid computation of decomposition if it is the same matrix
compute_jacobian(displacement, m_velocity, m_jacobian);
+ m_solver->decompose(m_jacobian);
}
get_perfs().nb_iterations += 1;
m_gradient = m_jacobian.transpose()*m_residuals;
- SparseSolverReturnCode retcode = solve_linear_system(update);
- if (retcode != SparseSolverReturnCode::Success)
+ sparse_solvers::SparseSolverReturnCode retcode = m_solver->solve_scaling(m_residuals, scaling(), update);
+ if (retcode != sparse_solvers::SparseSolverReturnCode::Success)
{
ERROR << "Error when solving linear system : " << (int) retcode << std::endl;
return ParabolicDriverReturnCode::ErrorLinearSystem;
}
//if (update.norm() < get_options().step_tolerance) return_code = ParabolicDriverReturnCode::ErrorMinimized;
//else
return_code = linesearch(update, displacement);
}
return return_code;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::check_convergence()
{
ParabolicDriverReturnCode termcode = ParabolicDriverReturnCode::NotConvergedYet;
const int nb_iterations = get_perfs().nb_iterations;
const scalar_t norm_residuals = get_perfs().current_residual;
const scalar_t norm_update = get_perfs().current_update;
//std::cout << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0 << std::endl;
DEBUG << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0;
if (norm_residuals < get_options().residuals_tolerance)
{
termcode = ParabolicDriverReturnCode::ResidualMinimized;
}
else if (nb_iterations > 0 and norm_update > 0 and norm_update < 1.01*get_options().step_tolerance)
{
if (norm_residuals > get_options().threshold_stationary_point)
{
ERROR << "Stationary point detected !";
termcode = ParabolicDriverReturnCode::StationaryPoint;
}
WARNING << "Error is minimized - may indicate a stationnary point";
termcode = ParabolicDriverReturnCode::ErrorMinimized;
}
else if (nb_iterations > get_options().maximum_iterations)
{
ERROR << "Maximum number of iteration reached (" << get_options().maximum_iterations << ")";
termcode = ParabolicDriverReturnCode::MaxIterations;
}
else if (get_perfs().maximum_step_taken)
{
get_perfs().nb_consecutive_max_step_taken += 1;
get_perfs().nb_max_step_taken += 1;
if (get_perfs().nb_consecutive_max_step_taken == get_options().max_iterations_at_max_length) {
ERROR << "Divergence detected - Maximum step length taken two many times";
termcode = ParabolicDriverReturnCode::MaxStepTakenTooManyTimes;
}
}
else
{
get_perfs().nb_consecutive_max_step_taken = 0;
}
get_perfs().return_code = termcode;
return termcode;
}
-template <class Program>
-SparseSolverReturnCode ParabolicDriver<Program>::solve_linear_system(Eigen::VectorXd& update)
-{
- return solve_sparse_linear_system(scaling().asDiagonal()*m_residuals,
- m_jacobian,
- update,
- get_options().sparse_solver);
-}
-
template <class Program>
double ParabolicDriver<Program>::compute_residuals_linesearch(
Vector &update,
scalar_t lambda,
Vector &displacement
)
{
Eigen::VectorXd velocity(m_velocity);
Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq());
program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity);
compute_residuals(displacement, velocity, residual);
return 0.5*residual.squaredNorm();
}
template <class Program>
double ParabolicDriver<Program>::compute_residuals_strang_linesearch(
Vector &update,
scalar_t lambda,
Vector &displacement
)
{
Eigen::VectorXd velocity(m_velocity);
Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq());
program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity);
compute_residuals(displacement, velocity, residual);
return update.dot(residual);
}
template <class Program>
void ParabolicDriver<Program>::update_variable(
const Vector& update,
scalar_t lambda,
Vector& displacement
)
{
program().update_solution(update, lambda, get_options().alpha*m_current_dt,
m_predictor, displacement, m_velocity);
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::linesearch(
Vector &update,
Vector &displacements
)
{
base::is_step_too_long(update);
get_perfs().maximum_step_taken = false;
scalar_t lambda;
ParabolicLinesearchReturnCode retcode;
switch (get_options().linesearch) {
case ParabolicLinesearch::Bactracking:
retcode = backtracking_linesearch(update, displacements, lambda);
break;
case ParabolicLinesearch::Strang:
retcode = strang_linesearch(update, displacements, lambda);
break;
default:
throw std::runtime_error("Linesearch type for Parabolic driver is not recognized");
break;
}
if (retcode != ParabolicLinesearchReturnCode::Success)
{
return ParabolicDriverReturnCode::LinesearchFailed;
}
update_variable(update, lambda, displacements);
update *= lambda;
return ParabolicDriverReturnCode::NotConvergedYet;
}
namespace internal
{
//! \brief Return true if a and b have the same sign
inline bool have_same_sign(double a, double b)
{
return (std::copysign(1., a) * std::copysign(1., b) > 0);
}
} // end namespace internal
template <class Program>
ParabolicLinesearchReturnCode ParabolicDriver<Program>::strang_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda
)
{
DEBUG << "Strang linesearch";
Eigen::VectorXd xp(displacements);
const scalar_t s_tol = 0.5;
const scalar_t s_max = 2.0;
const int lin_max = 10;
scalar_t s_b = 0.0;
scalar_t g_b = 0.5*m_residuals.squaredNorm();
scalar_t s_a = 1.0;
scalar_t g_a = compute_residuals_strang_linesearch(update, 1.0, xp);
scalar_t newtlen = (scaling().asDiagonal()*update).norm();
// const scalar_t s_r = s_a;
const scalar_t g_r = g_a;
if (std::abs(g_a) <= s_tol*std::abs(g_b))
{
DEBUG << "Skip linesearch ";
lambda = 1.0;
if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
return ParabolicLinesearchReturnCode::Success;
}
while (internal::have_same_sign(g_a, g_b) and s_a < s_max)
{
s_b = s_a;
s_a = 2*s_a;
g_b = g_a;
g_a = compute_residuals_strang_linesearch(update, s_a, xp);
}
scalar_t g_ = g_a;
scalar_t g_0 = g_a;
scalar_t s = s_a;
int l;
for (l=0; l<lin_max; ++l)
{
if (internal::have_same_sign(g_a, g_b) or
std::abs(g_) < s_tol*std::abs(g_0)
or std::abs(s_a - s_b) < s_tol*(s_a + s_b)/2 )
break;
s = s_a - g_a * ( s_a - s_b)/(g_a - g_b);
g_ = compute_residuals_strang_linesearch(update, s, xp);
if (internal::have_same_sign(g_, g_a))
g_b = g_b/2.;
else
{
s_b = s_a;
g_b = g_a;
}
s_a = s;
g_a = g_;
}
if (l == lin_max) return ParabolicLinesearchReturnCode::MaximumIterations;
if (g_a >= g_r) {
WARNING << "Failed to find better update in Strang linesearch";
//lambda = 0.1*s_r;
return backtracking_linesearch(update, displacements, lambda);
}
lambda = s_a;
if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
return ParabolicLinesearchReturnCode::Success;
}
template <class Program>
ParabolicLinesearchReturnCode ParabolicDriver<Program>::backtracking_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda_out
)
{
// References
// ----------
// - Algo A6.3.1 : Dennis and Schnabel (1983)
// - Nocedal & Wrigth (2006)
DEBUG << "Linesearch";
Eigen::VectorXd xp(displacements);
double fcp;
int retcode = 2; // 2 not converged, 1 problem, 0 success
const scalar_t alpha = 1e-4;
scalar_t newtlen = (scaling().asDiagonal()*update).norm();
scalar_t init_slope = m_gradient.dot(update);
const scalar_t rellength = update_norm(update);
const scalar_t minlambda = get_options().step_tolerance / rellength;
scalar_t lambda = 1.0;
scalar_t lambda_prev = lambda;
scalar_t merit_value = 0.5*m_residuals.squaredNorm();
// new residual
fcp = compute_residuals_linesearch(update, lambda, xp);
// Skip linesearch if enough progress is done
// ------------------------------------------
if (fcp < get_options().coeff_accept_newton_step *merit_value)
{
DEBUG << "Skip linesearch ";
lambda_out = 1.0;
return ParabolicLinesearchReturnCode::Success;
}
// The linesearch
// --------------
scalar_t fc = merit_value;
scalar_t fcp_prev;
int cnt = 0;
do
{
SPAM << "cnt : " <<cnt << " - lambda : " << lambda << " # " << fc << " # " << fcp << " # " << init_slope;
if (fcp <= fc + alpha*lambda*init_slope) //pg760 Fachinei2003
{
retcode = 0;
if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
break;
}
else if (lambda < minlambda)
{
retcode = 0;
lambda = minlambda;
//retcode = 1;
break;
}
else
{
//WARNING << "denom : " << fcp - fc -init_slope;
// Select a new step length
// - - - - - - - - - - - -
double lambdatmp;
if (cnt == 0) { // only a quadratic at the first
lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
}
else
{
const scalar_t factor = 1.0 /(lambda - lambda_prev);
const scalar_t x1 = fcp - fc - lambda*init_slope;
const scalar_t x2 = fcp_prev - fc - lambda_prev*init_slope;
const scalar_t a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
const scalar_t b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
if (a == 0)
{ // cubic interpolation is in fact a quadratic interpolation
lambdatmp = - init_slope/(2*b);
}
else
{
const scalar_t disc = b*b-3*a*init_slope;
lambdatmp = (-b+std::sqrt(disc))/(3*a);
}
if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
}
//WARNING << "lambdatmp : " << lambdatmp;
lambda_prev = lambda;
fcp_prev = fcp;
if (lambdatmp < 0.1*lambda) {
lambda = 0.1 * lambda;
} else {
lambda = lambdatmp;
}
}
if (not std::isfinite(lambda))
{
ERROR << "Lambda is non finite - we stop";
return ParabolicLinesearchReturnCode::Divergence;
}
fcp = compute_residuals_linesearch(update, lambda, xp);
//xp.velocity = x.velocity;
++cnt;
} while(retcode == 2 and cnt < 50);
//WARNING << "Lambda : " << lambda << " - iterations : " << cnt;
if (cnt == 50)
{
ERROR << "Too much linesearch iterations ! We stop";
return ParabolicLinesearchReturnCode::MaximumIterations;
}
lambda_out = lambda;
switch (retcode) {
case 0:
return ParabolicLinesearchReturnCode::Success;
break;
case 1:
return ParabolicLinesearchReturnCode::LambdaTooSmall;
default:
return ParabolicLinesearchReturnCode::NotSupposedToHappen;
break;
}
}
template <class Program>
void ParabolicDriver<Program>::set_velocity(Vector& velocity_vector)
{
m_velocity.resize(program().get_tot_ndf());
m_velocity.setZero();
for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
{
if (program().id_equation(dof) == no_equation)
{
m_velocity(dof) = velocity_vector(dof);
}
}
}
template <class Program>
void ParabolicDriver<Program>::reset_velocity()
{
for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
{
if (program().id_equation(dof) == no_equation) continue;
m_velocity(dof) = 0.0;
}
}
} // end namespace dfpmsolver
} // end namespace specmicp