线性共轭梯度算法库 |
您所在的位置：网站首页 › maple定义参数a是b的共轭 › 线性共轭梯度算法库 |
|
线性共轭梯度算法库（C++ Library of Linear Conjugate Gradient，LIBLCG）
简介
安装
使用说明
自定义Ax计算函数
自定义进程监控函数
调用 lcg() 函数求解
调用 lpcg() 函数求解
示例
View on Github 线性共轭梯度算法库（C++ Library of Linear Conjugate Gradient，LIBLCG）浙江大学地球科学学院·地球物理研究所 简介liblcg是一个简单的C++线性共轭梯度算法库，包含了一般形式的共轭梯度与预优共轭梯度算法。可用于求解如下形式的线性方程组： Ax = B （可用共轭梯度求解）或者 PAx = PB（可用预优共轭梯度求解） 其中，A是一个N阶的对称矩阵、x为N*1的待求解的模型向量，B为N*1需拟合的目标向量。P则为预优算法中的预优矩阵，一般是一个N阶的对角阵。共轭梯度法广泛应用于无约束的线性最优化问题，拥有优良收敛与计算效率。 安装GCC（MacOS） mkdir build cd build cmake .. make make installWindows 请自行拷贝代码新建项目并编译…嘿嘿😁！ 使用说明 自定义Ax计算函数通常我们在使用共轭梯度法求解线性方程组Ax=B时A的维度可能会很大，直接储存A将消耗大量的内存空间，因此一般并不直接计算并储存A而是在需要的时候计算Ax的乘积。因此用户在使用liblcg时需要定义Ax的计算函数，同时提供初始解x与共轭梯度的B项（即拟合的对象）。如果使用预优方法还需要提供预优矩阵P项。Ax计算函数的形式必须满足算法库定义的一般形式： typedef void (*lcg_axfunc_ptr)(void* instance, const lcg_float* x, lcg_float* prod_Ax, const int n_size);函数接受4个参数，分别为： void *instance 传入的实例对象； const lcg_float *input_array Ax计算中的x数组的指针； lcg_float *output_array Ax的乘积； const int n_size 矩阵的大小。 自定义进程监控函数用户可以以下面的模版创建函数来显示共轭梯度迭代中的参数，并可以在适当的情况下停止迭代的进程。 typedef int (*lcg_progress_ptr)(void* instance, const lcg_float* m, const lcg_float converge, const lcg_para* param, const int n_size, const int k);函数接收6个参数，分别为： void* instance 传入的实例对象； const lcg_float* m 当前迭代的模型参数数组； const lcg_float converge 当前迭代的目标值； const lcg_para* param 当前迭代过程使用的参数； const int n_size 模型数组的大小； const int k 当前迭代的次数。 调用 lcg() 函数求解Ax计算函数定义完成后，用户需要给出方程组的初始解（一般可直接赋0）与共轭梯度的B项。然后调用 lcg() 函数求解： int lcg(lcg_axfunc_ptr Afp, lcg_progress_ptr Pfp, lcg_float* m, const lcg_float* B, const int n_size, const lcg_para* param, void* instance);函数接受7个参数，分别为： lcg_axfunc_ptr Afp 计算Ax的回调函数； lcg_progress_ptr Pfp 监控迭代过程的回调函数（非必须，无需监控时使用NULL参数即可）； lcg_float* m 模型参数数组，解得线性方程组的解也为这个数组； const lcg_float* B Ax = B 中的B项； const int n_size 模型参数数组的大小； const lcg_para* param 此次迭代使用的参数； void* instance 传入的实例对象。 调用 lpcg() 函数求解Ax计算函数定义完成后，用户需要给出方程组的初始解（一般可直接赋0）、共轭梯度的B项与预优矩阵P项（P一般为一个N阶的对角阵）。然后调用 lpcg() 函数求解： int lpcg(lcg_axfunc_ptr Afp, lcg_progress_ptr Pfp, lcg_float* m, const lcg_float* B, const lcg_float* P, const int n_size, const lcg_para* param, void* instance);函数接受8个参数，分别为： lcg_axfunc_ptr Afp 计算Ax的回调函数； lcg_progress_ptr Pfp 监控迭代过程的回调函数（非必须，无需监控时使用NULL参数即可）； lcg_float* m 模型参数数组，解得线性方程组的解也为这个数组； const lcg_float* B Ax = B 中的B项； const lcg_float* P 预优矩阵，一般是一个N阶的对角阵，这里直接用一个一维数组表示； const int n_size 模型参数数组的大小； const lcg_para* param 此次迭代使用的参数； void* instance 传入的实例对象。 示例 #include "../lib/lcg.h" #include "iostream" using std::clog; using std::endl; class TESTFUNC { public: TESTFUNC(); ~TESTFUNC(); void Routine(); /** * 因为类的成员函数指针不能直接被调用，所以我们在这里定义一个静态的中转函数来辅助Ax函数的调用 * 这里我们利用reinterpret_cast将_Ax的指针转换到Ax上，需要注意的是成员函数的指针只能通过 * 实例对象进行调用，因此需要void* instance变量。 */ static void _Ax(void* instance, const lcg_float* a, lcg_float* b, const int num) { return reinterpret_cast(instance)->Ax(a, b, num); } void Ax(const lcg_float* a, lcg_float* b, const int num); //定义共轭梯度中Ax的算法 static int _Progress(void* instance, const lcg_float* m, const lcg_float converge, const lcg_para *param, const int n_size, const int k) { return reinterpret_cast(instance)->Progress(m, converge, param, n_size, k); } int Progress(const lcg_float* m, const lcg_float converge, const lcg_para *param, const int n_size, const int k); private: lcg_float* m_; lcg_float* b_; lcg_float* p_; lcg_float kernel_[3][3]; }; TESTFUNC::TESTFUNC() { // 测试线性方程组 // 6.3*x1 + 3.9*x2 + 2.5*x3 = -2.37 // 3.9*x1 + 1.2*x2 + 3.1*x3 = 5.82 // 2.5*x1 + 3.1*x2 + 7.6*x3 = 5.21 // 目标解 x1=1.2 x2=-3.7 x3=1.8 // 注意根据共轭梯度法的要求 kernel是一个N阶对称阵 kernel_[0][0] = 6.3; kernel_[0][1] = 3.9; kernel_[0][2] = 2.5; kernel_[1][0] = 3.9; kernel_[1][1] = 1.2; kernel_[1][2] = 3.1; kernel_[2][0] = 2.5; kernel_[2][1] = 3.1; kernel_[2][2] = 7.6; // 初始解 m_ = lcg_malloc(3); m_[0] = 0.0; m_[1] = 0.0; m_[2] = 0.0; // 拟合目标值（含有一定的噪声） b_ = lcg_malloc(3); b_[0] = -2.3723; b_[1] = 5.8221; b_[2] = 5.2165; // 测试预优矩阵 这里只是测试流程 预优矩阵值全为1 并没有什么作用 p_ = lcg_malloc(3); p_[0] = p_[1] = p_[2] = 1.0; } TESTFUNC::~TESTFUNC() { lcg_free(m_); lcg_free(b_); lcg_free(p_); } void TESTFUNC::Ax(const lcg_float* a, lcg_float* b, const int num) { for (int i = 0; i |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |