388 lines
14 KiB
C++
388 lines
14 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#ifndef _GCTL_LCG_H
|
|
#define _GCTL_LCG_H
|
|
|
|
#include "gctl/core.h"
|
|
#include "gctl/maths.h"
|
|
#include "gctl/algorithm.h"
|
|
|
|
#include "gctl_optimization_config.h"
|
|
#ifdef GCTL_OPTIMIZATION_TOML
|
|
#include "toml.hpp"
|
|
#endif // GCTL_OPTIMIZATION_TOML
|
|
|
|
#if defined _WINDOWS || __WIN32__
|
|
#include "windows.h"
|
|
#endif // _WINDOWS || __WIN32__
|
|
|
|
namespace gctl
|
|
{
|
|
/**
|
|
* @brief Types of method that could be recognized by the lcg_solver() function.
|
|
*/
|
|
enum lcg_solver_type
|
|
{
|
|
/**
|
|
* Conjugate gradient method.
|
|
*/
|
|
LCG_CG,
|
|
/**
|
|
* Preconditioned conjugate gradient method.
|
|
*/
|
|
LCG_PCG,
|
|
/**
|
|
* Conjugate gradient squared method.
|
|
*/
|
|
LCG_CGS,
|
|
/**
|
|
* Biconjugate gradient method.
|
|
*/
|
|
LCG_BICGSTAB,
|
|
/**
|
|
* Biconjugate gradient method with restart.
|
|
*/
|
|
LCG_BICGSTAB2,
|
|
/**
|
|
* Conjugate gradient method with projected gradient for inequality constraints.
|
|
* This algorithm comes without non-monotonic linear search for the step length.
|
|
*/
|
|
LCG_PG,
|
|
/**
|
|
* Conjugate gradient method with spectral projected gradient for inequality constraints.
|
|
* This algorithm comes with non-monotonic linear search for the step length.
|
|
*/
|
|
LCG_SPG,
|
|
};
|
|
|
|
/**
|
|
* @brief return value of the lcg_solver() function
|
|
*/
|
|
enum lcg_return_code
|
|
{
|
|
LCG_SUCCESS = 0, ///< The solver function terminated successfully.
|
|
LCG_CONVERGENCE = 0, ///< The iteration reached convergence.
|
|
LCG_STOP, ///< The iteration is stopped by the monitoring function.
|
|
LCG_ALREADY_OPTIMIZIED, ///< The initial solution is already optimized.
|
|
// A negative number means a error
|
|
LCG_UNKNOWN_ERROR = -1024, ///< Unknown error.
|
|
LCG_INVILAD_VARIABLE_SIZE, ///< The variable size is negative
|
|
LCG_INVILAD_MAX_ITERATIONS, ///< The maximal iteration times is negative.
|
|
LCG_INVILAD_EPSILON, ///< The epsilon is negative.
|
|
LCG_INVILAD_RESTART_EPSILON, ///< The restart epsilon is negative.
|
|
LCG_REACHED_MAX_ITERATIONS, ///< Iteration reached maximal limit.
|
|
LCG_NULL_PRECONDITION_MATRIX, ///< Null precondition matrix.
|
|
LCG_NAN_VALUE, ///< Nan value.
|
|
LCG_INVALID_POINTER, ///< Invalid pointer.
|
|
LCG_INVALID_LAMBDA, ///< Invalid range for lambda.
|
|
LCG_INVALID_SIGMA, ///< Invalid range for sigma.
|
|
LCG_INVALID_BETA, ///< Invalid range for beta.
|
|
LCG_INVALID_MAXIM, ///< Invalid range for maxi_m.
|
|
LCG_SIZE_NOT_MATCH, ///< Sizes of m and B do not match
|
|
};
|
|
|
|
/**
|
|
* @brief Message type of the LCG algorithms.
|
|
*
|
|
*/
|
|
enum lcg_message_type
|
|
{
|
|
LCG_THROW, ///< throw error only
|
|
LCG_ERROR, ///< display error only
|
|
LCG_SOLUTION, ///< display info for evry solution
|
|
LCG_ITERATION, ///< display info for every iteration
|
|
};
|
|
|
|
/**
|
|
* @brief Parameters of the conjugate gradient methods.
|
|
*/
|
|
struct lcg_para
|
|
{
|
|
/**
|
|
* Maximal iteration times. The process will continue till the convergence is met
|
|
* if this option is set to zero (default).
|
|
*/
|
|
int max_iterations;
|
|
|
|
/**
|
|
* Epsilon for convergence test.
|
|
* This parameter determines the accuracy with which the solution is to be
|
|
* found. A minimization terminates when ||g||/max(||g0||, 1.0) <= epsilon or
|
|
* sqrt(||g||)/N <= epsilon for the lcg_solver() function, where ||.|| denotes
|
|
* the Euclidean (L2) norm. The default value of epsilon is 1e-8.
|
|
*/
|
|
double epsilon;
|
|
|
|
/**
|
|
* Whether to use absolute mean differences (AMD) between |Ax - B| to evaluate the process.
|
|
* The default value is false which means the gradient based evaluating method is used.
|
|
* The AMD based method will be used if this variable is set to true. This parameter is only
|
|
* applied to the non-constrained methods.
|
|
*/
|
|
int abs_diff;
|
|
|
|
/**
|
|
* Restart epsilon for the LCG_BICGSTAB2 algorithm. The default value is 1e-6
|
|
*/
|
|
double restart_epsilon;
|
|
|
|
/**
|
|
* Initial step length for the project gradient method. The default is 1.0
|
|
*/
|
|
double step;
|
|
|
|
/**
|
|
* multiplier for updating solutions with the spectral projected gradient method. The range of
|
|
* this variable is (0, 1). The default is given as 0.95
|
|
*/
|
|
double sigma;
|
|
|
|
/**
|
|
* descending ratio for conducting the non-monotonic linear search. The range of
|
|
* this variable is (0, 1). The default is given as 0.9
|
|
*/
|
|
double beta;
|
|
|
|
/**
|
|
* The maximal record times of the objective values for the SPG method. The method use the
|
|
* objective values from the most recent maxi_m times to preform the non-monotonic linear search.
|
|
* The default value is 10.
|
|
*/
|
|
int maxi_m;
|
|
};
|
|
|
|
/**
|
|
* @brief This abstract class implements conjugate gradient algorithms for solving
|
|
* a linear system like Ax = B where A is a N*N matrix. As the actual variable
|
|
* we need is the product of 'Ax', the kernel matrix 'A' is not declared within
|
|
* the class definition. Instead, a pure virtual function is declared as the
|
|
* callback interface for calculating the product of 'Ax' as 'void LCG_Ax(const
|
|
* array<double> &x, array<double> &ax)'. A virtual function "int LCG_Progress(const
|
|
* array<double> &m, const double converge, const lcg_para ¶m, size_t t,
|
|
* std::ostream &ss)" could be reloaded for customed convergence tests.
|
|
*/
|
|
class lcg_solver
|
|
{
|
|
private:
|
|
lcg_para lcg_param_;
|
|
size_t lcg_inter_;
|
|
lcg_message_type lcg_msg_;
|
|
|
|
// make them class variables are more suitable for repetitively usages
|
|
array<double> zk, gk, dk, Adk;
|
|
array<double> rk, r0_T, pk, vk;
|
|
array<double> Apx, uk, qk, qk_m, wk;
|
|
array<double> m_new, gk_new;
|
|
array<double> sk, yk;
|
|
|
|
/**
|
|
* @brief Display info of a given return code. This is a private function
|
|
* and can only be called by other class functions.
|
|
*
|
|
* @param err_code Input retrun code
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lcg_error_str(lcg_return_code err_code, std::ostream &ss);
|
|
|
|
public:
|
|
lcg_solver(); ///< default constructor
|
|
virtual ~lcg_solver(); ///< default de-constructor
|
|
|
|
/**
|
|
* @brief Callback interface for calculating the product of 'A' multipled by an arbitrary vector 'x'.
|
|
*
|
|
* @param x Multipler
|
|
* @param ax Product of Ax
|
|
*/
|
|
virtual void LCG_Ax(const array<double> &x, array<double> &ax) = 0;
|
|
|
|
/**
|
|
* @brief Callback interface for calculating the product of 'M' multipled by an arbitrary vector 'x'.
|
|
* In which 'M' is the inverse of the pre-conditioning matrix. This function must be reloaded for the
|
|
* LCG_PCG algorithm.
|
|
*
|
|
* @param x Multipler
|
|
* @param mx Product of Ax
|
|
*/
|
|
virtual void LCG_Mx(const array<double> &x, array<double> &mx);
|
|
|
|
/**
|
|
* @brief Utility function for monitoring the solving process.
|
|
*
|
|
* @param m Current solution
|
|
* @param converge Current convergence
|
|
* @param param Employed parameters
|
|
* @param t Current iterative times
|
|
* @param ss Output stream of runtime info.
|
|
* @return Quit the solver if returned a non-zero value.
|
|
*/
|
|
virtual int LCG_Progress(const array<double> &m, const double converge,
|
|
const lcg_para ¶m, size_t t, std::ostream &ss);
|
|
|
|
/**
|
|
* @brief Set the lcg message object
|
|
*
|
|
* @param msg Input message type.
|
|
*/
|
|
void set_lcg_message(lcg_message_type msg);
|
|
|
|
/**
|
|
* @brief Set the lcg report intervals
|
|
*
|
|
* @param inter Input reprot intervals.
|
|
*/
|
|
void set_lcg_report_interval(size_t inter);
|
|
|
|
/**
|
|
* @brief Set the lcg para object
|
|
*
|
|
* @param param Input lcg parameters.
|
|
*/
|
|
void set_lcg_para(const lcg_para ¶m);
|
|
|
|
/**
|
|
* @brief Return a lcg_para object with default values
|
|
*
|
|
* @return lcg_para
|
|
*/
|
|
lcg_para default_lcg_para();
|
|
|
|
#ifdef GCTL_OPTIMIZATION_TOML
|
|
|
|
/**
|
|
* @brief Set parameters of the conjugate gradient algorithms using a toml file.
|
|
* All parameter options must be listed under a top-level table 'lcg'. Available options
|
|
* under the 'lcg' table are as declared in the lcg_para structure.
|
|
*
|
|
* @param filename Input configuration file
|
|
*/
|
|
void set_lcg_para(std::string filename);
|
|
|
|
/**
|
|
* @brief Set parameters of the conjugate gradient algorithms using a toml::value object.
|
|
* All parameter options must be listed under a top-level table 'lcg'. Available options
|
|
* under the 'lcg' table are as declared in the lcg_para structure.
|
|
*
|
|
* @param toml_data Input toml data
|
|
*/
|
|
void set_lcg_para(const toml::value &toml_data);
|
|
|
|
#endif // GCTL_OPTIMIZATION_TOML
|
|
|
|
/**
|
|
* @brief The allback interface for all CG algorithms. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param solver_id Selected solver type.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void LCG_Minimize(array<double> &m, const array<double> &b, lcg_solver_type solver_id = LCG_CG, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The allback interface for all CG algorithms. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param solver_id Selected solver type.
|
|
* @param low Lower bound of acceptable solutions.
|
|
* @param hig Higher bound of acceptable solutions.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void LCG_MinimizeConstrained(array<double> &m, const array<double> &b, const array<double> &low, const array<double> &hig, lcg_solver_type solver_id = LCG_PG, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the CG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lcg(array<double> &m, const array<double> &B, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the PCG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lpcg(array<double> &m, const array<double> &B, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the CGS algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lcgs(array<double> &m, const array<double> &B, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the BICGSTAB algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lbicgstab(array<double> &m, const array<double> &B, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the BICGSTAB2 algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lbicgstab2(array<double> &m, const array<double> &B, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the PG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param low Lower bound of acceptable solutions.
|
|
* @param hig Higher bound of acceptable solutions.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lpg(array<double> &m, const array<double> &B, const array<double> &low, const array<double> &hig, std::ostream &ss = std::clog);
|
|
|
|
/**
|
|
* @brief The standalone callback interface for the SPG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
|
|
*
|
|
* @param m Initial/Input solution.
|
|
* @param B Right hand term of the system system.
|
|
* @param low Lower bound of acceptable solutions.
|
|
* @param hig Higher bound of acceptable solutions.
|
|
* @param ss Output stream of runtime info.
|
|
*/
|
|
void lspg(array<double> &m, const array<double> &B, const array<double> &low, const array<double> &hig, std::ostream &ss = std::clog);
|
|
};
|
|
}
|
|
|
|
#endif // _GCTL_LCG_H
|