Documentation of CSL
libdiagonalization_cppdata.h
Go to the documentation of this file.
1 // This file is part of MARTY.
2 //
3 // MARTY is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 3 of the License, or
6 // (at your option) any later version.
7 //
8 // MARTY is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with MARTY. If not, see <https://www.gnu.org/licenses/>.
15 
23 #pragma once
24 
25 #include <iostream>
26 
27 namespace csl {
28 
29 void print_libdiagonalization_cppdata(std::ostream &out) {
30  out << "#include <cassert>\n";
31  out << "#include \"libdiagonalization.h\"\n";
32  out << "#include <gsl/gsl_complex.h>\n";
33  out << "#include <gsl/gsl_eigen.h>\n";
34  out << "#include <gsl/gsl_errno.h>\n";
35  out << "#include <gsl/gsl_complex_math.h>\n";
36  out << "#include <iostream>\n";
37  out << "\n";
38  out << "#ifdef QUAD\n";
39  out << " #include <quadmath.h>\n";
40  out << " #define ABS fabsq\n";
41  out << " #define SQRT sqrtq\n";
42  out << " #define REAL(arg) crealq(arg)\n";
43  out << " #define IMAG(arg) cimagq(arg)\n";
44  out << " #define CONJ(arg) conjq(arg)\n";
45  out << "#else\n";
46  out << " #define ABS std::abs\n";
47  out << " #define SQRT std::sqrt\n";
48  out << " #define REAL(arg) (arg).real()\n";
49  out << " #define IMAG(arg) (arg).imag()\n";
50  out << " #define CONJ(arg) std::conj(arg)\n";
51  out << "#endif\n";
52  out << "\n";
53  out << "Diagonalizer::real Diagonalizer::precision = 1e-5;\n";
54  out << "\n";
55  out << "void Diagonalizer::applyDiagonalization(\n";
56  out << " std::initializer_list<complex> massMatrix,\n";
57  out << " std::initializer_list<complex*> transfer,\n";
58  out << " std::initializer_list<real*> masses\n";
59  out << " )\n";
60  out << "{\n";
61  out << " Diagonalizer diag(massMatrix);\n";
62  out << "\n";
63  out << " std::vector<complex> const &t = diag.transfer();\n";
64  out << " size_t i = 0;\n";
65  out << " for (auto iter = transfer.begin(); iter != transfer.end(); ++iter)\n";
66  out << " **iter = t[i++];\n";
67  out << "\n";
68  out << " i = 0;\n";
69  out << " std::vector<real> const &m = diag.mass();\n";
70  out << " for (auto iter = masses.begin(); iter != masses.end(); ++iter)\n";
71  out << " **iter = m[i++];\n";
72  out << "}\n";
73  out << "\n";
74  out << "void Diagonalizer::applyBiDiagonalization(\n";
75  out << " std::initializer_list<complex> massMatrix,\n";
76  out << " std::initializer_list<complex*> transfer,\n";
77  out << " std::initializer_list<complex*> transfer2,\n";
78  out << " std::initializer_list<real*> masses\n";
79  out << " )\n";
80  out << "{\n";
81  out << " Diagonalizer diag(massMatrix, true);\n";
82  out << "\n";
83  out << " std::vector<complex> const &t = diag.transfer();\n";
84  out << " size_t i = 0;\n";
85  out << " for (auto iter = transfer.begin(); iter != transfer.end(); ++iter)\n";
86  out << " **iter = t[i++];\n";
87  out << "\n";
88  out << " i = 0;\n";
89  out << " std::vector<complex> const &t2 = diag.transfer2();\n";
90  out << " for (auto iter = transfer2.begin(); iter != transfer2.end(); ++iter)\n";
91  out << " **iter = t2[i++];\n";
92  out << "\n";
93  out << " i = 0;\n";
94  out << " std::vector<real> const &m = diag.mass();\n";
95  out << " for (auto iter = masses.begin(); iter != masses.end(); ++iter)\n";
96  out << " **iter = m[i++];\n";
97  out << "}\n";
98  out << "\n";
99  out << "Diagonalizer::Diagonalizer(\n";
100  out << " std::initializer_list<complex> massMatrix,\n";
101  out << " bool t_bidiag\n";
102  out << " )\n";
103  out << " :m_data(massMatrix),\n";
104  out << " m_N(std::sqrt(m_data.size())),\n";
105  out << " m_bidiag(t_bidiag),\n";
106  out << " m_computed(false)\n";
107  out << "{\n";
108  out << " assert(m_N * m_N == massMatrix.size());\n";
109  out << "}\n";
110  out << "\n";
111  out << "std::vector<Diagonalizer::complex> const &Diagonalizer::transfer()\n";
112  out << "{\n";
113  out << " updateDiagonalization();\n";
114  out << " return m_transfer;\n";
115  out << "}\n";
116  out << "\n";
117  out << "std::vector<Diagonalizer::complex> const &Diagonalizer::transfer2()\n";
118  out << "{\n";
119  out << " updateDiagonalization();\n";
120  out << " assert(m_bidiag);\n";
121  out << " return m_transfer2;\n";
122  out << "}\n";
123  out << "\n";
124  out << "std::vector<Diagonalizer::real> const &\n";
125  out << " Diagonalizer::mass()\n";
126  out << "{\n";
127  out << " updateDiagonalization();\n";
128  out << " return m_mass;\n";
129  out << "}\n";
130  out << "\n";
131  out << "\n";
132  out << "Diagonalizer::real Diagonalizer::det(\n";
133  out << " std::vector<Diagonalizer::complex> const &mat\n";
134  out << " ) const\n";
135  out << "{\n";
136  out << " return REAL(\n";
137  out << " mat[index(0, 0)]*mat[index(1, 1)]\n";
138  out << " - mat[index(1, 0)]*mat[index(0, 1)]\n";
139  out << " );\n";
140  out << "}\n";
141  out << "\n";
142  out << "void Diagonalizer::updateDiagonalization()\n";
143  out << "{\n";
144  out << " if (!m_computed) {\n";
145  out << " (m_bidiag) ? bidiagonalize() : diagonalize();\n";
146  out << " m_computed = true;\n";
147  out << " }\n";
148  out << "}\n";
149  out << "\n";
150  out << "void Diagonalizer::diagonalize(\n";
151  out << " std::vector<complex> massMatrix,\n";
152  out << " std::vector<complex> &transfer,\n";
153  out << " std::vector<real> &mass,\n";
154  out << " bool bidiagonalization\n";
155  out << " )\n";
156  out << "{\n";
157  out << " transfer = std::vector<complex>(m_N*m_N, 0);\n";
158  out << " mass = std::vector<real>(m_N, 0);\n";
159  out << " gsl_eigen_hermv_workspace* workspace = gsl_eigen_hermv_alloc(m_N);\n";
160  out << "\n";
161  out << " if (!bidiagonalization) {\n";
162  out << " // Symmetrize the mass matrix (if it has not been done yet)\n";
163  out << " const size_t N = std::sqrt(massMatrix.size());\n";
164  out << " for (size_t i = 0; i != N; ++i)\n";
165  out << " for (size_t j = i+1; j < N; ++j) {\n";
166  out << " massMatrix[index(i, j)] = \n";
167  out << " (massMatrix[index(i, j)] + massMatrix[index(j, i)]) / 2.;\n";
168  out << " massMatrix[index(j, i)] = massMatrix[index(i, j)];\n";
169  out << " }\n";
170  out << " }\n";
171  out << " gsl_matrix_complex* init_gsl = getGSLMassMatrix(massMatrix);\n";
172  out << " gsl_vector* eigenValues = gsl_vector_alloc(m_N);\n";
173  out << " gsl_matrix_complex* transfer_gsl = gsl_matrix_complex_alloc(m_N, m_N);\n";
174  out << "\n";
175  out << " // Diagonalization\n";
176  out << " gsl_eigen_hermv(init_gsl, eigenValues, transfer_gsl, workspace);\n";
177  out << "\n";
178  out << " // Gathering results\n";
179  out << " loadGSLMatrix(transfer_gsl, transfer);\n";
180  out << " loadGSLVector(eigenValues, mass);\n";
181  out << "\n";
182  out << " // Freeing allocated memory by gsl\n";
183  out << " gsl_matrix_complex_free(transfer_gsl);\n";
184  out << " gsl_vector_free (eigenValues);\n";
185  out << " gsl_matrix_complex_free(init_gsl);\n";
186  out << " gsl_eigen_hermv_free (workspace);\n";
187  out << "}\n";
188  out << "\n";
189  out << "void Diagonalizer::swap(\n";
190  out << " const size_t i,\n";
191  out << " const size_t j,\n";
192  out << " std::vector<complex> &transfer,\n";
193  out << " std::vector<real> &mass\n";
194  out << " )\n";
195  out << "{\n";
196  out << " const size_t N = mass.size();\n";
197  out << " for (size_t k = 0; k != N; ++k) {\n";
198  out << " std::swap(transfer[index(k, i)], transfer[index(k, j)]);\n";
199  out << " }\n";
200  out << "}\n";
201  out << "\n";
202  out << "void Diagonalizer::sort(\n";
203  out << " std::vector<complex> &transfer,\n";
204  out << " std::vector<real> &mass\n";
205  out << " )\n";
206  out << "{\n";
207  out << " for (size_t i = 0; i != mass.size(); ++i) {\n";
208  out << " size_t i_min = i;\n";
209  out << " for (size_t j = i+1; j < mass.size(); ++j) {\n";
210  out << " if (ABS(mass[j]) < ABS(mass[i_min]))\n";
211  out << " i_min = j;\n";
212  out << " }\n";
213  out << " if (i != i_min) {\n";
214  out << " swap(i, i_min, transfer, mass);\n";
215  out << " std::swap(mass[i_min], mass[i]);\n";
216  out << " }\n";
217  out << " }\n";
218  out << "}\n";
219  out << "\n";
220  out << "void Diagonalizer::sort(\n";
221  out << " std::vector<complex> &transfer,\n";
222  out << " std::vector<complex> &transfer2,\n";
223  out << " std::vector<real> &mass\n";
224  out << " )\n";
225  out << "{\n";
226  out << " for (size_t i = 0; i != mass.size(); ++i) {\n";
227  out << " size_t i_min = i;\n";
228  out << " for (size_t j = i+1; j < mass.size(); ++j) {\n";
229  out << " if (ABS(mass[j]) < ABS(mass[i_min]))\n";
230  out << " i_min = j;\n";
231  out << " }\n";
232  out << " if (i != i_min) {\n";
233  out << " swap(i, i_min, transfer, mass);\n";
234  out << " swap(i, i_min, transfer2, mass);\n";
235  out << " std::swap(mass[i_min], mass[i]);\n";
236  out << " }\n";
237  out << " }\n";
238  out << "}\n";
239  out << "\n";
240  out << "void Diagonalizer::diagonalize()\n";
241  out << "{\n";
242  out << " diagonalize(m_data, m_transfer, m_mass);\n";
243  out << " sort(m_transfer, m_mass);\n";
244  out << "}\n";
245  out << "\n";
246  out << "void Diagonalizer::bidiagonalize()\n";
247  out << "{\n";
248  out << " std::vector<complex> m_data_herm = hermitian(m_data);\n";
249  out << " std::vector<complex> m_m_dagger = dot(m_data, m_data_herm);\n";
250  out << " std::vector<complex> m_dagger_m = dot(m_data_herm, m_data);\n";
251  out << "\n";
252  out << " std::vector<real> mass2 = m_mass;\n";
253  out << " diagonalize(m_m_dagger, m_transfer2, m_mass, true);\n";
254  out << " diagonalize(m_dagger_m, m_transfer, mass2, true);\n";
255  out << " sort(m_transfer2, m_mass);\n";
256  out << " sort(m_transfer, mass2);\n";
257  out << "\n";
258  out << " for (size_t i = 0; i != m_N; ++i) {\n";
259  out << " assert(ABS(m_mass[i] - mass2[i]) / ((m_mass[i] == 0) ? 1 : m_mass[i])\n";
260  out << " < precision);\n";
261  out << " m_mass[i] = SQRT(m_mass[i]);\n";
262  out << " }\n";
263  out << "\n";
264  out << " auto diagMass = dot(dot(hermitian(m_transfer2), m_data), m_transfer);\n";
265  out << " for (size_t i = 0; i != m_N; ++i) \n";
266  out << " if (REAL(diagMass[index(i, i)]*m_mass[i]) < 0) \n";
267  out << " for (size_t j = 0; j != m_N; ++j)\n";
268  out << " m_transfer[index(j, i)] *= -1;\n";
269  out << "\n";
270  out << " for (auto &el : m_transfer2)\n";
271  out << " el = CONJ(el);\n";
272  out << "}\n";
273  out << "\n";
274  out << "gsl_matrix_complex *Diagonalizer::getGSLMassMatrix(\n";
275  out << " std::vector<complex> const &massMatrix\n";
276  out << " ) const \n";
277  out << "{\n";
278  out << " gsl_matrix_complex *matrix = gsl_matrix_complex_alloc(m_N, m_N);\n";
279  out << " for (size_t i = 0; i != m_N; ++i)\n";
280  out << " for (size_t j = 0; j != m_N; ++j)\n";
281  out << " gsl_matrix_complex_set(\n";
282  out << " matrix, \n";
283  out << " i, j, \n";
284  out << " gsl_complex_rect(REAL(massMatrix[index(i, j)]),\n";
285  out << " IMAG(massMatrix[index(i, j)]))\n";
286  out << " );\n";
287  out << "\n";
288  out << " return matrix;\n";
289  out << "}\n";
290  out << "\n";
291  out << "void Diagonalizer::loadGSLMatrix(\n";
292  out << " gsl_matrix_complex const *matrix,\n";
293  out << " std::vector<complex> &target\n";
294  out << " )\n";
295  out << "{\n";
296  out << " target = std::vector<complex>(m_N * m_N);\n";
297  out << " for (size_t i = 0; i != m_N; ++i)\n";
298  out << " for (size_t j = 0; j != m_N; ++j) {\n";
299  out << " auto gslcomplex = gsl_matrix_complex_get(matrix, i, j);\n";
300  out << " target[index(i, j)] = {GSL_REAL(gslcomplex),\n";
301  out << " GSL_IMAG(gslcomplex)};\n";
302  out << " }\n";
303  out << "}\n";
304  out << "\n";
305  out << "void Diagonalizer::loadGSLVector(\n";
306  out << " gsl_vector const *vect,\n";
307  out << " std::vector<real> &target\n";
308  out << " )\n";
309  out << "{\n";
310  out << " target = std::vector<real>(m_N);\n";
311  out << " for (size_t i = 0; i != m_N; ++i) \n";
312  out << " target[i] = gsl_vector_get(vect, i);\n";
313  out << "}\n";
314  out << "\n";
315  out << "void Diagonalizer::positiveDiagonal(\n";
316  out << " std::vector<complex> &transfer\n";
317  out << " )\n";
318  out << "{\n";
319  out << " for (size_t j = 0; j != m_N; ++j) {\n";
320  out << " if (REAL(transfer[index(j, j)]) < 0) {\n";
321  out << " for (size_t i = 0; i != m_N; ++i)\n";
322  out << " transfer[index(i, j)] *= -1;\n";
323  out << " }\n";
324  out << " }\n";
325  out << "}\n";
326  out << "\n";
327  out << "std::vector<Diagonalizer::complex> Diagonalizer::hermitian(\n";
328  out << " std::vector<complex> const &init\n";
329  out << " ) const\n";
330  out << "{\n";
331  out << " std::vector<complex> res(init);\n";
332  out << " for (size_t i = 0; i != m_N; ++i)\n";
333  out << " for (size_t j = 0; j != m_N; ++j)\n";
334  out << " res[index(i, j)] = CONJ(init[index(j, i)]);\n";
335  out << " return res;\n";
336  out << "}\n";
337  out << "\n";
338  out << "std::vector<Diagonalizer::complex> Diagonalizer::dot(\n";
339  out << " std::vector<complex> const &A,\n";
340  out << " std::vector<complex> const &B\n";
341  out << " ) const\n";
342  out << "{\n";
343  out << " std::vector<complex> C(A);\n";
344  out << " for (size_t i = 0; i != m_N; ++i)\n";
345  out << " for (size_t j = 0; j != m_N; ++j) {\n";
346  out << " C[index(i, j)] = 0;\n";
347  out << " for (size_t k = 0; k != m_N; ++k)\n";
348  out << " C[index(i, j)] += A[index(i, k)] * B[index(k, j)];\n";
349  out << " }\n";
350  out << " return C;\n";
351  out << "}\n";
352 }
353 
354 } // End of namespace csl
csl
Namespace for csl library.
Definition: abreviation.h:34
Generated by   doxygen 1.8.13