-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsample5.cpp
155 lines (136 loc) · 5.63 KB
/
sample5.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
/******************************************************
* C++ Library of the Linear Conjugate Gradient Methods (LibLCG)
*
* Copyright (C) 2022 Yi Zhang ([email protected])
*
* LibLCG 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 (LGPL) 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 LibLCG, please consider the option to obtain a commercial
* license for a fee. These licenses are offered by the LibLCG developing
* team. As a rule, licenses are provided "as-is", unlimited in time for
* a one time fee. Please send corresponding requests to: [email protected].
* 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.
******************************************************/
#include "../lib/lcg_eigen.h"
#include "iostream"
#include "Eigen/Dense"
#define M 1000
#define N 800
lcg_float max_diff(const Eigen::VectorXd &a, const Eigen::VectorXd &b)
{
lcg_float max = -1;
for (int i = 0; i < a.size(); i++)
{
max = lcg_max(sqrt((a[i] - b[i])*(a[i] - b[i])), max);
}
return max;
}
// 普通二维数组做核矩阵
Eigen::MatrixXd kernel = Eigen::MatrixXd::Random(M, N);
// 中间结果数组
Eigen::VectorXd tmp_arr(M);
Eigen::VectorXd p = Eigen::VectorXd::Constant(N, 1.0);
// 计算核矩阵乘向量的乘积
void CalAx(void* instance, const Eigen::VectorXd &x, Eigen::VectorXd &prod_Ax)
{
tmp_arr = kernel * x;
prod_Ax = kernel.transpose() * tmp_arr;
return;
}
void CalMx(void* instance, const Eigen::VectorXd &x, Eigen::VectorXd &prod_Mx)
{
prod_Mx = p.cwiseProduct(x);
return;
}
//定义共轭梯度监控函数
int Prog(void* instance, const Eigen::VectorXd *m, const lcg_float converge,
const lcg_para *param, const int k)
{
std::clog << "\rIteration-times: " << k << "\tconvergence: " << converge;
return 0;
}
int main(int argc, char const *argv[])
{
// 生成一组正演解
lcg_float LO = 1.0, HI = 2.0, Range = HI - LO;
Eigen::VectorXd fm = Eigen::VectorXd::Random(N);
fm = (fm + Eigen::VectorXd::Constant(N, 1.0))*0.5*Range;
fm = (fm + Eigen::VectorXd::Constant(N, LO));
// 计算共轭梯度B项
Eigen::VectorXd B(N);
tmp_arr = kernel * fm;
B = kernel.transpose() * tmp_arr;
/********************准备工作完成************************/
lcg_para self_para = lcg_default_parameters();
self_para.epsilon = 1e-5;
self_para.abs_diff = 0;
// 声明一组解
Eigen::VectorXd m = Eigen::VectorXd::Zero(N);
//Eigen::VectorXd p = Eigen::VectorXd::Constant(N, 1.0);
Eigen::VectorXd low = Eigen::VectorXd::Constant(N, LO);
Eigen::VectorXd hig = Eigen::VectorXd::Constant(N, HI);
std::clog << "solver: cg" << std::endl;
clock_t start = clock();
int ret = lcg_solver_eigen(CalAx, Prog, m, B, &self_para, NULL, LCG_CG);
clock_t end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
m.setZero();
std::clog << "solver: pcg" << std::endl;
start = clock();
ret = lcg_solver_preconditioned_eigen(CalAx, CalMx, Prog, m, B, &self_para, NULL, LCG_PCG);
end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
m.setZero();
std::clog << "solver: cgs" << std::endl;
start = clock();
ret = lcg_solver_eigen(CalAx, Prog, m, B, &self_para, NULL, LCG_CGS);
end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
m.setZero();
std::clog << "solver: bicgstab" << std::endl;
start = clock();
ret = lcg_solver_eigen(CalAx, Prog, m, B, &self_para, NULL, LCG_BICGSTAB);
end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
m.setZero();
std::clog << "solver: bicgstab2" << std::endl;
start = clock();
ret = lcg_solver_eigen(CalAx, Prog, m, B, &self_para, NULL, LCG_BICGSTAB2);
end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
m.setZero();
std::clog << "solver: pg" << std::endl;
start = clock();
ret = lcg_solver_constrained_eigen(CalAx, Prog, m, B, low, hig, &self_para, NULL, LCG_PG);
end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
m.setZero();
std::clog << "solver: spg" << std::endl;
start = clock();
ret = lcg_solver_constrained_eigen(CalAx, Prog, m, B, low, hig, &self_para, NULL, LCG_SPG);
end = clock();
std::clog << std::endl; lcg_error_str(ret);
std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
std::clog << "time use: "<<1000*(end-start)/(double)CLOCKS_PER_SEC<<" ms" << std::endl;
return 0;
}