Line data Source code
1 : /*
2 : * Copyright (c) 2016 Sebastian Weber, Henri Menke. All rights reserved.
3 : *
4 : * This file is part of the pairinteraction library.
5 : *
6 : * The pairinteraction library is free software: you can redistribute it and/or modify
7 : * it under the terms of the GNU Lesser General Public License as published by
8 : * the Free Software Foundation, either version 3 of the License, or
9 : * (at your option) any later version.
10 : *
11 : * The pairinteraction library is distributed in the hope that it will be useful,
12 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : * GNU Lesser General Public License for more details.
15 : *
16 : * You should have received a copy of the GNU Lesser General Public License
17 : * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
18 : */
19 :
20 : #include "Hamiltonianmatrix.hpp"
21 :
22 : #include "EigenCompat.hpp"
23 : #include <Eigen/Core>
24 : #include <Eigen/Eigenvalues>
25 : #include <fmt/format.h>
26 :
27 : #include <stdexcept>
28 :
29 : template <typename Scalar>
30 35 : Hamiltonianmatrix<Scalar>::Hamiltonianmatrix() = default;
31 :
32 : template <typename Scalar>
33 4 : Hamiltonianmatrix<Scalar>::Hamiltonianmatrix(const Eigen::SparseMatrix<Scalar> &entries,
34 : const Eigen::SparseMatrix<Scalar> &basis)
35 4 : : entries_(entries), basis_(basis) {}
36 :
37 : template <typename Scalar>
38 45 : Hamiltonianmatrix<Scalar>::Hamiltonianmatrix(size_t szBasis, size_t szEntries) {
39 45 : triplets_basis.reserve(szBasis);
40 45 : triplets_entries.reserve(szEntries);
41 45 : }
42 :
43 : template <typename Scalar>
44 638 : Eigen::SparseMatrix<Scalar> &Hamiltonianmatrix<Scalar>::entries() {
45 638 : bytes.clear();
46 638 : return entries_;
47 : }
48 :
49 : template <typename Scalar>
50 12 : const Eigen::SparseMatrix<Scalar> &Hamiltonianmatrix<Scalar>::entries() const {
51 12 : return entries_;
52 : }
53 :
54 : template <typename Scalar>
55 16 : Eigen::SparseMatrix<Scalar> &Hamiltonianmatrix<Scalar>::basis() {
56 16 : bytes.clear();
57 16 : return basis_;
58 : }
59 :
60 : template <typename Scalar>
61 160546 : const Eigen::SparseMatrix<Scalar> &Hamiltonianmatrix<Scalar>::basis() const {
62 160546 : return basis_;
63 : }
64 :
65 : template <typename Scalar>
66 198 : size_t Hamiltonianmatrix<Scalar>::num_basisvectors() const {
67 198 : return basis_.cols();
68 : }
69 :
70 : template <typename Scalar>
71 4692 : size_t Hamiltonianmatrix<Scalar>::num_coordinates() const {
72 4692 : return basis_.rows();
73 : }
74 :
75 : template <typename Scalar>
76 5222 : void Hamiltonianmatrix<Scalar>::addBasis(idx_t row, idx_t col, Scalar val) {
77 5222 : triplets_basis.emplace_back(row, col, val);
78 5222 : }
79 :
80 : template <typename Scalar>
81 6318 : void Hamiltonianmatrix<Scalar>::addEntries(idx_t row, idx_t col, Scalar val) {
82 6318 : triplets_entries.emplace_back(row, col, val);
83 6318 : }
84 :
85 : template <typename Scalar>
86 45 : void Hamiltonianmatrix<Scalar>::compress(size_t nBasis, size_t nCoordinates) {
87 45 : basis_.resize(nCoordinates, nBasis);
88 45 : entries_.resize(nBasis, nBasis);
89 45 : basis_.setFromTriplets(triplets_basis.begin(), triplets_basis.end());
90 45 : entries_.setFromTriplets(triplets_entries.begin(), triplets_entries.end());
91 45 : triplets_basis.clear();
92 45 : triplets_entries.clear();
93 45 : }
94 :
95 : template <typename Scalar>
96 0 : std::vector<Hamiltonianmatrix<Scalar>> Hamiltonianmatrix<Scalar>::findSubs() const { // TODO
97 0 : std::vector<Hamiltonianmatrix<Scalar>> submatrices;
98 0 : submatrices.push_back(*this);
99 0 : return submatrices;
100 : }
101 :
102 : template <typename Scalar>
103 0 : Hamiltonianmatrix<Scalar> Hamiltonianmatrix<Scalar>::abs() const {
104 0 : return Hamiltonianmatrix<Scalar>(entries_.cwiseAbs().template cast<Scalar>(), basis_);
105 : }
106 :
107 : template <typename Scalar>
108 : Hamiltonianmatrix<Scalar>
109 4 : Hamiltonianmatrix<Scalar>::changeBasis(const Eigen::SparseMatrix<Scalar> &basis) const {
110 4 : auto transformator = basis_.adjoint() * basis;
111 4 : auto entries = transformator.adjoint() * entries_ * transformator;
112 8 : return Hamiltonianmatrix<Scalar>(entries, basis);
113 : }
114 :
115 : template <typename Scalar>
116 0 : void Hamiltonianmatrix<Scalar>::applyCutoff(double cutoff) {
117 0 : bytes.clear();
118 :
119 : // build transformator
120 0 : Eigen::VectorX<Scalar> diag = entries_.diagonal();
121 :
122 0 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
123 0 : triplets_transformator.reserve(num_basisvectors());
124 :
125 0 : size_t idxBasis = 0;
126 0 : for (size_t idx = 0; idx < this->num_basisvectors(); ++idx) {
127 0 : if (std::abs(diag[idx]) < cutoff) {
128 0 : triplets_transformator.emplace_back(idx, idxBasis++, 1);
129 : }
130 : }
131 :
132 0 : Eigen::SparseMatrix<Scalar> transformator(this->num_basisvectors(), idxBasis);
133 0 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
134 :
135 : // apply transformator
136 0 : basis_ = basis_ * transformator;
137 0 : entries_ = transformator.adjoint() * entries_ * transformator;
138 0 : }
139 :
140 : template <typename Scalar>
141 0 : void Hamiltonianmatrix<Scalar>::findUnnecessaryStates(
142 : std::vector<bool> &isNecessaryCoordinate) const {
143 0 : std::vector<double> isNecessaryCoordinate_real(num_coordinates(), 0);
144 0 : for (int k = 0; k < basis_.outerSize(); ++k) {
145 0 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basis_, k); triple;
146 0 : ++triple) {
147 0 : isNecessaryCoordinate_real[triple.row()] += std::pow(std::abs(triple.value()), 2);
148 : }
149 : }
150 :
151 0 : for (size_t idx = 0; idx < this->num_coordinates(); ++idx) {
152 0 : if (isNecessaryCoordinate_real[idx] > 0.05) { // TODO
153 0 : isNecessaryCoordinate[idx] = true;
154 : }
155 : }
156 0 : }
157 :
158 : template <typename Scalar>
159 0 : void Hamiltonianmatrix<Scalar>::removeUnnecessaryBasisvectors(
160 : const std::vector<bool> &isNecessaryCoordinate) {
161 0 : bytes.clear();
162 :
163 : // build transformator
164 0 : std::vector<double> isNecessaryBasisvector(num_basisvectors(), 0);
165 0 : for (int k_1 = 0; k_1 < basis_.outerSize(); ++k_1) {
166 0 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basis_, k_1); triple;
167 0 : ++triple) {
168 0 : ptrdiff_t col = triple.col(); // basis vector
169 0 : ptrdiff_t row = triple.row(); // coordinate
170 0 : if (isNecessaryCoordinate[row]) {
171 0 : isNecessaryBasisvector[col] += std::pow(std::abs(triple.value()), 2);
172 : }
173 : }
174 : }
175 :
176 0 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
177 0 : triplets_transformator.reserve(num_basisvectors());
178 :
179 0 : size_t idxBasis = 0;
180 0 : for (size_t idx = 0; idx < this->num_basisvectors(); ++idx) {
181 0 : if (isNecessaryBasisvector[idx] > 0.05) { // TODO
182 0 : triplets_transformator.emplace_back(idx, idxBasis++, 1);
183 : }
184 : }
185 :
186 0 : Eigen::SparseMatrix<Scalar> transformator(this->num_basisvectors(), idxBasis);
187 0 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
188 :
189 : // apply transformator
190 0 : basis_ = basis_ * transformator;
191 0 : entries_ = transformator.adjoint() * entries_ * transformator;
192 0 : }
193 :
194 : template <typename Scalar>
195 4 : void Hamiltonianmatrix<Scalar>::removeUnnecessaryBasisvectors() {
196 4 : bytes.clear();
197 :
198 : // build transformator
199 8 : std::vector<double> isNecessaryBasisvector(num_basisvectors(), 0);
200 124 : for (int k_1 = 0; k_1 < basis_.outerSize(); ++k_1) {
201 360 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basis_, k_1); triple;
202 240 : ++triple) {
203 240 : ptrdiff_t col = triple.col(); // basis vector
204 240 : isNecessaryBasisvector[col] += std::pow(std::abs(triple.value()), 2);
205 : }
206 : }
207 :
208 8 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
209 4 : triplets_transformator.reserve(num_basisvectors());
210 :
211 4 : size_t idxBasis = 0;
212 124 : for (size_t idx = 0; idx < this->num_basisvectors(); ++idx) {
213 120 : if (isNecessaryBasisvector[idx] > 0.05) {
214 120 : triplets_transformator.emplace_back(idx, idxBasis++, 1);
215 : }
216 : }
217 :
218 4 : Eigen::SparseMatrix<Scalar> transformator(this->num_basisvectors(), idxBasis);
219 4 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
220 :
221 : // apply transformator
222 4 : basis_ = basis_ * transformator;
223 4 : entries_ = transformator.adjoint() * entries_ * transformator;
224 4 : }
225 :
226 : template <typename Scalar>
227 2 : void Hamiltonianmatrix<Scalar>::removeUnnecessaryStates(
228 : const std::vector<bool> &isNecessaryCoordinate) {
229 2 : bytes.clear();
230 :
231 : // build transformator
232 4 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
233 2 : triplets_transformator.reserve(num_coordinates());
234 :
235 2 : size_t idxCoordinate = 0;
236 226 : for (size_t idx = 0; idx < this->num_coordinates(); ++idx) {
237 224 : if (isNecessaryCoordinate[idx]) {
238 224 : triplets_transformator.emplace_back(idxCoordinate++, idx, 1);
239 : }
240 : }
241 :
242 2 : Eigen::SparseMatrix<Scalar> transformator(idxCoordinate, this->num_coordinates());
243 2 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
244 :
245 : // apply transformator
246 2 : basis_ = transformator * basis_;
247 2 : }
248 :
249 : template <typename Scalar>
250 : Hamiltonianmatrix<Scalar>
251 0 : Hamiltonianmatrix<Scalar>::getBlock(const std::vector<ptrdiff_t> &indices) {
252 0 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
253 0 : triplets_transformator.reserve(indices.size());
254 0 : for (size_t idx = 0; idx < indices.size(); ++idx) {
255 0 : triplets_transformator.emplace_back(indices[idx], idx, 1);
256 : }
257 0 : Eigen::SparseMatrix<Scalar> transformator(this->num_basisvectors(), indices.size());
258 0 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
259 :
260 0 : Eigen::SparseMatrix<Scalar> block_entries = transformator.adjoint() * entries_ * transformator;
261 0 : Eigen::SparseMatrix<Scalar> block_basis = basis_ * transformator;
262 :
263 0 : return Hamiltonianmatrix<Scalar>(block_entries, block_basis);
264 : }
265 :
266 : template <typename Scalar>
267 27 : void Hamiltonianmatrix<Scalar>::diagonalize() {
268 27 : if (this->num_basisvectors() > 1) { // NOLINT
269 : // diagonalization
270 12 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixX<Scalar>> eigensolver(
271 6 : Eigen::MatrixX<Scalar>(this->entries()));
272 :
273 : // eigenvalues and eigenvectors
274 12 : Eigen::VectorX<double> evals = eigensolver.eigenvalues();
275 6 : Eigen::SparseMatrix<Scalar> evecs = eigensolver.eigenvectors().sparseView(1e-4, 0.5);
276 :
277 6 : this->entries().setZero();
278 6 : this->entries().reserve(evals.size());
279 620 : for (int idx = 0; idx < evals.size(); ++idx) {
280 614 : this->entries().insert(idx, idx) = evals.coeffRef(idx);
281 : }
282 6 : this->entries().makeCompressed();
283 :
284 6 : this->basis() = (this->basis() * evecs).pruned(1e-4, 0.5);
285 : }
286 27 : }
287 :
288 : template <typename Scalar>
289 18 : Hamiltonianmatrix<Scalar> operator+(Hamiltonianmatrix<Scalar> lhs,
290 : const Hamiltonianmatrix<Scalar> &rhs) {
291 18 : lhs.bytes.clear();
292 18 : lhs.entries_ += rhs.entries_;
293 18 : return lhs;
294 : }
295 :
296 : template <typename Scalar>
297 18 : Hamiltonianmatrix<Scalar> operator-(Hamiltonianmatrix<Scalar> lhs,
298 : const Hamiltonianmatrix<Scalar> &rhs) {
299 18 : lhs.bytes.clear();
300 18 : lhs.entries_ -= rhs.entries_;
301 18 : return lhs;
302 : }
303 :
304 : template <typename Scalar, typename T>
305 12 : Hamiltonianmatrix<Scalar> operator*(const T &lhs, Hamiltonianmatrix<Scalar> rhs) {
306 12 : rhs.bytes.clear();
307 12 : rhs.entries_ *= lhs;
308 12 : return rhs;
309 : }
310 :
311 : template <typename Scalar, typename T>
312 72 : Hamiltonianmatrix<Scalar> operator*(Hamiltonianmatrix<Scalar> lhs, const T &rhs) {
313 72 : lhs.bytes.clear();
314 72 : lhs.entries_ *= rhs;
315 72 : return lhs;
316 : }
317 :
318 : template <typename Scalar>
319 : Hamiltonianmatrix<Scalar> &
320 24 : Hamiltonianmatrix<Scalar>::operator+=(const Hamiltonianmatrix<Scalar> &rhs) {
321 24 : bytes.clear();
322 24 : entries_ += rhs.entries_;
323 24 : return *this;
324 : }
325 :
326 : template <typename Scalar>
327 : Hamiltonianmatrix<Scalar> &
328 0 : Hamiltonianmatrix<Scalar>::operator-=(const Hamiltonianmatrix<Scalar> &rhs) {
329 0 : bytes.clear();
330 0 : entries_ -= rhs.entries_;
331 0 : return *this;
332 : }
333 :
334 : template <typename Scalar>
335 0 : bytes_t &Hamiltonianmatrix<Scalar>::serialize() {
336 0 : doSerialization();
337 0 : return bytes;
338 : }
339 :
340 : template <typename Scalar>
341 27 : void Hamiltonianmatrix<Scalar>::doSerialization() {
342 27 : if (bytes.empty()) {
343 27 : entries_.makeCompressed();
344 27 : basis_.makeCompressed();
345 :
346 : // convert matrix "entries" to vectors of primitive data types
347 27 : byte_t entries_flags = 0;
348 : if (entries_.IsRowMajor != 0) {
349 : entries_flags |= csr_not_csc;
350 : }
351 : if (utils::is_complex<Scalar>::value) {
352 1 : entries_flags |= complex_not_real;
353 : }
354 27 : storage_idx_t entries_rows = entries_.rows();
355 27 : storage_idx_t entries_cols = entries_.cols();
356 81 : std::vector<Scalar> entries_data(entries_.valuePtr(),
357 27 : entries_.valuePtr() + entries_.nonZeros());
358 54 : std::vector<storage_double> entries_data_real, entries_data_imag;
359 27 : splitComplex(entries_data_real, entries_data_imag, entries_data);
360 81 : std::vector<storage_idx_t> entries_indices(entries_.innerIndexPtr(),
361 27 : entries_.innerIndexPtr() + entries_.nonZeros());
362 81 : std::vector<storage_idx_t> entries_indptr(entries_.outerIndexPtr(),
363 27 : entries_.outerIndexPtr() + entries_.outerSize());
364 :
365 : // convert matrix "basis" to vectors of primitive data types
366 27 : byte_t basis_flags = 0;
367 : if (basis_.IsRowMajor != 0) {
368 : basis_flags |= csr_not_csc;
369 : }
370 : if (utils::is_complex<Scalar>::value) {
371 1 : basis_flags |= complex_not_real;
372 : }
373 27 : storage_idx_t basis_rows = basis_.rows();
374 27 : storage_idx_t basis_cols = basis_.cols();
375 54 : std::vector<Scalar> basis_data(basis_.valuePtr(), basis_.valuePtr() + basis_.nonZeros());
376 54 : std::vector<storage_double> basis_data_real, basis_data_imag;
377 27 : splitComplex(basis_data_real, basis_data_imag, basis_data);
378 81 : std::vector<storage_idx_t> basis_indices(basis_.innerIndexPtr(),
379 27 : basis_.innerIndexPtr() + basis_.nonZeros());
380 81 : std::vector<storage_idx_t> basis_indptr(basis_.outerIndexPtr(),
381 27 : basis_.outerIndexPtr() + basis_.outerSize());
382 :
383 : // serialize vectors of primitive data types
384 54 : Serializer s;
385 27 : s << entries_flags;
386 27 : s << entries_rows;
387 27 : s << entries_cols;
388 27 : s << entries_data_real;
389 27 : if ((entries_flags & complex_not_real) != 0) {
390 1 : s << entries_data_imag;
391 : }
392 27 : s << entries_indices;
393 27 : s << entries_indptr;
394 27 : s << basis_flags;
395 27 : s << basis_rows;
396 27 : s << basis_cols;
397 27 : s << basis_data_real;
398 27 : if ((basis_flags & complex_not_real) != 0) {
399 1 : s << basis_data_imag;
400 : }
401 27 : s << basis_indices;
402 27 : s << basis_indptr;
403 27 : s.save(bytes);
404 : }
405 27 : }
406 :
407 : template <typename Scalar>
408 0 : void Hamiltonianmatrix<Scalar>::deserialize(bytes_t &bytesin) {
409 0 : bytes = bytesin;
410 0 : doDeserialization();
411 0 : }
412 :
413 : template <typename Scalar>
414 0 : void Hamiltonianmatrix<Scalar>::doDeserialization() {
415 : // deserialize vectors of primitive data types
416 : byte_t entries_flags;
417 : storage_idx_t entries_rows, entries_cols;
418 0 : std::vector<storage_double> entries_data_real, entries_data_imag;
419 0 : std::vector<idx_t> entries_indices;
420 0 : std::vector<idx_t> entries_indptr;
421 : byte_t basis_flags;
422 : storage_idx_t basis_rows, basis_cols;
423 0 : std::vector<storage_double> basis_data_real, basis_data_imag;
424 0 : std::vector<idx_t> basis_indices;
425 0 : std::vector<idx_t> basis_indptr;
426 :
427 0 : Serializer s;
428 0 : s.load(bytes);
429 0 : s >> entries_flags;
430 0 : s >> entries_rows;
431 0 : s >> entries_cols;
432 0 : s >> entries_data_real;
433 0 : if ((entries_flags & complex_not_real) != 0) {
434 0 : s >> entries_data_imag;
435 : }
436 0 : s >> entries_indices;
437 0 : s >> entries_indptr;
438 0 : s >> basis_flags;
439 0 : s >> basis_rows;
440 0 : s >> basis_cols;
441 0 : s >> basis_data_real;
442 0 : if ((basis_flags & complex_not_real) != 0) {
443 0 : s >> basis_data_imag;
444 : }
445 0 : s >> basis_indices;
446 0 : s >> basis_indptr;
447 :
448 0 : if ((((entries_flags & complex_not_real) > 0) != utils::is_complex<Scalar>::value) ||
449 0 : (((basis_flags & complex_not_real) > 0) != utils::is_complex<Scalar>::value)) {
450 0 : std::string msg("The data type used in the program does not fit the data type used in the "
451 : "serialized objects.");
452 0 : std::cout << fmt::format(">>ERR{:s}", msg.c_str()) << std::endl;
453 0 : throw std::runtime_error(msg);
454 : }
455 :
456 : // build matrix "entries_"
457 0 : std::vector<Scalar> entries_data;
458 0 : mergeComplex(entries_data_real, entries_data_imag, entries_data);
459 0 : entries_ = Eigen::SparseMatrix<Scalar>(entries_rows, entries_cols);
460 0 : entries_.makeCompressed();
461 0 : entries_.resizeNonZeros(entries_data.size());
462 0 : std::copy(entries_data.begin(), entries_data.end(), entries_.valuePtr());
463 0 : std::copy(entries_indices.begin(), entries_indices.end(), entries_.innerIndexPtr());
464 0 : std::copy(entries_indptr.begin(), entries_indptr.end(), entries_.outerIndexPtr());
465 0 : entries_.finalize();
466 :
467 : // build matrix "basis_"
468 0 : std::vector<Scalar> basis_data;
469 0 : mergeComplex(basis_data_real, basis_data_imag, basis_data);
470 0 : basis_ = Eigen::SparseMatrix<Scalar>(basis_rows, basis_cols);
471 0 : basis_.makeCompressed();
472 0 : basis_.resizeNonZeros(basis_data.size());
473 0 : std::copy(basis_data.begin(), basis_data.end(), basis_.valuePtr());
474 0 : std::copy(basis_indices.begin(), basis_indices.end(), basis_.innerIndexPtr());
475 0 : std::copy(basis_indptr.begin(), basis_indptr.end(), basis_.outerIndexPtr());
476 0 : basis_.finalize();
477 0 : }
478 :
479 : template <typename Scalar>
480 0 : uint64_t Hamiltonianmatrix<Scalar>::hashEntries() {
481 : // TODO bring this functionality to the matrix class and use it for serialization, too
482 0 : doSerialization();
483 0 : return utils::FNV64(&bytes[0], bytes.size());
484 : }
485 :
486 : template <typename Scalar>
487 0 : uint64_t Hamiltonianmatrix<Scalar>::hashBasis() {
488 : // TODO bring this functionality to the matrix class and use it for serialization, too
489 0 : doSerialization();
490 0 : return utils::FNV64(&bytes[0], bytes.size());
491 : }
492 :
493 : template <typename Scalar>
494 27 : void Hamiltonianmatrix<Scalar>::save(const std::string &fname) {
495 27 : doSerialization();
496 :
497 : // open file
498 : FILE *pFile;
499 27 : pFile = fopen(fname.c_str(), "wb");
500 :
501 : // write
502 27 : fwrite(&bytes[0], 1, sizeof(byte_t) * bytes.size(), pFile);
503 :
504 : // close file
505 27 : fclose(pFile);
506 27 : }
507 :
508 : template <typename Scalar>
509 0 : bool Hamiltonianmatrix<Scalar>::load(const std::string &fname) {
510 : try {
511 : // open file
512 0 : if (FILE *pFile = fopen(fname.c_str(), "rb")) {
513 : // obtain file size:
514 0 : fseek(pFile, 0, SEEK_END);
515 0 : size_t size_file = ftell(pFile);
516 0 : rewind(pFile);
517 :
518 : // read
519 0 : bytes.resize(size_file / sizeof(byte_t));
520 0 : size_t size_result = fread(&bytes[0], 1, sizeof(byte_t) * bytes.size(), pFile);
521 0 : if (size_result != size_file) {
522 0 : throw std::runtime_error("Matrix could not be read from file.");
523 : }
524 :
525 : // close file
526 0 : fclose(pFile);
527 :
528 0 : doDeserialization();
529 :
530 0 : return true;
531 : }
532 0 : return false;
533 :
534 0 : } catch (std::exception &e) {
535 0 : #pragma omp critical(textoutput)
536 0 : std::cerr << e.what() << std::endl;
537 0 : return false;
538 : }
539 : }
540 :
541 : template <typename Scalar>
542 4 : Hamiltonianmatrix<Scalar> combine(const Hamiltonianmatrix<Scalar> &lhs,
543 : const Hamiltonianmatrix<Scalar> &rhs, const double &deltaE,
544 : const std::shared_ptr<BasisnamesTwo> &basis_two,
545 : const Symmetry &sym) {
546 : // TODO program a faster method for samebasis == true
547 :
548 4 : size_t num_basisvectors = lhs.num_basisvectors() * rhs.num_basisvectors();
549 4 : size_t num_coordinates = lhs.num_coordinates() * rhs.num_coordinates();
550 :
551 : ////////////////////////////////////////////////////////
552 : ////// Mapping used in case of reflection symmetry /////
553 : ////////////////////////////////////////////////////////
554 :
555 8 : std::vector<size_t> mapping(num_coordinates, -1);
556 4 : if (sym.reflection != NA) { // NOLINT
557 0 : std::unordered_map<StateTwoOld, size_t> buffer;
558 0 : for (auto state : *basis_two) {
559 0 : if (state.m[0] < 0) {
560 0 : continue;
561 : }
562 0 : state.m[0] *= -1;
563 0 : state.m[1] *= -1;
564 0 : buffer[state] = state.idx;
565 : }
566 0 : for (auto state : *basis_two) {
567 0 : if (state.m[0] > 0) {
568 0 : continue;
569 : }
570 0 : mapping[buffer[state]] = state.idx;
571 0 : if (sym.inversion != NA || sym.permutation != NA) {
572 0 : mapping[state.idx] = buffer[state];
573 : }
574 : }
575 : }
576 :
577 : ////////////////////////////////////////////////////////
578 : ////// Combine basis and entries ///////////////////////
579 : ////////////////////////////////////////////////////////
580 :
581 : // This approach does only work if lhs.entries() and rhs.entries() are diagonal // TODO assert
582 :
583 : // === Initialize matrices ===
584 8 : Eigen::VectorX<Scalar> diag1 = lhs.entries().diagonal();
585 8 : Eigen::VectorX<Scalar> diag2 = rhs.entries().diagonal();
586 :
587 : // Number of elements for which space sould be reserved
588 4 : size_t size_basis = num_basisvectors; // TODO estimate better
589 4 : size_t size_entries = num_basisvectors;
590 :
591 4 : Hamiltonianmatrix<Scalar> mat(size_basis, size_entries);
592 :
593 : // === Combine basis and entries ===
594 4 : size_t col = 0; // basis vector
595 452 : for (int col_1 = 0; col_1 < lhs.basis().outerSize();
596 : ++col_1) { // outerSize() == num_cols = num_basisvectors()
597 100800 : for (int col_2 = 0; col_2 < rhs.basis().outerSize(); ++col_2) {
598 :
599 : // In case of inversion symmetry: skip half of the basis vector pairs
600 100352 : if ((sym.inversion == EVEN && col_1 <= col_2) || // gerade
601 75152 : (sym.inversion == ODD && col_1 < col_2)) { // ungerade
602 50176 : continue;
603 : }
604 :
605 : // In case of permutation symmetry: skip half of the basis vector pairs
606 50176 : if ((sym.permutation == EVEN && col_1 <= col_2) || // sym
607 49952 : (sym.permutation == ODD && col_1 < col_2)) { // asym
608 224 : continue;
609 : }
610 :
611 : // TODO combine inversion and permutation symmetry (one variable "permuinversion"),
612 : // think about further simplifications
613 :
614 : // --- Combine diagonal elements for mat.entries() ---
615 49952 : Scalar val_entries = diag1[col_1] + diag2[col_2]; // diag(V) x I + I x diag(V)
616 :
617 : // Check whether the new diagonal element is within the energy cutoff
618 97742 : if (std::abs(val_entries) < deltaE + 1e-11 ||
619 47790 : deltaE < 0) { // TODO avoid the "+1e-11" hack
620 :
621 : // Variable that stores whether the combindes basis vector, that belongs to the
622 : // combined diagonal element, is valid
623 2162 : bool existing = false;
624 :
625 : // --- Combine basis vectors for mat.basis() ---
626 4324 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple_1(lhs.basis(),
627 : col_1);
628 2162 : triple_1; ++triple_1) {
629 4324 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple_2(rhs.basis(),
630 : col_2);
631 2162 : triple_2; ++triple_2) {
632 2162 : size_t row =
633 2162 : rhs.num_coordinates() * triple_1.row() + triple_2.row(); // coordinate
634 :
635 : // Get pair state that belongs to the current coordinate of the combined
636 : // basis vector
637 2162 : const auto &state = basis_two->get(row);
638 :
639 2162 : float M = state.m[0] + state.m[1];
640 2162 : int parityL = std::pow(-1, state.l[0] + state.l[1]);
641 2162 : int parityJ = std::pow(-1, state.j[0] + state.j[1]);
642 2162 : int parityM = std::pow(-1, state.m[0] + state.m[1]);
643 :
644 : // In case of inversion and reflection symmetry: check whether the inversion
645 : // symmetric state is already reflection symmetric
646 2162 : bool skip_reflection = false;
647 2162 : if (sym.inversion != NA && col_1 != col_2 && sym.reflection != NA &&
648 0 : mapping[row] ==
649 0 : rhs.num_coordinates() * triple_2.row() + triple_1.row()) {
650 0 : if (((sym.inversion == EVEN) ? -parityL : parityL) !=
651 0 : ((sym.reflection == EVEN) ? parityL * parityJ * parityM
652 0 : : -parityL * parityJ * parityM)) {
653 0 : continue; // the parity under inversion and reflection is different
654 : }
655 0 : skip_reflection =
656 : true; // the parity under inversion and reflection is the same
657 : }
658 :
659 : // In case of permutation and reflection symmetry: check whether the
660 : // permutation symmetric state is already reflection symmetric
661 2162 : if (sym.permutation != NA && col_1 != col_2 && sym.reflection != NA &&
662 0 : mapping[row] ==
663 0 : rhs.num_coordinates() * triple_2.row() + triple_1.row()) {
664 0 : if (((sym.permutation == EVEN) ? -1 : 1) !=
665 0 : ((sym.reflection == EVEN) ? parityL * parityJ * parityM
666 0 : : -parityL * parityJ * parityM)) {
667 0 : continue; // the parity under permutation and reflection is
668 : // different
669 : }
670 0 : skip_reflection =
671 : true; // the parity under permutation and reflection is the same
672 : }
673 :
674 : // In case of inversion and permutation symmetry: the inversion symmetric
675 : // state is already permutation symmetric
676 2162 : bool skip_permutation = false;
677 2162 : if (sym.inversion != NA && sym.permutation != NA && col_1 != col_2) {
678 4324 : if (((sym.inversion == EVEN) ? -parityL : parityL) !=
679 2162 : ((sym.permutation == EVEN) ? -1 : 1)) {
680 1202 : continue; // the parity under inversion and permutation is different
681 : }
682 960 : skip_permutation =
683 : true; // the parity under inversion and permutation is the same
684 : }
685 :
686 : // In case of rotation symmetry: skip coordinates with wrong total magnetic
687 : // momentum
688 960 : if (sym.rotation != NA && sym.rotation != M &&
689 840 : (sym.reflection == NA || sym.rotation != -M)) {
690 840 : continue;
691 : }
692 :
693 : // In case of reflection symmetry: skip half of the coordinates
694 120 : if (sym.reflection != NA && state.m[0] < 0 && !skip_reflection) {
695 0 : continue;
696 : }
697 :
698 : // Calculate coefficient that belongs to the current coordinate
699 120 : Scalar val_basis = triple_1.value() * triple_2.value(); // coefficient
700 120 : if (sym.reflection != NA && !skip_reflection) {
701 0 : val_basis /= std::sqrt(2);
702 : }
703 120 : if (sym.inversion != NA && col_1 != col_2) {
704 120 : val_basis /= std::sqrt(2);
705 : }
706 120 : if (sym.permutation != NA && col_1 != col_2 && !skip_permutation) {
707 0 : val_basis /= std::sqrt(2);
708 : }
709 :
710 : // Save the coefficient taking into account the symmetrization
711 120 : mat.addBasis(row, col, val_basis);
712 :
713 120 : if (sym.reflection != NA && !skip_reflection) {
714 0 : size_t r = mapping[row];
715 0 : Scalar v = val_basis;
716 0 : v *= (sym.reflection == EVEN) ? parityL * parityJ * parityM
717 0 : : -parityL * parityJ * parityM;
718 0 : mat.addBasis(r, col, v);
719 : }
720 :
721 120 : if (sym.inversion != NA && col_1 != col_2) {
722 120 : size_t r = rhs.num_coordinates() * triple_2.row() + triple_1.row();
723 120 : Scalar v = val_basis;
724 120 : v *= (sym.inversion == EVEN) ? -parityL : parityL;
725 120 : mat.addBasis(r, col, v);
726 : }
727 :
728 120 : if (sym.inversion != NA && col_1 != col_2 && sym.reflection != NA &&
729 0 : !skip_reflection) {
730 0 : size_t r = rhs.num_coordinates() * triple_2.row() + triple_1.row();
731 0 : r = mapping[r];
732 0 : Scalar v = val_basis;
733 0 : v *= (sym.reflection == EVEN) ? parityL * parityJ * parityM
734 0 : : -parityL * parityJ * parityM;
735 0 : v *= (sym.inversion == EVEN) ? -parityL : parityL;
736 0 : mat.addBasis(r, col, v);
737 : }
738 :
739 120 : if (sym.permutation != NA && col_1 != col_2 && !skip_permutation) {
740 0 : size_t r = rhs.num_coordinates() * triple_2.row() + triple_1.row();
741 0 : Scalar v = val_basis;
742 0 : v *= (sym.permutation == EVEN) ? -1 : 1;
743 0 : mat.addBasis(r, col, v);
744 : }
745 :
746 120 : if (sym.permutation != NA && col_1 != col_2 && !skip_permutation &&
747 0 : sym.reflection != NA && !skip_reflection) {
748 0 : size_t r = rhs.num_coordinates() * triple_2.row() + triple_1.row();
749 0 : r = mapping[r];
750 0 : Scalar v = val_basis;
751 0 : v *= (sym.reflection == EVEN) ? parityL * parityJ * parityM
752 0 : : -parityL * parityJ * parityM;
753 0 : v *= (sym.permutation == EVEN) ? -1 : 1;
754 0 : mat.addBasis(r, col, v);
755 : }
756 :
757 120 : existing = true;
758 : }
759 : }
760 :
761 : // Save the combined diagonal element if the corresponding combined basis vector is
762 : // valid
763 2162 : if (existing) {
764 120 : mat.addEntries(col, col, val_entries);
765 120 : ++col;
766 : }
767 : }
768 : }
769 : }
770 :
771 : // Finalize the Hamiltonian matrix
772 4 : num_basisvectors = col;
773 4 : mat.compress(num_basisvectors, num_coordinates);
774 8 : return mat;
775 : }
776 :
777 : template <typename Scalar>
778 2 : void energycutoff(const Hamiltonianmatrix<Scalar> &lhs, const Hamiltonianmatrix<Scalar> &rhs,
779 : const double &deltaE, std::vector<bool> &necessary) {
780 4 : Eigen::VectorX<Scalar> diag1 = lhs.entries().diagonal();
781 4 : Eigen::VectorX<Scalar> diag2 = rhs.entries().diagonal();
782 :
783 226 : for (int col_1 = 0; col_1 < lhs.basis().outerSize();
784 : ++col_1) { // outerSize() == num_cols = num_basisvectors()
785 50400 : for (int col_2 = 0; col_2 < rhs.basis().outerSize(); ++col_2) {
786 50176 : Scalar val_entries = diag1[col_1] + diag2[col_2]; // diag(V) x I + I x diag(V)
787 98180 : if (std::abs(val_entries) < deltaE + 1e-11 ||
788 48004 : deltaE < 0) { // TODO make +1e-11 unnecessary
789 4344 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple_1(lhs.basis(),
790 : col_1);
791 2172 : triple_1; ++triple_1) {
792 4344 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple_2(rhs.basis(),
793 : col_2);
794 2172 : triple_2; ++triple_2) {
795 2172 : size_t row =
796 2172 : rhs.num_coordinates() * triple_1.row() + triple_2.row(); // coordinate
797 2172 : necessary[row] = true;
798 : }
799 : }
800 : }
801 : }
802 : }
803 2 : }
804 :
805 : template class Hamiltonianmatrix<std::complex<double>>;
806 : template Hamiltonianmatrix<std::complex<double>>
807 : operator+(Hamiltonianmatrix<std::complex<double>> lhs,
808 : const Hamiltonianmatrix<std::complex<double>> &rhs);
809 : template Hamiltonianmatrix<std::complex<double>>
810 : operator-(Hamiltonianmatrix<std::complex<double>> lhs,
811 : const Hamiltonianmatrix<std::complex<double>> &rhs);
812 : template Hamiltonianmatrix<std::complex<double>>
813 : operator*(const std::complex<double> &lhs, Hamiltonianmatrix<std::complex<double>> rhs);
814 : template Hamiltonianmatrix<std::complex<double>>
815 : operator*(Hamiltonianmatrix<std::complex<double>> lhs, const std::complex<double> &rhs);
816 : template Hamiltonianmatrix<std::complex<double>>
817 : operator*(const double &lhs, Hamiltonianmatrix<std::complex<double>> rhs);
818 : template Hamiltonianmatrix<std::complex<double>>
819 : operator*(Hamiltonianmatrix<std::complex<double>> lhs, const double &rhs);
820 : template Hamiltonianmatrix<std::complex<double>>
821 : combine(const Hamiltonianmatrix<std::complex<double>> &lhs,
822 : const Hamiltonianmatrix<std::complex<double>> &rhs, const double &deltaE,
823 : const std::shared_ptr<BasisnamesTwo> &basis_two, const Symmetry &sym);
824 : template void energycutoff(const Hamiltonianmatrix<std::complex<double>> &lhs,
825 : const Hamiltonianmatrix<std::complex<double>> &rhs, const double &deltaE,
826 : std::vector<bool> &necessary);
827 : template class Hamiltonianmatrix<double>;
828 : template Hamiltonianmatrix<double> operator+(Hamiltonianmatrix<double> lhs,
829 : const Hamiltonianmatrix<double> &rhs);
830 : template Hamiltonianmatrix<double> operator-(Hamiltonianmatrix<double> lhs,
831 : const Hamiltonianmatrix<double> &rhs);
832 : template Hamiltonianmatrix<double> operator*(const double &lhs, Hamiltonianmatrix<double> rhs);
833 : template Hamiltonianmatrix<double> operator*(Hamiltonianmatrix<double> lhs, const double &rhs);
834 : template Hamiltonianmatrix<double>
835 : combine(const Hamiltonianmatrix<double> &lhs, const Hamiltonianmatrix<double> &rhs,
836 : const double &deltaE, const std::shared_ptr<BasisnamesTwo> &basis_two, const Symmetry &sym);
837 : template void energycutoff(const Hamiltonianmatrix<double> &lhs,
838 : const Hamiltonianmatrix<double> &rhs, const double &deltaE,
839 : std::vector<bool> &necessary);
|
