# Langevin Preconditioning Matrix

Feb 26th, 2024 (edited)
811
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. MatrixXd precond_3dim_new(const VectorXd &p, const VectorXd &params) {
2.     int n = p.size();
3.
4.     double r1 = params(0);
5.     double r2 = params(1);
6.     double r3 = params(2);
7.     double k  = params(3);
8.
9.     double x = p(0), y = p(1), z = p(2);
10.     double phi = atan2(y, x);
11.     double theta = phi * k;
12.
13.
14.     // ==== e1, e2, e3 ====
15.     VectorXd e1(3);
16.     e1(0) = -x * sin(phi) / cos(phi) + k * z * cos(phi);
17.     e1(1) = x * cos(phi) /  cos(phi) + k * z * sin(phi);
18.     e1(2) = -k * x / cos(phi) + k * r1;
19.     e1.normalize();
20.
21.     VectorXd e2(3);
22.     e2(0) = cos(theta) * cos(phi);
23.     e2(1) = cos(theta) * sin(phi);
24.     e2(2) = -sin(theta);
25.
26.     VectorXd e3(3);
27.     e3(0) = sin(theta) * cos(phi);
28.     e3(1) = sin(theta) * sin(phi);
29.     e3(2) = cos(theta);
30.
31.
32.     // ==== e1, e2, e3 coefs ====
33.     double de1 = sqrt(
34.         2 * r3 * (r1 * r1 + k*k * r2*r2) / sqrt(r1*r1 + k*k*k*k * r2*r2)
35.     ); // опасно r2 <-> r3
36.     double de2 = r2;
37.     double de3 = r3;
38.
39.
40.     // ==== create preconditioning matrix ====
41.     MatrixXd res = MatrixXd::Zero(3, 3);
42.     res.col(0) = e1 * de1;
43.     res.col(1) = e2 * de2;
44.     res.col(2) = e3 * de3;
45.
46.     return res;
47. }