/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* 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 .
*
* 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_CLCG_H
#define _GCTL_CLCG_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 clcg_solver() function.
*/
enum clcg_solver_type
{
/**
* Jacob's Bi-Conjugate Gradient Method
*/
CLCG_BICG,
/**
* Bi-Conjugate Gradient Method accelerated for complex symmetric A
*/
CLCG_BICG_SYM,
/**
* Conjugate Gradient Squared Method with real coefficients.
*/
CLCG_CGS,
/**
* Biconjugate gradient method.
*/
CLCG_BICGSTAB,
/**
* Transpose Free Quasi-Minimal Residual Method
*/
CLCG_TFQMR,
};
/**
* @brief return value of the clcg_solver() function
*/
enum clcg_return_code
{
CLCG_SUCCESS = 0, ///< The solver function terminated successfully.
CLCG_CONVERGENCE = 0, ///< The iteration reached convergence.
CLCG_STOP, ///< The iteration is stopped by the monitoring function.
CLCG_ALREADY_OPTIMIZIED, ///< The initial solution is already optimized.
// A negative number means a error
CLCG_UNKNOWN_ERROR = -1024, ///< Unknown error.
CLCG_INVILAD_VARIABLE_SIZE, ///< The variable size is negative
CLCG_INVILAD_MAX_ITERATIONS, ///< The maximal iteration times is negative.
CLCG_INVILAD_EPSILON, ///< The epsilon is negative.
CLCG_REACHED_MAX_ITERATIONS, ///< Iteration reached maximal limit.
CLCG_NAN_VALUE, ///< Nan value.
CLCG_INVALID_POINTER, ///< Invalid pointer.
CLCG_SIZE_NOT_MATCH, ///< Sizes of m and B do not match
CLCG_UNKNOWN_SOLVER, ///< Unknown solver
};
/**
* @brief Message type of the CLCG algorithms.
*
*/
enum clcg_message_type
{
CLCG_THROW, ///< throw error only
CLCG_ERROR, ///< display error only
CLCG_SOLUTION, ///< display info for evry solution
CLCG_ITERATION, ///< display info for every iteration
};
/**
* @brief Parameters of the conjugate gradient methods.
*/
struct clcg_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. For box-constrained methods,the convergence test
* is implemented using ||P(m-g) - m|| <= epsilon, in which P is the projector that
* transfers m into the constrained domain.
*/
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;
/**
* Minimal value for testing if a vector's module is bigger than the threshold.
* The default value is 1e-8
*
*/
double lamda;
};
class clcg_solver
{
private:
clcg_para clcg_param_;
size_t clcg_inter_;
clcg_message_type clcg_msg_;
// make them class variables are more suitable for repetitively usages
array > r1k, r2k, d1k, d2k;
array > uk, qk, wk;
array > Ax, Ad;
void clcg_error_str(clcg_return_code err_code, std::ostream &ss);
public:
clcg_solver();
virtual ~clcg_solver();
virtual void CLCG_Ax(const array > &x, array > &ax,
matrix_layout_e layout, conjugate_type_e conj) = 0;
virtual int CLCG_Progress(const array > &m, const double converge,
const clcg_para ¶m, size_t t, std::ostream &ss);
void set_clcg_message(clcg_message_type msg);
void set_clcg_report_interval(size_t inter);
void set_clcg_para(const clcg_para ¶m);
clcg_para default_clcg_para();
#ifdef GCTL_OPTIMIZATION_TOML
void set_clcg_para(std::string filename);
void set_clcg_para(const toml::value &toml_data);
#endif // GCTL_OPTIMIZATION_TOML
void CLCG_Minimize(array > &m, const array > &B,
clcg_solver_type solver_id = CLCG_CGS, std::ostream &ss = std::clog,
bool verbose = true, bool er_throw = false);
void clbicg(array > &m, const array > &B, std::ostream &ss = std::clog);
void clbicg_symmetric(array > &m, const array > &B, std::ostream &ss = std::clog);
void clcgs(array > &m, const array > &B, std::ostream &ss = std::clog);
void clbicgstab(array > &m, const array > &B, std::ostream &ss = std::clog);
void cltfqmr(array > &m, const array > &B, std::ostream &ss = std::clog);
};
}
#endif // _GCTL_CLCG_H