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 : #ifndef HAMILTONIANMATRIX_H
21 : #define HAMILTONIANMATRIX_H
22 :
23 : #include "Basisnames.hpp"
24 : #include "Serializable.hpp"
25 : #include "StateOld.hpp"
26 : #include "Symmetry.hpp"
27 : #include "utils.hpp"
28 :
29 : #include "EigenCompat.hpp"
30 : #include <Eigen/SparseCore>
31 :
32 : #include <cmath>
33 : #include <cstdint>
34 : #include <cstdio>
35 : #include <exception>
36 : #include <memory>
37 : #include <stdexcept>
38 : #include <vector>
39 :
40 : typedef double storage_double;
41 :
42 : const uint8_t csr_not_csc = 0x01; // xxx0: csc, xxx1: csr
43 : const uint8_t complex_not_real = 0x02; // xx0x: real, xx1x: complex
44 :
45 : template <typename Scalar>
46 : class Hamiltonianmatrix;
47 :
48 : template <typename Scalar>
49 : Hamiltonianmatrix<Scalar> operator+(Hamiltonianmatrix<Scalar> lhs,
50 : const Hamiltonianmatrix<Scalar> &rhs);
51 : template <typename Scalar>
52 : Hamiltonianmatrix<Scalar> operator-(Hamiltonianmatrix<Scalar> lhs,
53 : const Hamiltonianmatrix<Scalar> &rhs);
54 : template <typename Scalar, typename T>
55 : Hamiltonianmatrix<Scalar> operator*(const T &lhs, Hamiltonianmatrix<Scalar> rhs);
56 : template <typename Scalar, typename T>
57 : Hamiltonianmatrix<Scalar> operator*(Hamiltonianmatrix<Scalar> lhs, const T &rhs);
58 :
59 : template <typename Scalar>
60 : Hamiltonianmatrix<Scalar>
61 : combine(const Hamiltonianmatrix<Scalar> &lhs, const Hamiltonianmatrix<Scalar> &rhs,
62 : const double &deltaE, const std::shared_ptr<BasisnamesTwo> &basis_two, const Symmetry &sym);
63 : template <typename Scalar>
64 : void energycutoff(const Hamiltonianmatrix<Scalar> &lhs, const Hamiltonianmatrix<Scalar> &rhs,
65 : const double &deltaE, std::vector<bool> &necessary);
66 :
67 : template <typename Scalar>
68 : class Hamiltonianmatrix : public Serializable {
69 : public:
70 : Hamiltonianmatrix();
71 : Hamiltonianmatrix(const Eigen::SparseMatrix<Scalar> &entries,
72 : const Eigen::SparseMatrix<Scalar> &basis);
73 : Hamiltonianmatrix(size_t szBasis, size_t szEntries);
74 : Eigen::SparseMatrix<Scalar> &entries();
75 : const Eigen::SparseMatrix<Scalar> &entries() const;
76 : Eigen::SparseMatrix<Scalar> &basis();
77 : const Eigen::SparseMatrix<Scalar> &basis() const;
78 : size_t num_basisvectors() const;
79 : size_t num_coordinates() const;
80 : void addBasis(idx_t row, idx_t col, Scalar val);
81 : void addEntries(idx_t row, idx_t col, Scalar val);
82 : void compress(size_t nBasis, size_t nCoordinates);
83 : std::vector<Hamiltonianmatrix> findSubs() const;
84 : Hamiltonianmatrix abs() const;
85 : Hamiltonianmatrix changeBasis(const Eigen::SparseMatrix<Scalar> &basis) const;
86 : void applyCutoff(double cutoff);
87 : void findUnnecessaryStates(std::vector<bool> &isNecessaryCoordinate) const;
88 : void removeUnnecessaryBasisvectors(const std::vector<bool> &isNecessaryCoordinate);
89 : void removeUnnecessaryBasisvectors();
90 : void removeUnnecessaryStates(const std::vector<bool> &isNecessaryCoordinate);
91 : Hamiltonianmatrix getBlock(const std::vector<ptrdiff_t> &indices);
92 : void diagonalize();
93 : friend Hamiltonianmatrix operator+<>(Hamiltonianmatrix lhs, const Hamiltonianmatrix &rhs);
94 : friend Hamiltonianmatrix operator-<>(Hamiltonianmatrix lhs, const Hamiltonianmatrix &rhs);
95 : template <typename S, typename T>
96 : friend Hamiltonianmatrix<S> operator*(const T &lhs, Hamiltonianmatrix<S> rhs); // NOLINT
97 : template <typename S, typename T>
98 : friend Hamiltonianmatrix<S> operator*(Hamiltonianmatrix<S> lhs, const T &rhs); // NOLINT
99 : Hamiltonianmatrix &operator+=(const Hamiltonianmatrix &rhs);
100 : Hamiltonianmatrix &operator-=(const Hamiltonianmatrix &rhs);
101 : bytes_t &serialize() override;
102 : void doSerialization();
103 : void deserialize(bytes_t &bytesin) override;
104 : void doDeserialization();
105 : uint64_t hashEntries();
106 : uint64_t hashBasis();
107 : void save(const std::string &fname);
108 : bool load(const std::string &fname);
109 : friend Hamiltonianmatrix combine<>(const Hamiltonianmatrix &lhs, const Hamiltonianmatrix &rhs,
110 : const double &deltaE,
111 : const std::shared_ptr<BasisnamesTwo> &basis_two,
112 : const Symmetry &sym);
113 : friend void energycutoff<>(const Hamiltonianmatrix &lhs, const Hamiltonianmatrix &rhs,
114 : const double &deltaE, std::vector<bool> &necessary);
115 :
116 : template <typename T, typename std::enable_if<utils::is_complex<T>::value>::type * = nullptr>
117 0 : void mergeComplex(std::vector<storage_double> &real, std::vector<storage_double> &imag,
118 : std::vector<T> &complex) {
119 0 : std::vector<storage_double>::iterator real_it, imag_it;
120 0 : complex.reserve(real.size());
121 0 : for (real_it = real.begin(), imag_it = imag.begin(); real_it != real.end();
122 0 : ++real_it, ++imag_it) {
123 0 : complex.push_back(T(*real_it, *imag_it));
124 : }
125 0 : }
126 :
127 : template <typename T, typename std::enable_if<!utils::is_complex<T>::value>::type * = nullptr>
128 0 : void mergeComplex(std::vector<storage_double> &real, std::vector<storage_double> &imag,
129 : std::vector<T> &complex) {
130 : (void)imag;
131 0 : complex = real;
132 0 : }
133 :
134 : template <typename T, typename std::enable_if<utils::is_complex<T>::value>::type * = nullptr>
135 2 : void splitComplex(std::vector<storage_double> &real, std::vector<storage_double> &imag,
136 : std::vector<T> &complex) {
137 2 : real.reserve(complex.size());
138 2 : imag.reserve(imag.size());
139 18172 : for (auto complex_it = complex.begin(); complex_it != complex.end(); ++complex_it) {
140 18170 : real.push_back(complex_it->real());
141 18170 : imag.push_back(complex_it->imag());
142 : }
143 2 : }
144 :
145 : template <typename T, typename std::enable_if<!utils::is_complex<T>::value>::type * = nullptr>
146 52 : void splitComplex(std::vector<storage_double> &real, std::vector<storage_double> &imag,
147 : std::vector<T> &complex) {
148 52 : imag = std::vector<storage_double>();
149 52 : real = complex; // std::vector<storage_double>(complex.begin(),complex.end());
150 52 : }
151 :
152 : protected:
153 : Eigen::SparseMatrix<Scalar> entries_;
154 : Eigen::SparseMatrix<Scalar> basis_;
155 : bytes_t bytes;
156 : std::vector<Eigen::Triplet<Scalar>> triplets_basis;
157 : std::vector<Eigen::Triplet<Scalar>> triplets_entries;
158 : };
159 :
160 : extern template class Hamiltonianmatrix<std::complex<double>>;
161 : extern template Hamiltonianmatrix<std::complex<double>>
162 : operator+(Hamiltonianmatrix<std::complex<double>> lhs,
163 : const Hamiltonianmatrix<std::complex<double>> &rhs);
164 : extern template Hamiltonianmatrix<std::complex<double>>
165 : operator-(Hamiltonianmatrix<std::complex<double>> lhs,
166 : const Hamiltonianmatrix<std::complex<double>> &rhs);
167 : extern template Hamiltonianmatrix<std::complex<double>>
168 : operator*(const std::complex<double> &lhs, Hamiltonianmatrix<std::complex<double>> rhs);
169 : extern template Hamiltonianmatrix<std::complex<double>>
170 : operator*(Hamiltonianmatrix<std::complex<double>> lhs, const std::complex<double> &rhs);
171 : extern template Hamiltonianmatrix<std::complex<double>>
172 : operator*(const double &lhs, Hamiltonianmatrix<std::complex<double>> rhs);
173 : extern template Hamiltonianmatrix<std::complex<double>>
174 : operator*(Hamiltonianmatrix<std::complex<double>> lhs, const double &rhs);
175 : extern template Hamiltonianmatrix<std::complex<double>>
176 : combine(const Hamiltonianmatrix<std::complex<double>> &lhs,
177 : const Hamiltonianmatrix<std::complex<double>> &rhs, const double &deltaE,
178 : const std::shared_ptr<BasisnamesTwo> &basis_two, const Symmetry &sym);
179 : extern template void energycutoff(const Hamiltonianmatrix<std::complex<double>> &lhs,
180 : const Hamiltonianmatrix<std::complex<double>> &rhs,
181 : const double &deltaE, std::vector<bool> &necessary);
182 : extern template class Hamiltonianmatrix<double>;
183 : extern template Hamiltonianmatrix<double> operator+(Hamiltonianmatrix<double> lhs,
184 : const Hamiltonianmatrix<double> &rhs);
185 : extern template Hamiltonianmatrix<double> operator-(Hamiltonianmatrix<double> lhs,
186 : const Hamiltonianmatrix<double> &rhs);
187 : extern template Hamiltonianmatrix<double> operator*(const double &lhs,
188 : Hamiltonianmatrix<double> rhs);
189 : extern template Hamiltonianmatrix<double> operator*(Hamiltonianmatrix<double> lhs,
190 : const double &rhs);
191 : extern template Hamiltonianmatrix<double>
192 : combine(const Hamiltonianmatrix<double> &lhs, const Hamiltonianmatrix<double> &rhs,
193 : const double &deltaE, const std::shared_ptr<BasisnamesTwo> &basis_two, const Symmetry &sym);
194 : extern template void energycutoff(const Hamiltonianmatrix<double> &lhs,
195 : const Hamiltonianmatrix<double> &rhs, const double &deltaE,
196 : std::vector<bool> &necessary);
197 :
198 : #endif // HAMILTONIANMATRIX_H
|
