MatthiasQuintern/cpp-mps
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add dmrg

Browse Source
This commit is contained in:
matthias@quintern.xyz 2025-01-26 19:58:40 +01:00
parent cc805c0166
commit 8fcb948ecc
33 changed files with 1322 additions and 1362 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
.gitignore vendored Normal file
View File

@ -0,0 +1,17 @@
*.pdf
*.html
*.js
*benchmarks
*build*
docs
include
*.so
*.py
mps
src-alloc-debug/
src-eigen-opt/
src/.old
src/pch.hpp.gch
test-python/
web/

6
.gitmodules vendored Normal file
View File

@ -0,0 +1,6 @@
[submodule "lambda-lanczos"]
path = lambda-lanczos
url = https://github.com/mrcdr/lambda-lanczos
[submodule "src/spectra"]
path = spectra
url = https://github.com/yixuan/spectra

1
lambda-lanczos Submodule

@ -0,0 +1 @@
Subproject commit d45baddc4239dd1de3f0091ce495920b57a3fd5e

1
spectra Submodule

@ -0,0 +1 @@
Subproject commit a29f37fe3ed1c2895b07b449fcb8f09bc0940e40

2
src/.clangd
View File

@ -1,3 +1,3 @@
CompileFlags: # Tweak the parse settings
Add: [ -std=c++23, -Wall, -Wpedantic, -IFastor, -I., -I.., -I../include, -Imodels, -I../models, -Ialgorithms, -DEMSCRIPTEN ]
Add: [ -std=c++23, -Wall, -Wpedantic, -IFastor, -I., -I.., -Imodels, -I../models, -Ialgorithms, -iprefix ~/Projekte/science/cpp-mps/ -iwithprefix include -iwithprefix lambda-lanczos/include -iwithprefix spectra/include -DEMSCRIPTEN ]
# https://clangd.llvm.org/config

4
src/.doxygen_config
View File

@ -965,7 +965,7 @@ WARN_LOGFILE =
# spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
# Note: If this tag is empty the current directory is searched.
INPUT = .
INPUT = . models algorithms wasm
# This tag can be used to specify the character encoding of the source files
# that Doxygen parses. Internally Doxygen uses the UTF-8 encoding. Doxygen uses
@ -2087,7 +2087,7 @@ PAPER_TYPE = a4
# If left blank no extra packages will be included.
# This tag requires that the tag GENERATE_LATEX is set to YES.
EXTRA_PACKAGES = {amsmath}
EXTRA_PACKAGES = {amsmath} {bbm}
# The LATEX_HEADER tag can be used to specify a user-defined LaTeX header for
# the generated LaTeX document. The header should contain everything until the

70
src/Makefile
View File

@ -9,19 +9,34 @@ WEBEXEC = ~/Projekte/web/html-quintern/src/script/mps/mps-wasm.js
# LOG_LEVEL = LOG_LEVEL_0
CXX = /usr/bin/g++
CXXFLAGS = -std=c++23 -MMD -MP -Wpedantic -Wall #-O3: optimierungsstufe 3, -MMD .d files, -MP leeres Target, Wextra alle Warnungen
CXXDEBUGFLAGS = -g -fsanitize=address -fno-omit-frame-pointer -O1
CXXDEBUGFLAGS = -g -fno-omit-frame-pointer -O1
# CXXDEBUGFLAGS += -fsanitize=address
WEBCXX = em++
WEBCXXFLAGS = --cache "$(WEBOBJECT_DIR)/emscripten-cache" -Wno-c++11-narrowing -Wno-unused-main -DEMSCRIPTEN
WEBCXXDEBUGFLAGS = -gsource-map -g2 -O0
WEBLDFLAGS = -s EXPORTED_FUNCTIONS="['_free']" -s ASYNCIFY -s EXPORTED_RUNTIME_METHODS="['cwrap', 'getValue']" -lembind -s INITIAL_MEMORY=2GB -s ALLOW_MEMORY_GROWTH=1
WEBCXXFLAGS = --cache "$(WEBOBJECT_DIR)/emscripten-cache" -Wno-c++11-narrowing -Wno-unused-main -DEMSCRIPTEN
# the TEST___FLAGS will be made available with the getBuildOptions() function in webassembly for benchmarking purposes
# TESTCXXFLAGS = -g -ffast-math
# TESTCXXFLAGS = -g -fno-math-errno
TESTCXXFLAGS = -msimd128 # concurrency, TODO: check if does something
# TESTCXXFLAGS += -g
# TESTCXXFLAGS += -s INVOKE_RUN=0
# TESTLDFLAGS = -ffast-math
# TESTLDFLAGS += --closure 1
TESTLDFLAGS = -s ASSERTIONS=1 -s NO_DISABLE_EXCEPTION_CATCHING -g2
WEBCXXDEBUGFLAGS = -gsource-map -g3 -O0 -s ASSERTIONS=2
WEBLDFLAGS = -lembind -s ASYNCIFY -s EXPORTED_RUNTIME_METHODS="['cwrap', 'getValue']" -s INITIAL_MEMORY=2GB -s ALLOW_MEMORY_GROWTH=1
# WEBLDFLAGS += -s STACK_SIZE=131072 -sINVOKE_RUN=0 #-sSTACK_OVERFLOW_CHECK=1 #-s ASSERTIONS=2 # TODO remove
WEBLDFLAGS += -s STACK_SIZE=331072 -sINVOKE_RUN=0 #-sSTACK_OVERFLOW_CHECK=1 #-s ASSERTIONS=2 # TODO remove
WEBLDRELEASEFLAGS = -sDISABLE_EXCEPTION_CATCHING
# treated as error by em++
LDFLAGS = #-lgzutil
IFLAGS = -I../include/
LDFLAGS =
# IFLAGS = -I../include/ -I../../lambda-lanczos/include -I../spectra/include
IFLAGS = -iprefix $$HOME/Projekte/science/cpp-mps/ -iwithprefix include -iwithprefix lambda-lanczos/include -iwithprefix spectra/include
SRCDIRS = $(wildcard */)
IFLAGS += $(foreach dir,$(SRCDIRS), -I$(dir))
IFLAGS += $(foreach dir,$(SRCDIRS), -I../$(dir))
# IFLAGS += $(foreach dir,$(SRCDIRS), -I../$(dir))
IFLAGS += -I../ -I./
CXXFLAGS += $(IFLAGS)
@ -33,39 +48,43 @@ CXXFLAGS += $(IFLAGS)
SRC = $(wildcard *.cpp) $(wildcard */*.cpp)
# OBJECTS = $(SRC:%.cpp=$(OBJECT_DIR)/%.o)
OBJECTS = $($(notdir SRC):%.cpp=$(OBJECT_DIR)/%.o)
WEBOBJECTS = $($(notdir SRC):%.cpp=$(WEBOBJECT_DIR)/%.o)
OBJECT_DIRS = $(foreach dir,$(SRCDIRS), $(OBJECT_DIR)/$(dir))
WEBOBJECT_DIRS = $(foreach dir,$(SRCDIRS), $(WEBOBJECT_DIR)/$(dir))
DEPENDS = ${OBJECTS:.o=.d}
WEBOBJECTS = $($(notdir SRC):%.cpp=$(WEBOBJECT_DIR)/%.o)
WEBOBJECT_DIRS = $(foreach dir,$(SRCDIRS), $(WEBOBJECT_DIR)/$(dir))
WEBDEPENDS = ${WEBOBJECTS:.o=.d}
FMT_MESSAGE="\e[1;34m%s\e[0m %s %s %s\n"
# default: CXXFLAGS += -include pch.hpp
# include the makefiles generated by the -M flag
-include $(DEPENDS)
-include $(WEBDEPENDS)
default:
@printf $(FMT_MESSAGE) "CXXFLAGS = " "$(CXXFLAGS)"
@printf $(FMT_MESSAGE) "LDFLAGS = " "$(LDFLAGS)"
default: CXXFLAGS += -include pch.hpp
default: $(EXEC)
web: CXXFLAGS += $(WEBCXXFLAGS)
web: CXXFLAGS += $(WEBCXXFLAGS) $(TESTCXXFLAGS) $(TESTLDFLAGS) -D BUILD_OPTIONS_STRING="\"$(TESTCXXFLAGS) $(TESTLDFLAGS)\""
# web: LDFLAGS += $(WEBLDFLAGS) $(TESTLDFLAGS)
web: $(WEBEXEC)
webDebug: WEBLDFLAGS += -sASSERTIONS
webDebug: CXXFLAGS += $(WEBCXXDEBUGFLAGS)
webDebug: web
webO0: CXXFLAGS += $(WEBCXXFLAGS) -O0
webO0: TESTCXXFLAGS += -O0
webO0: web
webO1: CXXFLAGS += $(WEBCXXFLAGS) -O1
webO1: TESTCXXFLAGS += -O1
webO1: web
webO2: CXXFLAGS += $(WEBCXXFLAGS) -O2
webO2: TESTCXXFLAGS += -O2
webO2: web
webO3: CXXFLAGS += $(WEBCXXFLAGS) -O3
webO3: TESTCXXFLAGS += -O3
webO3: web
webrun: CXXFLAGS += $(WEBCXXFLAGS)
webrun: $(WEBHTML)
cd $(WEBOBJECT_DIR) && emrun $(WEBHTML)
.PHONY: release
release: CXXFLAGS += -O3
release: default
@ -84,9 +103,7 @@ $(EXEC): $(OBJECT_DIRS) $(OBJECTS)
$(CXX) $(OBJECTS) -o $@ $(CXXFLAGS) $(LDFLAGS)
$(WEBHTML) $(WEBEXEC): $(WEBOBJECT_DIRS) $(WEBOBJECTS)
$(WEBCXX) $(WEBOBJECTS) -o $@ $(CXXFLAGS) $(LDFLAGS) $(WEBLDFLAGS)
# include the makefiles generated by the -M flag
-include $(DEPENDS)
$(WEBCXX) $(WEBOBJECTS) -o $@ $(CXXFLAGS) $(LDFLAGS) $(WEBLDFLAGS)
# rule for all ../build/*.o files
$(OBJECT_DIR)/%.o: $(shell echo $<) %.cpp
@ -106,7 +123,13 @@ $(OBJECT_DIRS) $(WEBOBJECT_DIRS):
# with debug flags
.PHONY += debug run gdb pch clean docs
debug: CXXFLAGS += $(CXXDEBUGFLAGS)
MASSIF_OUT=../massif-mps.out
profile: $(EXEC)
valgrind --tool=massif --stacks=yes --massif-out-file=$(MASSIF_OUT) $(EXEC)
ms_print $(MASSIF_OUT) | less
debug: CXXFLAGS += $(CXXDEBUGFLAGS)
# debug: LDFLAGS += -static-libasan
debug: default
# make with debug flags and run afterwards
@ -126,6 +149,7 @@ pch:
# remove all object and dependecy files
clean:
-rm -r $(OBJECT_DIR)
webclean:
-rm -r $(WEBOBJECT_DIR)
docs:
doxygen .doxygen_config

3
src/algorithm.cpp Normal file
View File

@ -0,0 +1,3 @@
#include "algorithm.hpp"
template bool postUpdateNoop<dtype>(const MPS<dtype>&, unsigned);

20
src/algorithm.hpp Normal file
View File

@ -0,0 +1,20 @@
#pragma once
#include "mps.hpp"
#include "type.hpp"
#include <functional>
/**
* @brief Type for a function to after one algorithm iteration.
* @details
* The return value states wether the results have converged, thus
* if the function returns `true`, the simulation will stop.
*/
template<typename T>
using postUpdateFunction_t = std::function<bool (const MPS<T>&, unsigned)>;
/**
* @brief Default post update function that does nothing
*/
template<typename T>
bool postUpdateNoop(const MPS<T>&, unsigned) { return false; };

288
src/algorithms/dmrg.cpp Normal file
View File

@ -0,0 +1,288 @@
#include "dmrg.hpp"
#include "lanczos.hpp"
#include "type.hpp"
#include "util.hpp"
#include <lambda_lanczos/lambda_lanczos.hpp>
using lambda_lanczos::LambdaLanczos;
#include <print>
#include <vector>
#include <memory>
#include <cmath>
namespace dmrg {
void HeffTheta(
const e::Tensor<dtype, 3>& LP, // vL wL* vL*
const e::Tensor<dtype, 3>& RP, // vR* wR* vR
const e::Tensor<dtype, 4>& W1, // wL wC i i*
const e::Tensor<dtype, 4>& W2, // wC wR j j*
const std::vector<dtype>& thetaIn,
std::vector<dtype>& thetaOut
) {
assert(thetaIn.size() == LP.dimension(0) * W1.dimension(2) * W2.dimension(2) * RP.dimension(2));
assert(thetaOut.size() == LP.dimension(0) * W1.dimension(2) * W2.dimension(2) * RP.dimension(2));
e::array<e::Index, 4> thetaShape = {LP.dimension(0), W1.dimension(2), W2.dimension(2), RP.dimension(2)};
// Reshape theta into a 4D tensor
e::TensorMap<const e::Tensor<dtype, 4>> thetaTensor(thetaIn.data(), thetaShape); // vL i j vR
// e::Tensor<dtype, 4> x = theta_tensor;
auto LPT = LP.contract(thetaTensor, e::array<iP, 1>({iP{2, 0}})); // vL wL* [vL*], [vL] i j vR
auto LPTW1 = LPT.contract(W1, e::array<iP, 2>({iP{1, 0}, iP{2, 3}})); // vL [wL*] [i] j vR, [wL] wC i [i*]
auto LPTW1W2 = LPTW1.contract(W2, e::array<iP, 2>({iP{3, 0}, iP{1, 3}})); // vL [j] vR [wC] i, [wC] wR j [j*]
auto LPTW1W2RP = LPTW1W2.contract(RP, e::array<iP, 2>({iP{1, 0}, iP{3, 1}})); // vL [vR] i [wR] j, [vR*] [wR*] vR
// assigning to map will assign to underlying std::vector
e::TensorMap<e::Tensor<dtype, 4>> thetaOutMap(thetaOut.data(), thetaShape); // vL i j vR
thetaOutMap = LPTW1W2RP;
// std::println("thetaO: {}", thetaOut);
}
template<typename T>
DMRGEngine<T>::DMRGEngine(MPS<T>& psi, const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& H_mpo,
unsigned chiMax, T eps)
: H_mpo(H_mpo), chiMax(chiMax), eps(eps) {
int L = psi.getL();
LPs.resize(L);
RPs.resize(L);
initialize(psi);
}
template<typename T>
void DMRGEngine<T>::initialize(MPS<T>& psi) {
// Initialize LPs and RPs
int D = H_mpo[0]->dimension(0);
int chi = psi.B(0).dimension(0);
e::Tensor<T, 3> LP(chi, D, chi), RP(chi, D, chi);
LP.setZero();
RP.setZero();
LP.chip(0, 1) = identityTensor<T>(chi, chi);
RP.chip(D - 1, 1) = identityTensor<T>(chi, chi);
LPs[0] = LP;
RPs.back() = RP;
for (int i = psi.getL() - 1; i > 1; --i) {
update_RP(i, psi);
}
}
template<typename T>
void DMRGEngine<T>::sweep(MPS<T>& psi) {
// Sweep left to right
for (unsigned i = 0; i < psi.getL() - 2; ++i) {
update_bond(i, psi);
}
// Sweep right to left
for (unsigned i = psi.getL() - 2; i > 0; --i) {
update_bond(i, psi);
}
}
template<typename T>
void DMRGEngine<T>::update_bond(int i, MPS<T>& psi) {
int j = i + 1;
// std::println("update({})", i);
e::Index d1 = H_mpo[i]->dimension(2);
e::Index d2 = H_mpo[j]->dimension(2);
e::Index chi1 = LPs[i].dimension(0);
e::Index chi2 = RPs[j].dimension(2);
e::Index dim = chi1*d1*d2*chi2;
// 1) Find GS (smallest eigenvalue / vector)
// Approaches 1 and 2 are wrong, because we need a more complicated contraction series
// instead of a reshape + simple matrix product for Heff @ theta
// For i) and ii): create effective hamiltonian
// auto LW = LPs[i].contract(*H_mpo[i], e::array<iP, 1>{iP{1, 0}}); // vL [wL] vL*, [wL] wR i i*
// auto LWW = LW.contract(*H_mpo[j], e::array<iP, 1>{iP{2, 1}}); // vL vL* [wR] i i*, [wL] wR j j*
// auto LWWR = LWW.contract(RPs[j], e::array<iP, 1>{iP{4, 1}}); // vL vL* i i* [wR] j j*, vR* [wR*] vR
// e::Tensor<T, 8> Heff = LWWR;// vL vL* i i* j j*, vR* vR
// printTensorDims(Heff, "Heff");
// printTensor(LPs[i], "LP[i]");
// printTensor(RPs[j], "RP[j]");
// printTensor(*H_mpo[i], "H_mpo[i]");
// printTensor(*H_mpo[j], "H_mpo[j]");
// e::Map<e::MatrixX<T>> HeffMat(Heff.data(), dim, dim); // vL vL* i i*, j j* vR* vR
// std::cout << "HeffMat" << std::endl;
// std::cout << HeffMat << std::endl;
// i) Eigen
// Realllyyyy inefficient, it calculates all eigenvalues and we just need the eigenvector to the
// smallest absolute eigenvalue
// e::SelfAdjointEigenSolver<e::MatrixX<T>> es(HeffMat);
// const e::VectorX<T>& S = es.eigenvalues();
// const e::MatrixX<T>& U = es.eigenvectors();
// std::println("Eigenvalues with eigen:");
// std::cout << S << std::endl;
// std::println("Eigenvector with eigen:");
// std::cout << U.rows() << "x" << U.cols() << std::endl;
// e::Map<const e::MatrixX<T>> U1(U.col(0).data(), chi1*d1, d2*chi2);
// std::cout << "theta2 (eigen)" << std::endl;
// std::cout << U1 << std::endl;
// find the smallest eigenvalue
// e::TensorMap<const e::Tensor<T, 4>> theta(U.data(), chi1, d1, d2, chi2);
// auto [Ai, Sj, Bj] = split_truncate_theta<T>(theta, chiMax, eps);
// auto [Ai, Sj, Bj] = split_truncate_theta<T>(U1, chi1, d1, d2, chi2, chiMax, eps);
// ii) Lanczos
// std::println("Eigenvector with Lanczos:");
// auto [eval, theta2Vec] = lanczos(HeffMat, psi.getTheta2(i));
// e::Map<e::MatrixX<T>> theta2(theta2Vec.data(), chi1*d1, d2*chi2);
// std::cout << theta2.rows() << "x" << theta2.cols() << std::endl;
// std::cout << eval << std::endl;
// std::cout << theta2 << std::endl;
// auto [Ai, Sj, Bj] = split_truncate_theta<T>(theta2, chi1, d1, d2, chi2, chiMax, eps);
// auto theta0 = psi.getTheta2(i);
// printTensor(theta0, "theta0");
// e::Map<e::VectorX<T>> theta0v(theta0.data(), theta0.size());
// std::println("Heff @ theta0");
// e::VectorX<T> v = HeffMat * theta0v;
// std::cout << v << std::endl;
// iii) Lanczos with custom contraction function
e::Tensor<T, 4> theta0 = psi.getTheta2(i);
// Use Theta2 as the initial guess for the GS
auto initF = [&theta0](std::vector<T>& v) {
v.resize(theta0.size());
for (unsigned i = 0; i < theta0.size(); i++) {
v.at(i) = theta0.data()[i]; // TODO find nicer solution
}
// v = std::vector<T>(theta0.data(), theta0.size());
};
auto matmul = [this, i, j](const std::vector<T>& thetaIn, std::vector<T>& thetaOut) {
HeffTheta(this->LPs.at(i), this->RPs.at(j), *(this->H_mpo.at(i)), *(this->H_mpo.at(j)), thetaIn, thetaOut);
};
// TODO set eps
LambdaLanczos<T> engine(matmul, dim, false, 1); // Find 1 maximum eigenvalue
engine.init_vector = initF;
std::vector<double> eigenvalues;
std::vector<std::vector<T>> eigenvectors;
engine.run(eigenvalues, eigenvectors);
e::Map<e::MatrixX<T>> theta2(eigenvectors.front().data(), chi1*d1, d2*chi2);
// std::println("Lanczos ev: {}", eigenvalues);
// e::TensorMap<e::Tensor<T, 4>> theta2T(eigenvectors.front().data(), chi1, d1, d2, chi2);
// printTensor(theta2T, "theta2");
auto [Ai, Sj, Bj] = split_truncate_theta<T>(theta2, chi1, d1, d2, chi2, chiMax, eps);
// 2) Truncate
// printTensor(Ai, "Ai");
// printTensor(Sj, "Sj");
// printTensor(Bj, "Bj");
// std::cout << "a: " << a.dimensions() << "\nA: " << Ai.dimensions() << "\nS: " << Sj.dimensions() << "\nB: " << Bj.dimensions() << std::endl;
// renormalize to Bi
e::Tensor<T, 1> SiInv = 1. / psi.S(i);
e::Tensor<T, 2> SiDiv = diagonal(SiInv);
e::Tensor<T, 3> Gi = SiDiv.contract(Ai, e::array<iP, 1>{iP{1, 0}}); // vL [vL], [vL] i vC
// printTensor(SiDiv, "SiDiv");
// printTensor(Gi, "Gi");
// e::Tensor<T, 0> SS = Sj.contract(Sj.conjugate(), e::array<iP, 1>{iP{0,0}});
// std::println("Bond {:03}, Sj contracted with its conjugate = {}", i, SS(0));
// psi.B(i) = Gi.contract(diagonal(Sj), e::array<iP, 1>{iP{2, 0}}); // vL i [vC], [vC] vC
// psi.B(i+1) = std::move(Bj);
// psi.S(i+1) = std::move(Sj);
typename MPS<T>::BTensor_t Bi = Gi.contract(diagonal(Sj), e::array<iP, 1>{iP{2, 0}}); // vL i [vC], [vC] vC
// e::Tensor<bool, 0> equal = (psi.B(i) == Bi).all();
// printTensor(psi.B(i), "old Bi");
// printTensor(Bi, "new Bi");
psi.B(i) = std::move(Bi);
psi.B(j) = std::move(Bj);
psi.S(j) = std::move(Sj);
// 3) Update
update_LP(i, psi);
update_RP(j, psi);
}
template<typename T>
void DMRGEngine<T>::update_RP(int i, MPS<T>& psi) {
// std::println("update_RP({})", i);
int j = i - 1;
const e::Tensor<T, 3>& RP = RPs[i]; // vR* wR* vR
const e::Tensor<T, 3>& B = psi.B(i); // vL i vR
e::Tensor<T, 3> Bc = B.conjugate(); // vL* i* vR*
const e::Tensor<T, 4>& W = *H_mpo[i]; // wL wR i i*
#ifdef DMRGDEBUG
printTensorDims(RP, "RP");
printTensorDims(B, "B");
printTensorDims(Bc, "Bc");
printTensorDims(W, "W");
auto RP2 = B.contract(RP, e::array<iP, 1>{iP{2, 0}}); // vL i [vR], [vR*] wR* vR
e::Tensor<T, 4> RP2T = RP2;
printTensorDims(RP2T, "RP2T");
auto RP3 = RP2T.contract(W, e::array<iP, 2>{iP{1, 3}, iP{2, 1}}); // vL [i] [wR*] vR, wL [wR] i [i*]
e::Tensor<T, 4> RP3T = RP3;
printTensorDims(RP3T, "RP3T");
auto RP4 = RP3.contract(Bc, e::array<iP, 2>{iP{1, 2}, iP{3, 1}}); // vL [vR] wL [i], vL* [i*] [vR*]
e::Tensor<T, 3> RP4T = RP4;
printTensor(RP4T, "RP4T");
#else
auto RP2 = B.contract(RP, e::array<iP, 1>{iP{2, 0}}); // vL i [vR], [vR*] wR* vR
auto RP3 = RP2.contract(W, e::array<iP, 2>{iP{1, 3}, iP{2, 1}}); // vL [i] [wR*] vR, wL [wR] i [i*]
auto RP4 = RP3.contract(Bc, e::array<iP, 2>{iP{1, 2}, iP{3, 1}}); // vL [vR] wL [i], vL* [i*] [vR*]
#endif
RPs[j] = RP4;
}
template<typename T>
void DMRGEngine<T>::update_LP(int i, MPS<T>& psi) {
// std::println("update_LP({})", i);
int j = i + 1;
const e::Tensor<T, 3>& LP = LPs[i]; // vL wL vL*
const e::Tensor<T, 3>& B = psi.B(i); // vL i vR
const e::Tensor<T, 4>& W = *H_mpo[i]; // wL wR i i*
e::Tensor<T, 1> SjInvV = 1. / psi.S(j);
e::Tensor<T, 2> SjInv = diagonal(SjInvV);
auto G = B.contract(SjInv, e::array<iP, 1>{iP{2, 0}}); // vL i [vR], [vR*] vR
e::Tensor<T, 2> Si = diagonal(psi.S(i));
auto A = Si.contract(G, e::array<iP, 1>{iP{1, 0}}); // vL [vL*], [vL] i vR
auto Ac = A.conjugate(); // vL* i* vR*
#ifdef DMRGDEBUG
auto LP2 = LP.contract(A, e::array<iP, 1>{iP{2, 0}}); // vL wL* [vL*], [vL] i vR
e::Tensor<T, 4> LP2T = LP2;
printTensorDims(LP2T, "LP2T");
printTensorDims(W, "W");
auto LP3 = W.contract(LP2T, e::array<iP, 2>{iP{0, 1}, iP{3, 2}}); // [wL] wR i [i*], vL [wL*] [i] vR
e::Tensor<T, 4> LP3T = LP3;
printTensorDims(LP3T, "LP3T");
e::Tensor<T, 3> AcT = Ac;
printTensorDims(AcT, "AcT");
auto LP4 = AcT.contract(LP3, e::array<iP, 2>{iP{0, 2}, iP{1, 1}}); // [vL*] [i*] vR*, wR [i] [vL] vR
e::Tensor<T, 3> LP4T = LP4;
printTensor(LP4T, "LP4T");
LPs[j] = LP4T;
#else
auto LP2 = LP.contract(A, e::array<iP, 1>{iP{2, 0}}); // vL wL* [vL*], [vL] i vR
auto LP3 = W.contract(LP2, e::array<iP, 2>{iP{0, 1}, iP{3, 2}}); // [wL] wR i [i*], vL [wL*] [i] vR
auto LP4 = Ac.contract(LP3, e::array<iP, 2>{iP{0, 2}, iP{1, 1}}); // [vL*] [i*] vR*, wR [i] [vL] vR
LPs[j] = LP4;
#endif
}
template<typename T>
void run(MPS<T>& psi, const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& H_mpos, unsigned chiMax, T eps, unsigned nSweeps, postUpdateFunction_t<T> postUpdateFunc) {
auto engine = DMRGEngine(psi, H_mpos, chiMax, eps);
for (unsigned i = 0; i < nSweeps; i++) {
engine.sweep(psi);
if (postUpdateFunc(psi, i)) break;
}
}
template class dmrg::DMRGEngine<dtype>;
template void run<dtype>(MPS<dtype>& psi, const std::vector<std::shared_ptr<e::Tensor<dtype, 4>>>& H_mpos, unsigned chiMax, dtype eps, unsigned nSweeps, postUpdateFunction_t<dtype> postUpdateFunc);
} // namespace dmrg

52
src/algorithms/dmrg.hpp Normal file
View File

@ -0,0 +1,52 @@
#pragma once
#include "algorithm.hpp"
#include "mps.hpp"
#include <eigen3/Eigen/Eigen>
/**
* @brief Returns a \f$L\times L \f$ identity matrix in tensor form, where \f$L = \text{min}(N, M)\f$
*/
template <typename T>
e::Tensor<T, 2> identityTensor(e::Index N, e::Index M) {
e::Tensor<dtype, 2> t(N, M);
t.setZero();
e::Index min = std::min(N, M);
for (e::Index i = 0; i < min; i++) {
t(i, i) = 1;
}
return t;
}
/**
* @brief Implementation of the Density Matrix Renormalization Group (DMRG) algorithm
*/
namespace dmrg {
/**
* @brief Helper class that keeps the left and right blocks of the MPS
* @note Always call with the same mps (psi) that it was initialized with!
* It would have been better if the engine would take ownership of the MPS during sweeps,
* but this would make calling an update function in between more difficult.
*/
template<typename T>
class DMRGEngine {
public:
std::vector<std::shared_ptr<e::Tensor<T, 4>>> H_mpo;
std::vector<e::Tensor<T, 3>> LPs, RPs;
int chiMax;
T eps;
DMRGEngine(MPS<T>& psi, const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& H_mpo,
unsigned chiMax = 100, T eps = 1e-12);
void initialize(MPS<T>& psi);
void sweep(MPS<T>& psi);
void update_bond(int i, MPS<T>& psi);
void update_RP(int i, MPS<T>& psi);
void update_LP(int i, MPS<T>& psi);
};
template<typename T>
void run(MPS<T>& psi, const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& H_mpos, unsigned chiMax, T eps, unsigned nSweeps, postUpdateFunction_t<T> postUpdateFunc);
}

201
src/algorithms/tebd.cpp
View File

@ -1,10 +1,11 @@
#include "tebd.hpp"
#include "type.hpp"
#include "util.hpp"
#include "model.hpp"
#include <cstdio>
#include <ostream>
#include <unordered_map>
#include <print>
#include <eigen3/unsupported/Eigen/MatrixFunctions>
template<typename InType, typename OutType>
std::vector<std::shared_ptr<OutType>> forEachUnique(std::function<OutType(const InType&)> f, const std::vector<std::shared_ptr<InType>>& inRange) {
@ -25,107 +26,109 @@ std::vector<std::shared_ptr<OutType>> forEachUnique(std::function<OutType(const
}
template<typename T>
e::Tensor<T, 4> getTimeEvolutionMatrix(const e::MatrixX<T>& bond, T dt, unsigned localDim) {
e::MatrixX<T> U = (-bond * dt).exp();
// std::cout << "U = \n" << U << "\n\n";
return bondMatrix2Tensor<T>(U, localDim);
}
template<typename T>
Bonds<T> getTimeEvolutionMatrices(const std::vector<std::shared_ptr<e::MatrixX<T>>>& bonds, T dt, unsigned int localDim) {
std::function<e::Tensor<T, 4>(const e::MatrixX<T>&)> f = std::bind(getTimeEvolutionMatrix<T>, std::placeholders::_1, dt, localDim);
std::vector<std::shared_ptr<e::Tensor<T, 4>>> U_Bonds = forEachUnique(f, bonds);
return U_Bonds;
}
template<typename T>
void updateBond(MPS<T>& psi, size_t i, const e::Tensor<T, 4>& U, unsigned chiMax, T eps) {
// std::println("BOND {}", i);
auto theta = psi.getTheta2(i); // vL i j vR
// apply U: i j i* j*
// contraction gives vL vR i j -> shuffle to vL i j vR
// std::cout << theta.dimension(2) << " - " << U.dimension(2) << " / " << theta.dimensions() << std::endl;
e::Tensor<T, 4> a = theta.contract(U, e::array<iP, 2>({iP{1, 2}, iP{2, 3}}));
e::Tensor<T, 4> b = a.shuffle(e::array<int, 4>{0, 2, 3, 1});
// e::Tensor<T, 4> b = a.shuffle(e::array<int, 4>{2, 0, 1, 3});
// printTensor(b, "Utheta");
// SVD and truncation
auto [Ai, Sj, Bj] = split_truncate_theta<T>(b, chiMax, eps);
// printTensor(Ai, "Ai");
// printTensor(Sj, "Sj");
// printTensor(Bj, "Bj");
// std::cout << "a: " << a.dimensions() << "\nA: " << Ai.dimensions() << "\nS: " << Sj.dimensions() << "\nB: " << Bj.dimensions() << std::endl;
// renormalize to Bi
e::Tensor<T, 1> SiInv = 1. / psi.S(i);
e::Tensor<T, 2> SiDiv = diagonal(SiInv);
e::Tensor<T, 3> Gi = SiDiv.contract(Ai, e::array<iP, 1>{iP{1, 0}}); // vL [vL], [vL] i vC
// printTensor(SiDiv, "SiDiv");
// printTensor(Gi, "Gi");
// e::Tensor<T, 0> SS = Sj.contract(Sj.conjugate(), e::array<iP, 1>{iP{0,0}});
// std::println("Bond {:03}, Sj contracted with its conjugate = {}", i, SS(0));
// psi.B(i) = Gi.contract(diagonal(Sj), e::array<iP, 1>{iP{2, 0}}); // vL i [vC], [vC] vC
// psi.B(i+1) = std::move(Bj);
// psi.S(i+1) = std::move(Sj);
MPS<dtype>::BTensor_t Bi = Gi.contract(diagonal(Sj), e::array<iP, 1>{iP{2, 0}}); // vL i [vC], [vC] vC
// e::Tensor<bool, 0> equal = (psi.B(i) == Bi).all();
// printTensor(psi.B(i), "old Bi");
// printTensor(Bi, "new Bi");
psi.B(i) = Bi;
psi.B(i+1) = Bj;
psi.S(i+1) = Sj;
}
template<typename T>
void runTEBD_brickwall(MPS<T>& psi, const Bonds<T>& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc) {
for (unsigned n = 0; n < nSteps; n++) {
// std::print("i = {:03}\r", n);
// std::flush(std::cout);
for (unsigned row = 0; row < 2; row++) {
for (unsigned i = row; i < bonds.size(); i+=2) {
// std::cout << "i = " << i << std::endl;
updateBond(psi, i, *(bonds.at(i)), chiMax, eps);
}
}
postUpdateFunc(psi, n);
namespace tebd {
template<typename T>
e::Tensor<T, 4> getTimeEvolutionMatrix(const e::MatrixX<T>& bond, T dt, unsigned localDim) {
e::MatrixX<T> U = (-bond * dt).exp();
// std::cout << "U = \n" << U << "\n\n";
return bondMatrix2Tensor<T>(U, localDim);
}
}
// TODO
/**
* @brief Run TEBD using a second-order decomposition
* @todo I think this is currently only correct for even system sizes => check decomposition
*/
template<typename T>
void runTEBD_secondOrder(MPS<T>& psi, const Bonds<T>& bonds, const Bonds<T>& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc) {
for (unsigned n = 0; n < nSteps; n++) {
// std::print("i = {:03}\r", n);
// std::flush(std::cout);
for (unsigned row = 0; row < 2; row++) {
const Bonds<T>* useBonds = &bonds;
if (row == 0 && (n != 0 || n != nSteps-1)) {
useBonds = &bondsHalfDT;
}
for (unsigned i = row; i < bonds.size(); i+=2) {
// std::cout << "i = " << i << std::endl;
updateBond(psi, i, *(useBonds->at(i)), chiMax, eps);
}
}
postUpdateFunc(psi, n);
template<typename T>
Bonds<T> getTimeEvolutionMatrices(const std::vector<std::shared_ptr<e::MatrixX<T>>>& bonds, T dt, unsigned int localDim) {
std::function<e::Tensor<T, 4>(const e::MatrixX<T>&)> f = std::bind(getTimeEvolutionMatrix<T>, std::placeholders::_1, dt, localDim);
std::vector<std::shared_ptr<e::Tensor<T, 4>>> U_Bonds = forEachUnique(f, bonds);
return U_Bonds;
}
}
// Instantiations
template<typename T>
void updateBond(MPS<T>& psi, size_t i, const e::Tensor<T, 4>& U, unsigned chiMax, T eps) {
// std::println("BOND {}", i);
auto theta = psi.getTheta2(i); // vL i j vR
// apply U: i j i* j*
// contraction gives vL vR i j -> shuffle to vL i j vR
// std::cout << theta.dimension(2) << " - " << U.dimension(2) << " / " << theta.dimensions() << std::endl;
e::Tensor<T, 4> a = theta.contract(U, e::array<iP, 2>({iP{1, 2}, iP{2, 3}}));
e::Tensor<T, 4> b = a.shuffle(e::array<int, 4>{0, 2, 3, 1});
// e::Tensor<T, 4> b = a.shuffle(e::array<int, 4>{2, 0, 1, 3});
// printTensor(b, "Utheta");
// SVD and truncation
auto [Ai, Sj, Bj] = split_truncate_theta<T>(b, chiMax, eps);
// printTensor(Ai, "Ai");
// printTensor(Sj, "Sj");
// printTensor(Bj, "Bj");
// std::cout << "a: " << a.dimensions() << "\nA: " << Ai.dimensions() << "\nS: " << Sj.dimensions() << "\nB: " << Bj.dimensions() << std::endl;
// renormalize to Bi
e::Tensor<T, 1> SiInv = 1. / psi.S(i);
e::Tensor<T, 2> SiDiv = diagonal(SiInv);
e::Tensor<T, 3> Gi = SiDiv.contract(Ai, e::array<iP, 1>{iP{1, 0}}); // vL [vL], [vL] i vC
// printTensor(SiDiv, "SiDiv");
// printTensor(Gi, "Gi");
// e::Tensor<T, 0> SS = Sj.contract(Sj.conjugate(), e::array<iP, 1>{iP{0,0}});
// std::println("Bond {:03}, Sj contracted with its conjugate = {}", i, SS(0));
// psi.B(i) = Gi.contract(diagonal(Sj), e::array<iP, 1>{iP{2, 0}}); // vL i [vC], [vC] vC
// psi.B(i+1) = std::move(Bj);
// psi.S(i+1) = std::move(Sj);
MPS<dtype>::BTensor_t Bi = Gi.contract(diagonal(Sj), e::array<iP, 1>{iP{2, 0}}); // vL i [vC], [vC] vC
// e::Tensor<bool, 0> equal = (psi.B(i) == Bi).all();
// printTensor(psi.B(i), "old Bi");
// printTensor(Bi, "new Bi");
psi.B(i) = Bi;
psi.B(i+1) = Bj;
psi.S(i+1) = Sj;
}
template<typename T>
void runBrickwall(MPS<T>& psi, const Bonds<T>& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc) {
for (unsigned n = 0; n < nSteps; n++) {
// std::print("i = {:03}\r", n);
// std::flush(std::cout);
for (unsigned row = 0; row < 2; row++) {
for (unsigned i = row; i < bonds.size(); i+=2) {
// std::cout << "i = " << i << std::endl;
// std::println("Updating bond {:02}", i);
updateBond(psi, i, *(bonds.at(i)), chiMax, eps);
}
}
if (postUpdateFunc(psi, n)) break;
}
}
// TODO
/**
* @brief Run TEBD using a second-order decomposition
* @todo I think this is currently only correct for even system sizes => check decomposition
*/
template<typename T>
void runSecondOrder(MPS<T>& psi, const Bonds<T>& bonds, const Bonds<T>& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc) {
for (unsigned n = 0; n < nSteps; n++) {
// std::print("i = {:03}\r", n);
// std::flush(std::cout);
for (unsigned row = 0; row < 2; row++) {
const Bonds<T>* useBonds = &bonds;
if (row == 0 && (n != 0 || n != nSteps-1)) {
useBonds = &bondsHalfDT;
}
for (unsigned i = row; i < bonds.size(); i+=2) {
// std::cout << "i = " << i << std::endl;
updateBond(psi, i, *(useBonds->at(i)), chiMax, eps);
}
}
if (postUpdateFunc(psi, n)) break;
}
}
// Instantiations
template e::Tensor<dtype, 4> getTimeEvolutionMatrix<dtype>(const e::MatrixX<dtype>& bond, dtype dt, unsigned localDim);
template Bonds<dtype> getTimeEvolutionMatrices<dtype>(const std::vector<std::shared_ptr<e::MatrixX<dtype>>>& bonds, dtype dt, unsigned int localDim);
template void updateBond<dtype>(MPS<dtype>& psi, size_t i, const e::Tensor<dtype, 4>& U, unsigned chiMax, dtype eps);
template void runBrickwall<dtype>(MPS<dtype>& psi, const Bonds<dtype>& bonds, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t<dtype> postUpdateFunc);
template void runSecondOrder<dtype>(MPS<dtype>& psi, const Bonds<dtype>& bonds, const Bonds<dtype>&, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t<dtype> postUpdateFunc);
}; // namespace tebd
template std::vector<std::shared_ptr<e::Tensor<dtype, 4>>> forEachUnique(std::function<e::Tensor<dtype, 4> (const e::MatrixX<dtype> &)> f, const std::vector<std::shared_ptr<e::MatrixX<dtype>>>&);
template e::Tensor<dtype, 4> getTimeEvolutionMatrix<dtype>(const e::MatrixX<dtype>& bond, dtype dt, unsigned localDim);
template Bonds<dtype> getTimeEvolutionMatrices<dtype>(const std::vector<std::shared_ptr<e::MatrixX<dtype>>>& bonds, dtype dt, unsigned int localDim);
template void updateBond<dtype>(MPS<dtype>& psi, size_t i, const e::Tensor<dtype, 4>& U, unsigned chiMax, dtype eps);
template void postUpdateNoop<dtype>(const MPS<dtype>&, unsigned);
template void runTEBD_brickwall<dtype>(MPS<dtype>& psi, const Bonds<dtype>& bonds, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t<dtype> postUpdateFunc);
template void runTEBD_secondOrder<dtype>(MPS<dtype>& psi, const Bonds<dtype>& bonds, const Bonds<dtype>&, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t<dtype> postUpdateFunc);

65
src/algorithms/tebd.hpp
View File

@ -1,9 +1,8 @@
#pragma once
#include "model.hpp"
#include "algorithm.hpp"
#include "mps.hpp"
#include <eigen3/Eigen/Eigen>
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
#include <eigen3/unsupported/Eigen/MatrixFunctions>
/**
* @brief Apply the function f on every element if inRange, but reuse a cached result when two pointers point to the same object
@ -15,40 +14,40 @@ template<typename T>
using Bonds = std::vector<std::shared_ptr<e::Tensor<T, 4>>>;
/**
* @brief Apply time evolution on a bond
* @details
* \f[U = \exp{ B \, \del t} \f]
* @brief Implementation of the Time Evolved block Decimation (TEBD) algorithm
*/
template<typename T>
e::Tensor<T, 4> getTimeEvolutionMatrix(const e::MatrixX<T>& bond, T dt, unsigned localDim);
namespace tebd {
/**
* @brief Apply time evolution on a bond
* @details
* \f[U = \exp{ B \, \partial t} \f]
*/
template<typename T>
e::Tensor<T, 4> getTimeEvolutionMatrix(const e::MatrixX<T>& bond, T dt, unsigned localDim);
template<typename T>
Bonds<T> getTimeEvolutionMatrices(const std::vector<std::shared_ptr<e::MatrixX<T>>>& bonds, T dt, unsigned int localDim);
template<typename T>
Bonds<T> getTimeEvolutionMatrices(const std::vector<std::shared_ptr<e::MatrixX<T>>>& bonds, T dt, unsigned int localDim);
/**
* @brief Apply U on bond i, perform SVD and truncate, update MPS
*/
template<typename T>
void updateBond(MPS<T>& psi, size_t i, const e::Tensor<T, 4>& U, unsigned chiMax, T eps);
/**
* @brief Apply U on bond i, perform SVD and truncate, update MPS
*/
template<typename T>
void updateBond(MPS<T>& psi, size_t i, const e::Tensor<T, 4>& U, unsigned chiMax, T eps);
template<typename T>
using postUpdateFunction_t = std::function<void (const MPS<T>&, unsigned)>;
template<typename T>
void postUpdateNoop(const MPS<T>&, unsigned) {};
/**
* @brief Brickwall TEBS scheme
* @details
* Apply U on every second bond starting from the first,
* then apply U on every second bond starting from the second
*
* First order scheme
*
* @param postUpdateFunc function that can be used to track expectation values after each update
*/
template<typename T>
void runBrickwall(MPS<T>& psi, const Bonds<T>& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc);
/**
* @brief Brickwall TEBS scheme
* @details
* Apply U on every second bond starting from the first,
* then apply U on every second bond starting from the second
*
* First order scheme
*
* @param postUpdateFunc function that can be used to track expectation values after each update
*/
template<typename T>
void runTEBD_brickwall(MPS<T>& psi, const Bonds<T>& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc);
template<typename T>
void runTEBD_secondOrder(MPS<T>& psi, const Bonds<T>& bonds, const Bonds<T>& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc);
template<typename T>
void runSecondOrder(MPS<T>& psi, const Bonds<T>& bonds, const Bonds<T>& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t<T> postUpdateFunc);
}; // namespace tebd

5
src/benchmark.cpp
View File

@ -41,12 +41,13 @@ Results _runTEBD(unsigned L, double J, double g, unsigned nSteps, dtype dt, dtyp
result.M.push_back(psi.getSiteExpectationValueSum(sigma_z).real());
result.E.push_back(model.energy(psi).real());
result.S.push_back(psi.getEntanglementEntropy(middleBond).real());
return false;
});
// initial values now
F(psi, 0);
const Bonds<dtype> bonds = getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim);
runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, F);
const Bonds<dtype> bonds = tebd::getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim);
tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, F);
return result;
};

122
src/main.cpp
View File

@ -7,11 +7,29 @@
#include "tebd.hpp"
#include "util.hpp"
#include "type.hpp"
#include "dmrg.hpp"
#include "wasm/run.hpp"
#include "benchmark.hpp"
using namespace std::complex_literals;
void testWasmCode() {
unsigned L = 6;
unsigned nMax = 3;
auto model = std::make_shared<BoseHubbardModel1D<dtype>>(L, dtype(1.), dtype(0.05), dtype(0.5), nMax);
auto psi = initMPS_nParticles<dtype>(L, nMax, 1);
double eps = 1e-10;
unsigned chiMax = 20;
double dtReal = 0.00;
double dtImag = 0.05;
unsigned nSteps = 300;
ems::val callback;
BoseHubbardResults result(callback, 0);
runTEBD(psi, model, nSteps, dtReal, dtImag, eps, chiMax, TEBD_TYPE::BRICKWALL, result);
}
void testSVD() {
e::Tensor<dtype, 4> t4(2, 3, 4, 5);
auto data = t4.data();
@ -30,19 +48,20 @@ void testSVD() {
void testBoseHubbardTEBD() {
unsigned L = 8;
unsigned nMax = 5;
auto model = BoseHubbardModel1D<dtype>(L, dtype(1.), dtype(2.50), dtype(0.5), nMax);
unsigned L = 6;
unsigned nMax = 3;
auto model = BoseHubbardModel1D<dtype>(L, dtype(1.), dtype(0.05), dtype(0.5), nMax);
dtype eps = 1e-10;
unsigned chiMax = 20;
dtype dt = 0.05;
dtype dt = 0.05i;
unsigned nSteps = 300;
std::println("dt={}", dt);
// std::cout << "leftT:\n" << *(model.getH_Bonds().front()) << std::endl;
// std::cout << "leftT2:\n" << bondMatrix2Tensor(*(model.getH_BondsMatrices().front()), model.localDim) << std::endl;
std::cout << "middleT:\n" << *(model.getH_Bonds().at(1)) << std::endl;
//
auto bonds = getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim);
auto bonds = tebd::getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim);
std::cout << "N:\n" << model.getLadderOperators().getN() << std::endl;
std::cout << "a:\n" << model.getLadderOperators().getA() << std::endl;
std::cout << "aDagger:\n" << model.getLadderOperators().getADagger() << std::endl;
@ -50,36 +69,36 @@ void testBoseHubbardTEBD() {
// std::cout << "H_Bond:\n" << *(bonds.front()) << std::endl;
// MPS<dtype> psi = initMPS_spinUp<dtype>(L);
MPS<dtype> psi = initMPS_nParticles<dtype>(L, nMax, 2);
MPS<dtype> psi = initMPS_nParticles<dtype>(L, nMax, 1);
unsigned nSteps = 300;
// postUpdateFunction_t<dtype> fUpdate = postUpdateNoop<dtype>;
e::TensorMap<const e::Tensor<dtype, 2>> particleNumberOperator(model.getLadderOperators().getN().data(), model.localDim, model.localDim);
postUpdateFunction_t<dtype> fUpdate = [&model, &particleNumberOperator](const MPS<dtype>& psi, unsigned i) {
std::println("{:03}: N={:.2f} {}", i, psi.getSiteExpectationValueSum(particleNumberOperator).real() / model.L, psi.getChi());
return false;
};
// auto fUpdate = postUpdateFunction_t<dtype>([](const MPS<dtype>& psi, unsigned i){
// std::println("i={:03}: <psi|psi>={}", i, psi.collapse(0));
// });
std::println("E = {}", model.energy(psi).real());
// std::println("E = {}", model.energy(psi).real());
// std::println("m = {}", psi.getSiteExpectationValues(sigma_z));
std::println("chiS: {:b}", psi.getChi());
psi.print();
std::println("<psi|psi>={}", psi.collapse(0));
runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, fUpdate);
std::println("E = {}", model.energy(psi).real());
std::println("<N> = {}", psi.getSiteExpectationValueSum(particleNumberOperator).real() / model.L);
// std::println("m = {}", psi.getSiteExpectationValues(sigma_z));
std::println("chiS: {:b}", psi.getChi());
// std::println("chiS: {:b}", psi.getChi());
// psi.print();
std::println("<psi|psi>={}", psi.collapse(0));
// std::println("<psi|psi>={}", psi.collapse(0));
tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, fUpdate);
// std::println("E = {}", model.energy(psi).real());
// std::println("<N> = {}", psi.getSiteExpectationValueSum(particleNumberOperator).real() / model.L);
// std::println("m = {}", psi.getSiteExpectationValues(sigma_z));
// std::println("chiS: {:b}", psi.getChi());
// psi.print();
// std::println("<psi|psi>={}", psi.collapse(0));
}
void testTFI_TEBD() {
unsigned L = 8;
unsigned L = 40;
std::println("Make model");
auto model = TransverseFieldIsingModel1D<dtype>(L, dtype(1.), dtype(1.5));
// auto& mpo = *(model.getH_MPO().front());
@ -109,21 +128,21 @@ void testTFI_TEBD() {
e::Tensor<dtype, 2> sigma_z(2, 2);
sigma_z.setValues({{1, 0}, {0, -1}});
dtype eps = 1e-10;
unsigned chiMax = 20;
dtype dt = 0.05i;
dtype eps = 1e-20;
unsigned chiMax = 30;
dtype dt = 0.05;
std::println("dt={}", dt);
std::cout << "leftT:\n" << *(model.getH_Bonds().front()) << std::endl;
std::cout << "leftT2:\n" << bondMatrix2Tensor(*(model.getH_BondsMatrices().front()), model.localDim) << std::endl;
std::cout << "middleT:\n" << *(model.getH_Bonds().at(1)) << std::endl;
auto bonds = getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim);
auto bonds = tebd::getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim);
std::cout << "U_Bond:\n" << *(bonds.front()) << std::endl;
// std::cout << "H_Bond:\n" << *(bonds.front()) << std::endl;
// MPS<dtype> psi = initMPS_spinUp<dtype>(L);
MPS<dtype> psi = initMPS_spinRight<dtype>(L);
unsigned nSteps = 100;
unsigned nSteps = 200;
postUpdateFunction_t<dtype> fUpdate = postUpdateNoop<dtype>;
// auto fUpdate = postUpdateFunction_t<dtype>([](const MPS<dtype>& psi, unsigned i){
// std::println("i={:03}: <psi|psi>={}", i, psi.collapse(0));
@ -134,7 +153,7 @@ void testTFI_TEBD() {
psi.print();
std::println("<psi|psi>={}", psi.collapse(0));
runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, fUpdate);
tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, fUpdate);
// updateBond(psi, 0, *(bonds.at(0)), chiMax, eps);
// updateBond(psi, 2, *(bonds.at(2)), chiMax, eps);
// updateBond(psi, 4, *(bonds.at(4)), chiMax, eps);
@ -148,7 +167,7 @@ void testTFI_TEBD() {
std::println("E = {}", model.energy(psi).real());
std::println("m = {}", psi.getSiteExpectationValues(sigma_z));
std::println("chiS: {:b}", psi.getChi());
psi.print();
// psi.print();
std::println("<psi|psi>={}", psi.collapse(0));
// for (unsigned i = 0; i < L; i++) {
@ -157,10 +176,59 @@ void testTFI_TEBD() {
}
void testTFI_DMRG() {
unsigned L = 8;
std::println("Make model");
auto model = TransverseFieldIsingModel1D<dtype>(L, dtype(1.), dtype(1));
e::Tensor<dtype, 2> sigma_z(2, 2);
sigma_z.setValues({{1, 0}, {0, -1}});
std::cout << *(model.getH_BondsMatrices().at(0)) << std::endl;
printTensor(*model.getH_Bonds().at(0), "H0");
std::cout << *(model.getH_BondsMatrices().at(1)) << std::endl;
printTensor(*model.getH_Bonds().at(1), "H1");
std::cout << *(model.getH_BondsMatrices().at(6)) << std::endl;
printTensor(*model.getH_Bonds().at(6), "H6");
dtype eps = 1e-20;
unsigned chiMax = 10;
MPS<dtype> psi = initMPS_spinUp<dtype>(L);
for (unsigned i = 0; i < 7; i++) {
std::println("{}", psi.getBondExpectationValue(*(model.getH_Bonds()[i]), i));
}
// MPS<dtype> psi = initMPS_spinRight<dtype>(L);
// MPS<dtype> psi = initMPS_spinDown<dtype>(L);
dtype E = model.energy(psi);
std::println("E{} = {}", -1, E.real());
std::println("M{} = {}", -1, psi.getSiteExpectationValueSum(sigma_z));
printTensor(*(model.getH_MPO().at(0)), "H_mpo[0]");
printTensor(*(model.getH_MPO().at(1)), "H_mpo[1]");
auto engine = dmrg::DMRGEngine(psi, model.getH_MPO(), chiMax, eps);
std::println("Initialized");
E = model.energy(psi);
std::println("E{} = {}", 0, E);
// engine.update_bond(0);
for (unsigned i = 0; i < 4; i++) {
engine.sweep(psi);
std::println("\nAfter sweep #{}", i+1);
E = model.energy(psi);
std::println("E = {:.2f}", E);
std::println("M = {:.2f}", psi.getSiteExpectationValueSum(sigma_z));
}
// engine.psi.print();
}
int main() {
testTFI_DMRG();
// testSVD(); return 0;
testBoseHubbardTEBD();
// testBoseHubbardTEBD();
// testTFI_TEBD();
// testWasmCode();
// return runBenchmark_oneAtATime();
// return runBenchmark_all();

2
src/model.cpp
View File

@ -12,6 +12,7 @@ T Model1D<T>::energy(const MPS<T>& psi) const {
return E;
}
template<typename T>
e::Tensor<T, 4> constructMPO(const std::initializer_list<std::initializer_list<e::MatrixX<T>>>& matrices) {
size_t rows = matrices.size();
@ -42,7 +43,6 @@ e::Tensor<T, 4> constructMPO(const std::initializer_list<std::initializer_list<e
}
template<typename T>
e::MatrixX<T> MPO_Tensor2Matrix(const e::Tensor<T, 4>& mpo) {
size_t dim = mpo.dimension(0);

5
src/model.hpp
View File

@ -17,12 +17,13 @@ using namespace std::complex_literals;
template<typename T>
e::Tensor<T, 4> constructMPO(const std::initializer_list<std::initializer_list<e::MatrixX<T>>>& matrices);
template<typename T>
e::MatrixX<T> MPO_Tensor2Matrix(const e::Tensor<T, 4>& mpo);
/**
* @brief Reshape a N*localDim x N*localDim matrix into a N x localDim x N x localDim tensor
* @brief Reshape a \f$N \, d \times N \, d\f$ matrix into a \f$N \times d \times N \times d\f] tensor
*/
template<typename T>
e::Tensor<T, 4> bondMatrix2Tensor(const e::MatrixX<T>& bondMatrix, unsigned localDim);
@ -39,7 +40,7 @@ class Model1D {
virtual ~Model1D() {};
inline const std::vector<std::shared_ptr<e::MatrixX<T>>>& getH_BondsMatrices() const { return H_BondsMatrices; }
inline const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& getH_Bonds() const { return H_Bonds; }
inline const std::vector<std::shared_ptr<e::Tensor<T, 4>>>& getH_MPO() const { return H_MPO; }
inline std::vector<std::shared_ptr<e::Tensor<T, 4>>>& getH_MPO() { return H_MPO; }
const size_t L;
const unsigned localDim;
T energy(const MPS<T>& psi) const;

44
src/models/bose_hubbard.cpp
View File

@ -1,6 +1,7 @@
#include "bose_hubbard.hpp"
#include <eigen3/unsupported/Eigen/KroneckerProduct>
template<typename T>
LadderOperators<T>::LadderOperators(unsigned dim) : identity(dim, dim), aDagger(dim, dim), a(dim, dim), n(dim, dim) {
identity.setZero();
@ -22,34 +23,37 @@ LadderOperators<T>::LadderOperators(unsigned dim) : identity(dim, dim), aDagger(
template<typename T>
void BoseHubbardModel1D<T>::initH_Bonds() {
T tHalf = this->t * static_cast<T>(0.5);
T tFull = this->t;
T UHalf = this->U * static_cast<T>(0.5);
T UQuarter = this->U * static_cast<T>(0.25);
T muHalf = this->mu * static_cast<T>(0.5);
/// @todo TODO figure out which is right. Doesnt t need to be halfed at the edges and not the middle?
e::MatrixX<T> N_NminusI = this->lo.getN() * (this->lo.getN() - this->lo.getIdentity());
e::MatrixX<T> singleParticlePart(
(
e::kroneckerProduct(this->lo.getIdentity(), N_NminusI)
- e::kroneckerProduct(N_NminusI, this->lo.getIdentity())
) * (this->U/static_cast<T>(2))
- e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (this->mu/static_cast<T>(2))
- e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (this->mu/static_cast<T>(2))
e::MatrixX<T> twoParticlePart =
- this->t * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA())
- this->t * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger()
);
// TODO this is wrong, needs A and Adagger; check which t or t/half or identities are required
auto left = std::make_shared<e::MatrixX<T>>(
singleParticlePart
- tFull * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA())
- tFull * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger())
twoParticlePart
+ e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) * UHalf
+ e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) * UQuarter
- e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (this->mu)
- e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (muHalf)
);
auto right = std::make_shared<e::MatrixX<T>>(
singleParticlePart
- tFull * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA())
- tFull * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger())
twoParticlePart
+ e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) * UQuarter
+ e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) * UHalf
- e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (muHalf)
- e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (this->mu)
);
auto middle = std::make_shared<e::MatrixX<T>>(
singleParticlePart
- tFull * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA())
- tFull * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger())
twoParticlePart
+ e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) * UQuarter
+ e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) * UQuarter
- e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (muHalf)
- e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (muHalf)
);
// The matrices are N*N x N*N -> reshape to N x N x N x N tensors
unsigned N = this->localDim;
@ -90,7 +94,7 @@ void BoseHubbardModel1D<T>::initH_MPO() {
zero.setZero();
auto W = std::make_shared<e::Tensor<T, 4>>(std::move(constructMPO<T>({
{this->lo.getIdentity(), this->lo.getADagger(), this->lo.getA(), (this->U / 2. * (this->lo.getN() * this->lo.getN() - this->lo.getN()) - this->mu * this->lo.getN())},
{this->lo.getIdentity(), this->lo.getADagger(), this->lo.getA(), this->U / 2. * (this->lo.getN() * this->lo.getN() - this->lo.getN()) - this->mu * this->lo.getN()},
{zero, zero, zero, (-this->t * this->lo.getA())},
{zero, zero, zero, (-this->t * this->lo.getADagger())},
{zero, zero, zero, this->lo.getIdentity()}

38
src/models/bose_hubbard.hpp
View File

@ -1,6 +1,7 @@
#pragma once
#include "model.hpp"
/**
* @brief Helper class constructing identity, particle number and both ladder operators in one loop.
*/
@ -19,14 +20,44 @@ class LadderOperators {
e::MatrixX<T> n;
};
/**
* @brief Bose-Hubbard Model in 1 dimension
* Latex test \f$ 1 +2 \otimes 3\f$
* @details
* The Hamiltonian:
* \f[
* H =
* -t \sum_{i}\left(
* \hat{a}^\dagger_i \hat{a}_{i+1}
* + \hat{a}_i \hat{a}^\dagger_{i+1}
* \right)
* + \frac{U}{2} \sum_i \hat{n}_i \left(\hat{n}_i-\mathbbm{1}\right)
* - \mu \sum_i \hat{n}_i
* \f]
*
* The middle bonds are:
* TODO \f[ B_i = J \, \sigma^x_i \otimes \sigma^x_{i+1} + \frac{g}{2} \,\sigma^z_i \otimes I + \frac{g}{2} \, I \otimes \sigma^z_{i+1} \f]
* \f[
* -t \left(
* \hat{a}^\dagger_i \otimes \hat{a}_{i+1}
* + \hat{a}_i \otimes \hat{a}^\dagger_{i+1}
* \right)
* + \frac{U}{4} \Big(
* \hat{n}_{i} \left(\hat{n}_{i}-1_{i}\right) \otimes \mathbbm{1}_{i+1}
* + \mathbbm{1}_{i} \otimes \hat{n}_{i+1} \left(\hat{n}_{i+1}-1_{i+1} \right)
* \Big)
* - \frac{\mu}{2} \left(
* \hat{n}_i \otimes \mathbbm{1}_{i+1}
* \mathbbm{1}_{i} \otimes \hat{n}_{i+1}
* \right)
* \f]
*
* The MPO is:
* \f[ W = \begin{pmatrix} I & \hat{a}^\dagger & \hat{a} & \frac{U}{2}\hat{n}(\hat{n} - I) - \mu\hat{n} \\ 0 & 0 & 0 & -t \hat{a} \\ 0 & 0 & 0 & -t \hat{a}^\dagger\\ 0 & 0 & 0 & I \end{pmatrix} \f]
* \f[ W = \begin{pmatrix}
* \mathbbm{1} & \hat{a}^\dagger & \hat{a} & \frac{U}{2}\hat{n}(\hat{n} - \mathbbm{1}) - \mu\hat{n} \\
* 0 & 0 & 0 & -t \hat{a} \\ 0 & 0 & 0 & -t \hat{a}^\dagger\\
* 0 & 0 & 0 & \mathbbm{1} \end{pmatrix}
* \f]
*/
template<typename T>
class BoseHubbardModel1D : public Model1D<T> {
@ -43,7 +74,6 @@ class BoseHubbardModel1D : public Model1D<T> {
virtual std::string format() const override {
return std::format("BoseHubbardModel1D(L={},U={},t={},mu={})", this->L, this->U, this->t, this->mu);
}
protected:
void initH_Bonds() override;
void initH_MPO() override;
@ -53,5 +83,3 @@ class BoseHubbardModel1D : public Model1D<T> {
T mu;
LadderOperators<T> lo;
};

22
src/models/transverse_field_ising.cpp
View File

@ -8,22 +8,22 @@ PauliMatrices<T> TransverseFieldIsingModel1D<T>::pm;
template<typename T>
void TransverseFieldIsingModel1D<T>::initH_Bonds() {
T gHalf = this->g * static_cast<T>(0.5);
T BHalf = this->B * static_cast<T>(0.5);
// TODO: improve code or directly initialize the tensors
auto left = std::make_shared<e::MatrixX<T>>(
e::kroneckerProduct(pm.getSigmaX(), pm.getSigmaX()) * (-J) -
this->g * e::kroneckerProduct(pm.getSigmaZ(), pm.getIdentity()) -
gHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaZ())
e::kroneckerProduct(pm.getSigmaZ(), pm.getSigmaZ()) * (-J)
- this->B * e::kroneckerProduct(pm.getSigmaX(), pm.getIdentity())
- BHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaX())
);
auto right = std::make_shared<e::MatrixX<T>>(
e::kroneckerProduct(pm.getSigmaX(), pm.getSigmaX()) * (-J) -
gHalf * e::kroneckerProduct(pm.getSigmaZ(), pm.getIdentity()) -
this->g * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaZ())
e::kroneckerProduct(pm.getSigmaZ(), pm.getSigmaZ()) * (-J)
- BHalf * e::kroneckerProduct(pm.getSigmaX(), pm.getIdentity())
- this->B * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaX())
);
auto middle = std::make_shared<e::MatrixX<T>>(
e::kroneckerProduct(pm.getSigmaX(), pm.getSigmaX()) * (-J) -
gHalf * e::kroneckerProduct(pm.getSigmaZ(), pm.getIdentity()) -
gHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaZ())
e::kroneckerProduct(pm.getSigmaZ(), pm.getSigmaZ()) * (-J)
- BHalf * e::kroneckerProduct(pm.getSigmaX(), pm.getIdentity())
- BHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaX())
);
auto leftT = std::make_shared<e::Tensor<T, 4>>(2, 2, 2, 2);
auto middleT = std::make_shared<e::Tensor<T, 4>>(2, 2, 2, 2);
@ -63,7 +63,7 @@ void TransverseFieldIsingModel1D<T>::initH_MPO() {
e::Matrix2<T> zero(2, 2);
zero.setZero();
auto W = std::make_shared<e::Tensor<T, 4>>(std::move(constructMPO<T>({
{this->pm.getIdentity(), this->pm.getSigmaZ(), (this->g * this->pm.getSigmaX())},
{this->pm.getIdentity(), this->pm.getSigmaZ(), (-this->B) * this->pm.getSigmaX()},
{zero, zero, (-this->J * this->pm.getSigmaZ())},
{zero, zero, this->pm.getIdentity()}
})));

22
src/models/transverse_field_ising.hpp
View File

@ -4,25 +4,33 @@
/**
* @brief Transverse Field Ising Model in 1 dimension
* @details
* The nearest neighbor interaction is in z direction, the magnetic field in x direction.
* The Hamiltonian is:
* \f[ H = -J \sum_{\langle i,j\rangle} \sigma_{z,i} \sigma_{z,j} - B \sum_i \sigma_{x,j} \f]
*
* The middle bonds are:
* \f[ B_i = J \, \sigma^x_i \otimes \sigma^x_{i+1} + \frac{g}{2} \,\sigma^z_i \otimes I + \frac{g}{2} \, I \otimes \sigma^z_{i+1} \f]
* \f[ B_i = J \, \sigma^z_i \otimes \sigma^z_{i+1} + \frac{B}{2} \,\sigma^x_i \otimes I + \frac{B}{2} \, I \otimes \sigma^x_{i+1} \f]
*
* The MPO is:
* \f[ W = \begin{pmatrix} I & \sigma_z & g \, \sigma_x \\ 0 & 0 & J\,\sigma_z \\ 0 & 0 & I \end{pmatrix} \f]
* \f[ W = \begin{pmatrix}
* \mathbbm{1} & \sigma_z & -B \, \sigma_x \\
* 0 & 0 & -J\,\sigma_z \\
* 0 & 0 & \mathbbm{1}
* \end{pmatrix}
* \f]
*/
template<typename T>
class TransverseFieldIsingModel1D : public Model1D<T> {
public:
TransverseFieldIsingModel1D(size_t L, T J, T g) : Model1D<T>(L, 2), J(J), g(g) {
// std::println("TFIModel: L={}, J={}, g={}", L, J, g);
TransverseFieldIsingModel1D(size_t L, T J, T B) : Model1D<T>(L, 2), J(J), B(B) {
this->initH_Bonds();
this->initH_MPO();
}
inline const T& getJ() const { return J; }
inline const T& getg() const { return g; }
inline const T& getB() const { return B; }
virtual std::string format() const override {
return std::format("TransverseFieldIsingModel1D(L={},J={},g={})", this->L, this->J, this->g);
return std::format("TransverseFieldIsingModel1D(L={},J={},B={})", this->L, this->J, this->B);
}
protected:
void initH_Bonds() override;
@ -30,7 +38,7 @@ class TransverseFieldIsingModel1D : public Model1D<T> {
void initH_MPO() override;
private:
T J;
T g;
T B;
static PauliMatrices<T> pm;
};

154
src/mpl_colors.py
View File

@ -1,154 +0,0 @@
import matplotlib as mpl
import matplotlib.pyplot as plt
from cycler import cycler
# Palette
P = {
"bg0": "#282828",
"bg0-hard": "#1d2021",
"bg0-soft": "#32302f",
"bg1": "#3c3836",
"bg2": "#504945",
"bg3": "#665c54",
"bg4": "#7c6f64",
"fg0": "#fbf1c7",
"fg0-hard": "#f9f5d7",
"fg0-soft": "#f2e5bc",
"fg1": "#ebdbb2",
"fg2": "#d5c4a1",
"fg3": "#bdae93",
"fg4": "#a89984",
"dark-red": "#cc241d",
"dark-green": "#98971a",
"dark-yellow": "#d79921",
"dark-blue": "#458588",
"dark-purple": "#b16286",
"dark-aqua": "#689d6a",
"dark-orange": "#d65d0e",
"dark-gray": "#928374",
"light-red": "#fb4934",
"light-green": "#b8bb26",
"light-yellow": "#fabd2f",
"light-blue": "#83a598",
"light-purple": "#d3869b",
"light-aqua": "#8ec07c",
"light-orange": "#f38019",
"light-gray": "#a89984",
"alt-red": "#9d0006",
"alt-green": "#79740e",
"alt-yellow": "#b57614",
"alt-blue": "#076678",
"alt-purple": "#8f3f71",
"alt-aqua": "#427b58",
"alt-orange": "#af3a03",
"alt-gray": "#7c6f64",
}
COLORS_LIGHT = [
P["light-blue"],
P["light-orange"],
P["light-green"],
P["light-red"],
P["light-purple"],
P["light-yellow"],
P["light-aqua"],
P["light-gray"]]
COLORS_DARK = [
P["dark-blue"],
P["dark-orange"],
P["dark-green"],
P["dark-red"],
P["dark-purple"],
P["dark-yellow"],
P["dark-aqua"],
P["dark-gray"]]
COLORS_ALT = [
P["alt-blue"],
P["alt-orange"],
P["alt-green"],
P["alt-red"],
P["alt-purple"],
P["alt-yellow"],
P["alt-aqua"],
P["alt-gray"]]
COLORS = COLORS_LIGHT
def set_colorscheme(variant="dark"):
global COLORS
if variant == "dark":
FIG_BG_COLOR = P["bg0"]
PLT_FG_COLOR = P["fg0"]
PLT_BG_COLOR = P["bg0"]
PLT_GRID_COLOR = P["bg2"]
LEGEND_FG_COLOR = PLT_FG_COLOR
LEGEND_BG_COLOR = P["bg1"]
LEGEND_BORDER_COLOR = P["bg2"]
COLORS = COLORS_LIGHT
else:
FIG_BG_COLOR = P["fg0"]
PLT_FG_COLOR = P["bg0"]
PLT_BG_COLOR = P["fg0"]
PLT_GRID_COLOR = P["fg2"]
LEGEND_FG_COLOR = PLT_FG_COLOR
LEGEND_BG_COLOR = P["fg1"]
LEGEND_BORDER_COLOR = P["fg2"]
COLORS = COLORS_DARK
if variant == "alt":
COLORS = COLORS_ALT
color_rcParams = {
'axes.edgecolor': PLT_FG_COLOR,
'axes.facecolor': PLT_BG_COLOR,
'axes.labelcolor': PLT_FG_COLOR,
'axes.titlecolor': 'auto',
# 'axes.prop_cycle': cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']),
'axes.prop_cycle': cycler('color', COLORS),
# 'axes3d.xaxis.panecolor': (0.95, 0.95, 0.95, 0.5),
# 'axes3d.yaxis.panecolor': (0.9, 0.9, 0.9, 0.5),
# 'axes3d.zaxis.panecolor': (0.925, 0.925, 0.925, 0.5),
# 'boxplot.boxprops.color': 'black',
# 'boxplot.capprops.color': 'black',
# 'boxplot.flierprops.color': 'black',
# 'boxplot.flierprops.markeredgecolor': 'black',
# 'boxplot.flierprops.markeredgewidth': 1.0,
# 'boxplot.flierprops.markerfacecolor': 'none',
# 'boxplot.meanprops.color': 'C2',
# 'boxplot.meanprops.markeredgecolor': 'C2',
# 'boxplot.meanprops.markerfacecolor': 'C2',
# 'boxplot.meanprops.markersize': 6.0,
# 'boxplot.medianprops.color': 'C1',
# 'boxplot.whiskerprops.color': 'black',
'figure.edgecolor': PLT_BG_COLOR,
'figure.facecolor': PLT_BG_COLOR,
# 'figure.figsize': [6.4, 4.8],
# 'figure.frameon': True,
# 'figure.labelsize': 'large',
'grid.color': PLT_GRID_COLOR,
# 'hatch.color': 'black',
'legend.edgecolor': LEGEND_BORDER_COLOR,
'legend.facecolor': LEGEND_BG_COLOR,
'xtick.color': PLT_FG_COLOR,
'ytick.color': PLT_FG_COLOR,
'xtick.labelcolor': PLT_FG_COLOR,
'ytick.labelcolor': PLT_FG_COLOR,
# 'lines.color': 'C0',
'text.color': PLT_FG_COLOR,
}
for k, v in color_rcParams.items():
plt.rcParams[k] = v
# override single char codes
mpl.colors.get_named_colors_mapping()["b"] = COLORS[0]
mpl.colors.get_named_colors_mapping()["o"] = COLORS[1]
mpl.colors.get_named_colors_mapping()["g"] = COLORS[2]
mpl.colors.get_named_colors_mapping()["r"] = COLORS[3]
mpl.colors.get_named_colors_mapping()["m"] = COLORS[4]
mpl.colors.get_named_colors_mapping()["y"] = COLORS[5]
mpl.colors.get_named_colors_mapping()["c"] = COLORS[6]
mpl.colors.get_named_colors_mapping()["k"] = P["fg0"]
mpl.colors.get_named_colors_mapping()["w"] = P["bg0"]

32
src/mps.cpp
View File

@ -36,6 +36,24 @@ std::vector<unsigned> MPS<T>::getChi() const {
}
return chis;
}
template<typename T>
unsigned MPS<T>::getChiMax() const {
unsigned chiMax = this->B(0).dimension(2);
for (size_t i = 1; i < (L - 1); i++) {
if (this->B(i).dimension(2) > chiMax) {
chiMax = this->B(i).dimension(2);
}
}
return chiMax;
}
template<typename T>
double MPS<T>::getChiAverage() const {
double chiAvg = 0;
for (size_t i = 0; i < (L - 1); i++) {
chiAvg += this->B(i).dimension(2);
}
return chiAvg / (L-1);
}
template<typename T>
T MPS<T>::getSiteExpectationValue(const e::Tensor<T, 2>& op, size_t siteIdx) const {
@ -98,6 +116,7 @@ T MPS<T>::collapse(size_t siteIndex) const {
e::Tensor<dtype, 2> thetaTheta = theta.contract(theta.conjugate(), e::array<iP, 2>{iP{0, 0}, iP{1, 1}});
for (size_t i = siteIndex+1; i < this->L; i++) {
// vL vR vL* vR*
// TODO optimize by not assigning to tensor
e::Tensor<dtype, 4> BB = this->B(i).contract(this->B(i).conjugate(), e::array<iP, 1>{iP{1, 1}});
e::Tensor<dtype, 2> thetaTheta2 = thetaTheta.contract(BB, e::array<iP, 2>{iP{0, 0}, iP{1, 2}});
std::swap(thetaTheta, thetaTheta2);
@ -139,15 +158,20 @@ split_truncate_theta(const e::Tensor<dtype, 4>& theta, size_t chiMax, dtype eps)
e::Index dL = theta.dimension(1);
e::Index dR = theta.dimension(2);
e::Index chivR = theta.dimension(3);
// Reshape theta to 2D matrix for the SVD
using MatX = e::Matrix<dtype, -1, -1>;
/// todo (DONE): Possible optimization: either use a Map and later adjust the matrix to tensor loops, or globally use RowMajor memory layout?
// The map approach produces the same Matrix is tensor2Matrix with Variant 2
MatX theta_mat = e::Map<const MatX>(theta.data(), chivL * dL, dR * chivR);
// MatX theta_mat = tensor2Matrix(theta);
return split_truncate_theta(theta_mat, chivL, dL, dR, chivR, chiMax, eps);
}
template<typename dtype>
std::tuple<e::Tensor<dtype, 3>, e::Tensor<dtype, 1>, e::Tensor<dtype, 3>>
split_truncate_theta(const e::MatrixX<dtype>& theta, e::Index chivL, e::Index dL, e::Index dR, e::Index chivR, size_t chiMax, dtype eps) {
using MatX = e::Matrix<dtype, -1, -1>;
// std::cout << "theta:\n" << theta << "theta_mat\n": << theta_mat << std::endl;
// std::cout << "theta_mat = " << std::endl;
@ -155,7 +179,7 @@ split_truncate_theta(const e::Tensor<dtype, 4>& theta, size_t chiMax, dtype eps)
// perform SVD: theta = U * diag(S) * V^T
constexpr int SVDOptions = e::ComputeThinU | e::ComputeThinV;
e::JacobiSVD<MatX, SVDOptions> svd(theta_mat);
e::BDCSVD<MatX, SVDOptions> svd(theta);
// e::BDCSVD<MatX, SVDOptions> svd(theta_mat);
const MatX& U = svd.matrixU();
const e::VectorX<dtype>& S = svd.singularValues();

85
src/mps.hpp
View File

@ -1,9 +1,11 @@
#pragma once
#include "util.hpp"
#include <cmath>
#include <eigen3/Eigen/Eigen>
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
#include <vector>
#include <cassert>
#include <random>
#include <iostream>
@ -49,6 +51,10 @@ class MPS {
* @brief Return the dimensions of all bonds
*/
std::vector<unsigned> getChi() const;
/// @brief Return the maximum bond dimension
unsigned getChiMax() const;
/// @brief Return the average bond dimension
double getChiAverage() const;
/**
* @brief Compute the expectation value of a single-site operator on a single site
@ -75,6 +81,7 @@ class MPS {
STensor_t& S(size_t i) { return this->Ss.at(i); }
const BTensor_t& B(size_t i) const { return this->Bs.at(i); }
const STensor_t& S(size_t i) const { return this->Ss.at(i); }
// size_t getPhysicalDimension() const { return this->Bs.at(0).dimension(1); }
size_t getL() const { return this->L; }
void print() const {
for (size_t i = 0; i < this->L; i++) {
@ -98,24 +105,49 @@ class MPS {
* @returns (A, S, B) where A is the left-canonical matrix on site i (legs `vL, i, vC`),
* S are the singular values and
* B is the right-canonical matrix on site i+1, with legs `vC, j, vR`
* @note This functions maps the tensor and calls the Matrux overload of split_truncate_theta internally.
*/
template<typename dtype>
std::tuple<e::Tensor<dtype, 3>, e::Tensor<dtype, 1>, e::Tensor<dtype, 3>>
split_truncate_theta(const e::Tensor<dtype, 4>& theta, size_t chiMax, dtype eps);
template<typename dtype>
std::tuple<e::Tensor<dtype, 3>, e::Tensor<dtype, 1>, e::Tensor<dtype, 3>>
split_truncate_theta(const e::MatrixX<dtype>& theta, e::Index chivL, e::Index dL, e::Index dR, e::Index chivR, size_t chiMax, dtype eps);
template<typename T>
void applyPertubationOnB(std::vector<e::Tensor<T, 3>>& Bs, double pertubationStrength) {
std::mt19937 generator(493852934); // Mersenne twister rng
std::uniform_real_distribution<double> distribution(0.0, pertubationStrength);
for (auto it = Bs.begin(); it != Bs.end(); it++) {
double abs = 0;
for (unsigned i = 0; i < it->dimension(1); i++) {
(*it)(0, i, 0) += distribution(generator);
abs += std::pow((*it)(0, i, 0).real(), 2);
}
abs = std::sqrt(abs);
for (unsigned i = 0; i < it->dimension(1); i++) {
(*it)(0, i, 0) /= abs;
}
}
}
/**
* @brief Initialize a MPS in spin-up state
* @brief Initialize a MPS using copies of the B tensor
* @details
* Return a MPS spin-up MPS with length L and bond dimension 1
* A random pertubation up to `pertubationStrength` can be applied to each occupation number in the B tensor,
* the tensors are normalized after the pertubations are applied.
* This can help break symmetries, which can make get the simulations stuck.
*/
template<typename T>
MPS<T> initMPS_spinUp(size_t L) {
assert(L > 0);
e::Tensor<T, 3> B(1, 2, 1);
B.setZero();
B(0, 0, 0) = 1;
MPS<T> makeMPS(const e::Tensor<T, 3>& B, size_t L, double pertubationStrength) {
std::vector<e::Tensor<T, 3>> Bs(L, B);
if (pertubationStrength > 0) {
applyPertubationOnB(Bs, pertubationStrength);
}
e::Tensor<T, 1> S(1);
S(0) = 1;
std::vector<e::Tensor<T, 1>> Ss(L, S);
@ -123,24 +155,43 @@ MPS<T> initMPS_spinUp(size_t L) {
}
/**
* @brief Initialize a MPS in spin-up state
* @details
* Return a MPS spin-up MPS with length L and bond dimension 1
*/
template<typename T>
MPS<T> initMPS_spinUp(size_t L, double pertubationStrength=0) {
assert(L > 0);
e::Tensor<T, 3> B(1, 2, 1);
B.setZero();
B(0, 0, 0) = 1;
return makeMPS(B, L, pertubationStrength);
}
template<typename T>
MPS<T> initMPS_spinDown(size_t L, double pertubationStrength=0) {
assert(L > 0);
e::Tensor<T, 3> B(1, 2, 1);
B.setZero();
B(0, 1, 0) = 1;
return makeMPS(B, L, pertubationStrength);
}
/**
* @brief Initialize a MPS in spin-right state
* @details
* Return a MPS spin-right MPS with length L and bond dimension 1
*/
template<typename T>
MPS<T> initMPS_spinRight(size_t L) {
MPS<T> initMPS_spinRight(size_t L, double pertubationStrength=0) {
assert(L > 0);
e::Tensor<T, 3> B(1, 2, 1);
B.setZero();
T value = std::sqrt(0.5);
B(0, 0, 0) = value;
B(0, 1, 0) = value;
std::vector<e::Tensor<T, 3>> Bs(L, B);
e::Tensor<T, 1> S(1);
S(0) = 1;
std::vector<e::Tensor<T, 1>> Ss(L, S);
return MPS<T>(std::move(Bs), std::move(Ss));
return makeMPS(B, L, pertubationStrength);
}
@ -150,15 +201,11 @@ MPS<T> initMPS_spinRight(size_t L) {
* Return a MPS of dimension `nMax+1`, where each site is occupied by `nParticles`.
*/
template<typename T>
MPS<T> initMPS_nParticles(size_t L, size_t nMax, size_t nParticles) {
MPS<T> initMPS_nParticles(size_t L, size_t nMax, size_t nParticles, double pertubationStrength=0) {
assert(L > 0);
assert(nParticles <= nMax);
e::Tensor<T, 3> B(1, nMax+1, 1);
B.setZero();
B(0, nParticles, 0) = 1;
std::vector<e::Tensor<T, 3>> Bs(L, B);
e::Tensor<T, 1> S(1);
S(0) = 1;
std::vector<e::Tensor<T, 1>> Ss(L, S);
return MPS<T>(std::move(Bs), std::move(Ss));
return makeMPS(B, L, pertubationStrength);
}

562
src/output
View File

@ -1,562 +0,0 @@
Building: main.cpp -> ../build/main.o
In file included from main.cpp:3:
mps.hpp: In instantiation of 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]':
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
mps.hpp:140:83: error: no match for 'operator>' (operand types are 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' and 'std::complex<float>')
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/bits/stl_algobase.h:64,
from /usr/include/c++/14.2.1/vector:62,
from ./pch.hpp:1,
from <command-line>:
/usr/include/c++/14.2.1/bits/stl_pair.h:1023:5: note: candidate: 'template<class _T1, class _T2, class _U1, class _U2> constexpr std::common_comparison_category_t<decltype (std::__detail::__synth3way(declval<_T1&>(), declval<_U1&>())), decltype (std::__detail::__synth3way(declval<_T2&>(), declval<_U2&>()))> std::operator<=>(const pair<_T1, _T2>&, const pair<_U1, _U2>&)' (reversed)
1023 | operator<=>(const pair<_T1, _T2>& __x, const pair<_U1, _U2>& __y)
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_pair.h:1023:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::pair<_T1, _T2>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/bits/stl_algobase.h:67:
/usr/include/c++/14.2.1/bits/stl_iterator.h:576:5: note: candidate: 'template<class _IteratorL, class _IteratorR> requires three_way_comparable_with<_IteratorR, _IteratorL, std::partial_ordering> constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const reverse_iterator<_IteratorL>&, const reverse_iterator<_IteratorR>&)' (reversed)
576 | operator<=>(const reverse_iterator<_IteratorL>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:576:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::reverse_iterator<_IteratorL>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:1679:5: note: candidate: 'template<class _IteratorL, class _IteratorR> requires three_way_comparable_with<_IteratorR, _IteratorL, std::partial_ordering> constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const move_iterator<_IteratorL>&, const move_iterator<_IteratorR>&)' (reversed)
1679 | operator<=>(const move_iterator<_IteratorL>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:1679:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::move_iterator<_IteratorL>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/bits/uses_allocator_args.h:39,
from /usr/include/c++/14.2.1/bits/memory_resource.h:41,
from /usr/include/c++/14.2.1/vector:87:
/usr/include/c++/14.2.1/tuple:2589:5: note: candidate: 'template<class ... _Tps, class ... _Ups> constexpr std::common_comparison_category_t<decltype (std::__detail::__synth3way(declval<_Tps&>(), declval<_Ups&>()))...> std::operator<=>(const tuple<_Ts ...>&, const tuple<_Elements ...>&)' (reversed)
2589 | operator<=>(const tuple<_Tps...>& __t, const tuple<_Ups...>& __u)
| ^~~~~~~~
/usr/include/c++/14.2.1/tuple:2589:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::tuple<_Ts ...>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/bits/basic_string.h:47,
from /usr/include/c++/14.2.1/string:54,
from /usr/include/c++/14.2.1/bits/locale_classes.h:40,
from /usr/include/c++/14.2.1/locale:41,
from /usr/include/c++/14.2.1/format:47,
from /usr/include/c++/14.2.1/print:41,
from ./pch.hpp:2:
/usr/include/c++/14.2.1/string_view:617:5: note: candidate: 'template<class _CharT, class _Traits> constexpr decltype (__char_traits_cmp_cat<_Traits>(0)) std::operator<=>(basic_string_view<_CharT, _Traits>, __type_identity_t<basic_string_view<_CharT, _Traits> >)' (reversed)
617 | operator<=>(basic_string_view<_CharT, _Traits> __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/string_view:617:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'std::basic_string_view<_CharT, _Traits>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/basic_string.h:3805:5: note: candidate: 'template<class _CharT, class _Traits, class _Alloc> constexpr decltype (__char_traits_cmp_cat<_Traits>(0)) std::operator<=>(const __cxx11::basic_string<_CharT, _Traits, _Allocator>&, const _CharT*)' (reversed)
3805 | operator<=>(const basic_string<_CharT, _Traits, _Alloc>& __lhs,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/basic_string.h:3805:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::__cxx11::basic_string<_CharT, _Traits, _Allocator>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/format:48:
/usr/include/c++/14.2.1/optional:1283:5: note: candidate: 'template<class _Tp, class _Up> requires three_way_comparable_with<_Up, _Tp, std::partial_ordering> constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const optional<_Tp>&, const optional<_Up>&)' (reversed)
1283 | operator<=>(const optional<_Tp>& __x, const optional<_Up>& __y)
| ^~~~~~~~
/usr/include/c++/14.2.1/optional:1283:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::optional<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/optional:1298:5: note: candidate: 'template<class _Tp> constexpr std::strong_ordering std::operator<=>(const optional<_Tp>&, nullopt_t)' (reversed)
1298 | operator<=>(const optional<_Tp>& __x, nullopt_t) noexcept
| ^~~~~~~~
/usr/include/c++/14.2.1/optional:1298:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::optional<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/optional:1435:5: note: candidate: 'template<class _Tp, class _Up> requires !(__is_optional_v<_Up>) && (three_way_comparable_with<_Up, _Tp, std::partial_ordering>) constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const optional<_Tp>&, const _Up&)' (reversed)
1435 | operator<=>(const optional<_Tp>& __x, const _Up& __v)
| ^~~~~~~~
/usr/include/c++/14.2.1/optional:1435:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::optional<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/bits/shared_ptr_base.h:59,
from /usr/include/c++/14.2.1/bits/shared_ptr.h:53,
from /usr/include/c++/14.2.1/chrono:49,
from /usr/include/eigen3/unsupported/Eigen/CXX11/Tensor:38,
from ./pch.hpp:4:
/usr/include/c++/14.2.1/bits/unique_ptr.h:996:5: note: candidate: 'template<class _Tp, class _Dp, class _Up, class _Ep> requires three_way_comparable_with<typename std::unique_ptr<_Tp, _Dp>::pointer, typename std::unique_ptr<_Up, _Ep>::pointer, std::partial_ordering> constexpr std::compare_three_way_result_t<typename std::unique_ptr<_Tp, _Dp>::pointer, typename std::unique_ptr<_Up, _Ep>::pointer> std::operator<=>(const unique_ptr<_Tp, _Dp>&, const unique_ptr<_Up, _Ep>&)' (reversed)
996 | operator<=>(const unique_ptr<_Tp, _Dp>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:996:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::unique_ptr<_Tp, _Dp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:1005:5: note: candidate: 'template<class _Tp, class _Dp> requires three_way_comparable<typename std::unique_ptr<_Tp, _Dp>::pointer, std::partial_ordering> constexpr std::compare_three_way_result_t<typename std::unique_ptr<_Tp, _Dp>::pointer> std::operator<=>(const unique_ptr<_Tp, _Dp>&, nullptr_t)' (reversed)
1005 | operator<=>(const unique_ptr<_Tp, _Dp>& __x, nullptr_t)
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:1005:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::unique_ptr<_Tp, _Dp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1809:5: note: candidate: 'template<class _Tp, class _Up, __gnu_cxx::_Lock_policy _Lp> std::strong_ordering std::operator<=>(const __shared_ptr<_Tp1, _Lp>&, const __shared_ptr<_Tp2, _Lp>&)' (reversed)
1809 | operator<=>(const __shared_ptr<_Tp, _Lp>& __a,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1809:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::__shared_ptr<_Tp1, _Lp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1815:5: note: candidate: 'template<class _Tp, __gnu_cxx::_Lock_policy _Lp> std::strong_ordering std::operator<=>(const __shared_ptr<_Tp, _Lp>&, nullptr_t)' (reversed)
1815 | operator<=>(const __shared_ptr<_Tp, _Lp>& __a, nullptr_t) noexcept
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1815:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::__shared_ptr<_Tp, _Lp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr.h:566:5: note: candidate: 'template<class _Tp, class _Up> std::strong_ordering std::operator<=>(const shared_ptr<_Tp>&, const shared_ptr<_Tp>&)' (reversed)
566 | operator<=>(const shared_ptr<_Tp>& __a,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr.h:566:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::shared_ptr<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr.h:572:5: note: candidate: 'template<class _Tp> std::strong_ordering std::operator<=>(const shared_ptr<_Tp>&, nullptr_t)' (reversed)
572 | operator<=>(const shared_ptr<_Tp>& __a, nullptr_t) noexcept
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/shared_ptr.h:572:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::shared_ptr<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:594:5: note: candidate: 'template<class _Iterator> requires three_way_comparable<_Iterator, std::partial_ordering> constexpr std::compare_three_way_result_t<_Iterator, _Iterator> std::operator<=>(const reverse_iterator<_IteratorL>&, const reverse_iterator<_IteratorL>&)' (rewritten)
594 | operator<=>(const reverse_iterator<_Iterator>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:594:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::reverse_iterator<_IteratorL>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:1745:5: note: candidate: 'template<class _Iterator> requires three_way_comparable<_Iterator, std::partial_ordering> constexpr std::compare_three_way_result_t<_Iterator, _Iterator> std::operator<=>(const move_iterator<_IteratorL>&, const move_iterator<_IteratorL>&)' (rewritten)
1745 | operator<=>(const move_iterator<_Iterator>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:1745:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::move_iterator<_IteratorL>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/vector:66:
/usr/include/c++/14.2.1/bits/stl_vector.h:2069:5: note: candidate: 'template<class _Tp, class _Alloc> constexpr std::__detail::__synth3way_t<_Iterator> std::operator<=>(const vector<_Tp, _Alloc>&, const vector<_Tp, _Alloc>&)' (rewritten)
2069 | operator<=>(const vector<_Tp, _Alloc>& __x, const vector<_Tp, _Alloc>& __y)
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_vector.h:2069:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::vector<_Tp, _Alloc>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/format:43:
/usr/include/c++/14.2.1/array:310:5: note: candidate: 'template<class _Tp, long unsigned int _Nm> constexpr std::__detail::__synth3way_t<_Iterator> std::operator<=>(const array<_Tp, _Nm>&, const array<_Tp, _Nm>&)' (rewritten)
310 | operator<=>(const array<_Tp, _Nm>& __a, const array<_Tp, _Nm>& __b)
| ^~~~~~~~
/usr/include/c++/14.2.1/array:310:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::array<_Tp, _Nm>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/basic_string.h:3790:5: note: candidate: 'template<class _CharT, class _Traits, class _Alloc> constexpr decltype (__char_traits_cmp_cat<_Traits>(0)) std::operator<=>(const __cxx11::basic_string<_CharT, _Traits, _Allocator>&, const __cxx11::basic_string<_CharT, _Traits, _Allocator>&)' (rewritten)
3790 | operator<=>(const basic_string<_CharT, _Traits, _Alloc>& __lhs,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/basic_string.h:3790:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::__cxx11::basic_string<_CharT, _Traits, _Allocator>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/format:52:
/usr/include/c++/14.2.1/variant:1280:5: note: candidate: 'template<class ... _Types> requires (three_way_comparable<_Types, std::partial_ordering> && ...) constexpr std::common_comparison_category_t<typename std::__detail::__cmp3way_res_impl<_Types, _Types>::type ...> std::operator<=>(const variant<_Types ...>&, const variant<_Types ...>&)' (rewritten)
1280 | operator<=>(const variant<_Types...>& __v, const variant<_Types...>& __w)
| ^~~~~~~~
/usr/include/c++/14.2.1/variant:1280:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::variant<_Types ...>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/map:63,
from /usr/include/eigen3/Eigen/SparseCore:16,
from /usr/include/eigen3/Eigen/Sparse:26,
from /usr/include/eigen3/Eigen/Eigen:2,
from ./pch.hpp:5:
/usr/include/c++/14.2.1/bits/stl_map.h:1533:5: note: candidate: 'template<class _Key, class _Tp, class _Compare, class _Alloc> std::__detail::__synth3way_t<std::pair<const _Key, _Val> > std::operator<=>(const map<_Key, _Tp, _Compare, _Allocator>&, const map<_Key, _Tp, _Compare, _Allocator>&)' (rewritten)
1533 | operator<=>(const map<_Key, _Tp, _Compare, _Alloc>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_map.h:1533:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::map<_Key, _Tp, _Compare, _Allocator>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/map:64:
/usr/include/c++/14.2.1/bits/stl_multimap.h:1155:5: note: candidate: 'template<class _Key, class _Tp, class _Compare, class _Alloc> std::__detail::__synth3way_t<std::pair<const _Key, _Val> > std::operator<=>(const multimap<_Key, _Tp, _Compare, _Allocator>&, const multimap<_Key, _Tp, _Compare, _Allocator>&)' (rewritten)
1155 | operator<=>(const multimap<_Key, _Tp, _Compare, _Alloc>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_multimap.h:1155:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::multimap<_Key, _Tp, _Compare, _Allocator>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/list:65,
from /usr/include/eigen3/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h:15,
from /usr/include/eigen3/Eigen/IterativeLinearSolvers:47,
from /usr/include/eigen3/Eigen/Sparse:31:
/usr/include/c++/14.2.1/bits/stl_list.h:2158:5: note: candidate: 'template<class _Tp, class _Alloc> std::__detail::__synth3way_t<_Iterator> std::operator<=>(const __cxx11::list<_Tp, _Alloc>&, const __cxx11::list<_Tp, _Alloc>&)' (rewritten)
2158 | operator<=>(const list<_Tp, _Alloc>& __x, const list<_Tp, _Alloc>& __y)
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_list.h:2158:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::__cxx11::list<_Tp, _Alloc>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
In file included from /usr/include/c++/14.2.1/bits/ios_base.h:46,
from /usr/include/c++/14.2.1/bits/locale_facets.h:43,
from /usr/include/c++/14.2.1/locale:42:
/usr/include/c++/14.2.1/system_error:316:3: note: candidate: 'std::strong_ordering std::operator<=>(const error_code&, const error_code&)' (rewritten)
316 | operator<=>(const error_code& __lhs, const error_code& __rhs) noexcept
| ^~~~~~~~
/usr/include/c++/14.2.1/system_error:316:33: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' to 'const std::error_code&'
316 | operator<=>(const error_code& __lhs, const error_code& __rhs) noexcept
| ~~~~~~~~~~~~~~~~~~^~~~~
/usr/include/c++/14.2.1/system_error:498:3: note: candidate: 'std::strong_ordering std::operator<=>(const error_condition&, const error_condition&)' (rewritten)
498 | operator<=>(const error_condition& __lhs,
| ^~~~~~~~
/usr/include/c++/14.2.1/system_error:498:38: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' to 'const std::error_condition&'
498 | operator<=>(const error_condition& __lhs,
| ~~~~~~~~~~~~~~~~~~~~~~~^~~~~
/usr/include/c++/14.2.1/variant:1303:3: note: candidate: 'constexpr std::strong_ordering std::operator<=>(monostate, monostate)' (rewritten)
1303 | operator<=>(monostate, monostate) noexcept { return strong_ordering::equal; }
| ^~~~~~~~
/usr/include/c++/14.2.1/variant:1303:15: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' to 'std::monostate'
1303 | operator<=>(monostate, monostate) noexcept { return strong_ordering::equal; }
| ^~~~~~~~~
In file included from /usr/include/eigen3/unsupported/Eigen/CXX11/Tensor:45:
/usr/include/c++/14.2.1/thread:75:3: note: candidate: 'std::strong_ordering std::operator<=>(thread::id, thread::id)' (rewritten)
75 | operator<=>(thread::id __x, thread::id __y) noexcept
| ^~~~~~~~
/usr/include/c++/14.2.1/thread:75:26: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' to 'std::thread::id'
75 | operator<=>(thread::id __x, thread::id __y) noexcept
| ~~~~~~~~~~~^~~
In file included from /usr/include/eigen3/Eigen/src/Core/ArrayBase.h:99,
from /usr/include/eigen3/Eigen/Core:309,
from /usr/include/eigen3/unsupported/Eigen/CXX11/Tensor:14:
/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: candidate: 'template<class OtherDerived> const Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<typename OtherDerived::Scalar, typename Eigen::internal::traits<T>::Scalar, Eigen::internal::cmp_LT>, const OtherDerived, const Derived> Eigen::ArrayBase<Derived>::operator>(const Eigen::ArrayBase<OtherDerived>&) const [with Derived = Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >]'
204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const Eigen::ArrayBase<Derived>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: candidate: 'const Eigen::ArrayBase<Derived>::RCmpLTReturnType Eigen::ArrayBase<Derived>::operator>(const Scalar&) const [with Derived = Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >; RCmpLTReturnType = Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<double, double, Eigen::internal::cmp_LT, false>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Array<double, -1, 1> >, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> > >; Scalar = double]'
204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: no known conversion for argument 1 from 'std::complex<float>' to 'const Eigen::ArrayBase<Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> > >::Scalar&' {aka 'const double&'}
204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:551:5: note: candidate: 'template<class _IteratorL, class _IteratorR> constexpr bool std::operator>(const reverse_iterator<_IteratorL>&, const reverse_iterator<_IteratorR>&) requires requires{{std::operator>::__x->base() < std::operator>::__y->base()} -> decltype(auto) [requires std::convertible_to<<placeholder>, bool>];}'
551 | operator>(const reverse_iterator<_IteratorL>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:551:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::reverse_iterator<_IteratorL>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:1714:5: note: candidate: 'template<class _IteratorL, class _IteratorR> constexpr bool std::operator>(const move_iterator<_IteratorL>&, const move_iterator<_IteratorR>&) requires requires{{std::operator>::__y->base() < std::operator>::__x->base()} -> decltype(auto) [requires std::convertible_to<<placeholder>, bool>];}'
1714 | operator>(const move_iterator<_IteratorL>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/stl_iterator.h:1714:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::move_iterator<_IteratorL>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/optional:1258:5: note: candidate: 'template<class _Tp, class _Up> constexpr std::__optional_gt_t<_Tp, _Up> std::operator>(const optional<_Tp>&, const optional<_Up>&)'
1258 | operator>(const optional<_Tp>& __lhs, const optional<_Up>& __rhs)
| ^~~~~~~~
/usr/include/c++/14.2.1/optional:1258:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::optional<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/optional:1396:5: note: candidate: 'template<class _Tp, class _Up> constexpr std::__optional_gt_t<_Tp, _Up> std::operator>(const optional<_Tp>&, const _Up&)'
1396 | operator>(const optional<_Tp>& __lhs, const _Up& __rhs)
| ^~~~~~~~
/usr/include/c++/14.2.1/optional:1396:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::optional<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/optional:1402:5: note: candidate: 'template<class _Tp, class _Up> constexpr std::__optional_gt_t<_Up, _Tp> std::operator>(const _Up&, const optional<_Tp>&)'
1402 | operator>(const _Up& __lhs, const optional<_Tp>& __rhs)
| ^~~~~~~~
/usr/include/c++/14.2.1/optional:1402:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::optional<_Tp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/variant:1269:3: note: candidate: 'template<class ... _Types> constexpr bool std::operator>(const variant<_Types ...>&, const variant<_Types ...>&)'
1269 | _VARIANT_RELATION_FUNCTION_TEMPLATE(>, greater)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/c++/14.2.1/variant:1269:3: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::variant<_Types ...>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:942:5: note: candidate: 'template<class _Tp, class _Dp, class _Up, class _Ep> constexpr bool std::operator>(const unique_ptr<_Tp, _Dp>&, const unique_ptr<_Up, _Ep>&)'
942 | operator>(const unique_ptr<_Tp, _Dp>& __x,
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:942:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::unique_ptr<_Tp, _Dp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:950:5: note: candidate: 'template<class _Tp, class _Dp> constexpr bool std::operator>(const unique_ptr<_Tp, _Dp>&, nullptr_t)'
950 | operator>(const unique_ptr<_Tp, _Dp>& __x, nullptr_t)
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:950:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' is not derived from 'const std::unique_ptr<_Tp, _Dp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:960:5: note: candidate: 'template<class _Tp, class _Dp> constexpr bool std::operator>(nullptr_t, const unique_ptr<_Tp, _Dp>&)'
960 | operator>(nullptr_t, const unique_ptr<_Tp, _Dp>& __x)
| ^~~~~~~~
/usr/include/c++/14.2.1/bits/unique_ptr.h:960:5: note: template argument deduction/substitution failed:
mps.hpp:140:83: note: 'std::complex<float>' is not derived from 'const std::unique_ptr<_Tp, _Dp>'
140 | int chivC = static_cast<int>(std::min(chi_max, static_cast<size_t>((S.array() > eps).count())));
| ~~~~~~~~~~~^~~~~~
/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: candidate: 'const Eigen::ArrayBase<Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> > >::CmpLTReturnType Eigen::operator>(const ArrayBase<ArrayWrapper<const Matrix<double, -1, 1> > >::Scalar&, const ArrayWrapper<const Matrix<double, -1, 1> >&)'
204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> >' to 'const Eigen::ArrayBase<Eigen::ArrayWrapper<const Eigen::Matrix<double, -1, 1> > >::Scalar&' {aka 'const double&'}
204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/eigen3/Eigen/Core:300:
/usr/include/eigen3/Eigen/src/Core/Product.h: In instantiation of 'struct Eigen::internal::traits<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >':
/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1, 8>'
248 | class dense_product_base : public internal::dense_xpr_base<Product<Lhs, Rhs, Option>>::type {};
| ^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1, Eigen::Dense>'
276 | class ProductImpl<Lhs, Rhs, Option, Dense> : public internal::dense_product_base<Lhs, Rhs, Option> {
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>'
198 | class Product
| ^~~~~~~
mps.hpp:159:21: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
159 | MatX theta2 = U * S.asDiagonal() * V;
| ~~^~~~~~~~~~~~~~~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:33:90: error: no type named 'ReturnType' in 'struct Eigen::ScalarBinaryOpTraits<std::complex<float>, double, Eigen::internal::scalar_product_op<std::complex<float>, double> >'
33 | typename traits<RhsCleaned>::Scalar>::ReturnType Scalar;
| ^~~~~~
In file included from /usr/include/eigen3/Eigen/Core:296:
/usr/include/eigen3/Eigen/src/Core/DenseBase.h: In instantiation of 'class Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >':
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:52:7: required from 'class Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >'
52 | class MatrixBase : public DenseBase<Derived> {
| ^~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1, 8>'
248 | class dense_product_base : public internal::dense_xpr_base<Product<Lhs, Rhs, Option>>::type {};
| ^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1, Eigen::Dense>'
276 | class ProductImpl<Lhs, Rhs, Option, Dense> : public internal::dense_product_base<Lhs, Rhs, Option> {
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>'
198 | class Product
| ^~~~~~~
mps.hpp:159:21: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
159 | MatX theta2 = U * S.asDiagonal() * V;
| ~~^~~~~~~~~~~~~~~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
72 | using Base::coeff;
| ^~~~~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:73:15: error: 'coeffByOuterInner' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
73 | using Base::coeffByOuterInner;
| ^~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:81:24: error: 'operator()' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
81 | using Base::operator();
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:82:24: error: 'operator[]' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
82 | using Base::operator[];
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:88:15: error: 'w' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
88 | using Base::w;
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:89:15: error: 'x' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
89 | using Base::x;
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:90:15: error: 'y' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
90 | using Base::y;
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:91:15: error: 'z' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
91 | using Base::z;
| ^
In file included from /usr/include/eigen3/Eigen/Core:297:
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h: In instantiation of 'class Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >':
/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1, 8>'
248 | class dense_product_base : public internal::dense_xpr_base<Product<Lhs, Rhs, Option>>::type {};
| ^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1, Eigen::Dense>'
276 | class ProductImpl<Lhs, Rhs, Option, Dense> : public internal::dense_product_base<Lhs, Rhs, Option> {
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>'
198 | class Product
| ^~~~~~~
mps.hpp:159:21: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
159 | MatX theta2 = U * S.asDiagonal() * V;
| ~~^~~~~~~~~~~~~~~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
72 | using Base::coeff;
| ^~~~~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:77:15: error: 'eval' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
77 | using Base::eval;
| ^~~~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:81:23: error: 'operator-' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
81 | using Base::operator-;
| ^
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:84:23: error: 'operator*=' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
84 | using Base::operator*=;
| ^~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:85:23: error: 'operator/=' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1> >::Base'
85 | using Base::operator/=;
| ^~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h: In instantiation of 'class Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >':
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:52:7: required from 'class Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >'
52 | class MatrixBase : public DenseBase<Derived> {
| ^~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0, 8>'
248 | class dense_product_base : public internal::dense_xpr_base<Product<Lhs, Rhs, Option>>::type {};
| ^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0, Eigen::Dense>'
276 | class ProductImpl<Lhs, Rhs, Option, Dense> : public internal::dense_product_base<Lhs, Rhs, Option> {
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0>'
198 | class Product
| ^~~~~~~
mps.hpp:159:38: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
159 | MatX theta2 = U * S.asDiagonal() * V;
| ~~~~~~~~~~~~~~~~~~~^~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
72 | using Base::coeff;
| ^~~~~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:73:15: error: 'coeffByOuterInner' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
73 | using Base::coeffByOuterInner;
| ^~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:81:24: error: 'operator()' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
81 | using Base::operator();
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:82:24: error: 'operator[]' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
82 | using Base::operator[];
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:88:15: error: 'w' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
88 | using Base::w;
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:89:15: error: 'x' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
89 | using Base::x;
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:90:15: error: 'y' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
90 | using Base::y;
| ^
/usr/include/eigen3/Eigen/src/Core/DenseBase.h:91:15: error: 'z' has not been declared in 'Eigen::DenseBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
91 | using Base::z;
| ^
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h: In instantiation of 'class Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >':
/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0, 8>'
248 | class dense_product_base : public internal::dense_xpr_base<Product<Lhs, Rhs, Option>>::type {};
| ^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0, Eigen::Dense>'
276 | class ProductImpl<Lhs, Rhs, Option, Dense> : public internal::dense_product_base<Lhs, Rhs, Option> {
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0>'
198 | class Product
| ^~~~~~~
mps.hpp:159:38: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
159 | MatX theta2 = U * S.asDiagonal() * V;
| ~~~~~~~~~~~~~~~~~~~^~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
72 | using Base::coeff;
| ^~~~~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:77:15: error: 'eval' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
77 | using Base::eval;
| ^~~~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:81:23: error: 'operator-' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
81 | using Base::operator-;
| ^
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:84:23: error: 'operator*=' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
84 | using Base::operator*=;
| ^~
/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:85:23: error: 'operator/=' has not been declared in 'Eigen::MatrixBase<Eigen::Product<Eigen::Product<Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper<const Eigen::Matrix<double, -1, 1> >, 1>, Eigen::Matrix<std::complex<float>, -1, -1, 1, -1, -1>, 0> >::Base'
85 | using Base::operator/=;
| ^~
/usr/include/c++/14.2.1/format: In instantiation of 'class std::__format::_Checking_scanner<char, std::complex<float> >':
/usr/include/c++/14.2.1/format:4210:4: required from 'consteval std::basic_format_string<_CharT, _Args>::basic_format_string(const _Tp&) [with _Tp = char [4]; _CharT = char; _Args = {std::complex<float>&}]'
4210 | __scanner(_M_str);
| ^~~~~~~~~
mps.hpp:162:19: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
162 | std::print("{} ", theta2.data()[i]);
| ~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/c++/14.2.1/format:4065:10: error: static assertion failed: std::formatter must be specialized for each type being formatted
4065 | (is_default_constructible_v<formatter<_Args, _CharT>> && ...),
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/c++/14.2.1/format:4065:10: note: 'std::is_default_constructible_v<std::formatter<std::complex<float>, char> >' evaluates to false
/usr/include/c++/14.2.1/format: In instantiation of 'constexpr void std::__format::_Checking_scanner<_CharT, _Args>::_M_parse_format_spec(std::size_t) [with _Tp = std::complex<float>; _OtherArgs = {}; _CharT = char; _Args = {std::complex<float>}; std::size_t = long unsigned int]':
/usr/include/c++/14.2.1/format:4082:33: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
4082 | _M_parse_format_spec<_Args...>(__id);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
mps.hpp:162:19: in 'constexpr' expansion of 'std::basic_format_string<char, std::complex<float>&>("{} ")'
/usr/include/c++/14.2.1/format:4211:19: in 'constexpr' expansion of '__scanner.std::__format::_Checking_scanner<char, std::complex<float> >::std::__format::_Scanner<char>.std::__format::_Scanner<char>::_M_scan()'
/usr/include/c++/14.2.1/format:3952:30: in 'constexpr' expansion of '((std::__format::_Scanner<char>*)this)->std::__format::_Scanner<char>::_M_on_replacement_field()'
/usr/include/c++/14.2.1/format:4004:15: in 'constexpr' expansion of '((std::__format::_Scanner<char>*)this)->std::__format::_Scanner<char>::_M_format_arg(__id)'
/usr/include/c++/14.2.1/format:4095:38: error: use of deleted function 'std::formatter<_Tp, _CharT>::formatter() [with _Tp = std::complex<float>; _CharT = char]'
4095 | formatter<_Tp, _CharT> __f;
| ^~~
/usr/include/c++/14.2.1/format:178:7: note: declared here
178 | formatter() = delete; // No std::formatter specialization for this type.
| ^~~~~~~~~
/usr/include/c++/14.2.1/format:4095:38: note: use '-fdiagnostics-all-candidates' to display considered candidates
4095 | formatter<_Tp, _CharT> __f;
| ^~~
/usr/include/c++/14.2.1/format:4096:42: error: 'class std::formatter<std::complex<float>, char>' has no member named 'parse'
4096 | this->_M_pc.advance_to(__f.parse(this->_M_pc));
| ~~~~^~~~~
mps.hpp: In instantiation of 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]':
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
mps.hpp:162:19: error: call to consteval function 'std::basic_format_string<char, std::complex<float>&>("{} ")' is not a constant expression
162 | std::print("{} ", theta2.data()[i]);
| ~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/eigen3/Eigen/Core:175:
/usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h: In instantiation of 'constexpr void Eigen::internal::call_assignment_no_alias(Dst&, const Src&, const Func&) [with Dst = Eigen::Matrix<double, -1, 1>; Src = Eigen::Matrix<float, -1, 1>; Func = assign_op<double, float>]':
/usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:770:39: required from 'constexpr Derived& Eigen::PlainObjectBase<Derived>::_set_noalias(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Matrix<float, -1, 1>; Derived = Eigen::Matrix<double, -1, 1>]'
770 | internal::call_assignment_no_alias(this->derived(), other.derived(),
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
771 | internal::assign_op<Scalar, typename OtherDerived::Scalar>());
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:565:17: required from 'Eigen::PlainObjectBase<Derived>::PlainObjectBase(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Matrix<float, -1, 1>; Derived = Eigen::Matrix<double, -1, 1>]'
565 | _set_noalias(other);
| ~~~~~~~~~~~~^~~~~~~
/usr/include/eigen3/Eigen/src/Core/Matrix.h:386:108: required from 'Eigen::Matrix<Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_>::Matrix(const Eigen::EigenBase<OtherDerived>&) [with OtherDerived = Eigen::Matrix<float, -1, 1>; Scalar_ = double; int Rows_ = -1; int Cols_ = 1; int Options_ = 0; int MaxRows_ = -1; int MaxCols_ = 1]'
386 | EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Matrix(const EigenBase<OtherDerived>& other) : Base(other.derived()) {}
| ^
mps.hpp:136:24: required from 'std::tuple<Eigen::Tensor<dtype, 3, Storage, long int>, Eigen::Tensor<dtype, 1, Storage, long int>, Eigen::Tensor<dtype, 3, Storage, long int> > split_truncate_theta(const Eigen::Tensor<dtype, 4, Storage>&, size_t, dtype) [with dtype = std::complex<float>; unsigned int Storage = 1; size_t = long unsigned int]'
136 | const e::VectorXd& S = svd.singularValues();
| ^
main.cpp:37:62: required from here
37 | auto [A, S, B] = split_truncate_theta<dtype, e::RowMajor>(t4, chi_max, eps);
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h:839:3: error: static assertion failed: YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY
839 | EIGEN_CHECK_BINARY_COMPATIBILIY(Func, typename ActualDstTypeCleaned::Scalar, typename Src::Scalar);
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h:839:3: note: 'Eigen::internal::has_ReturnType<Eigen::ScalarBinaryOpTraits<double, float, Eigen::internal::assign_op<double, float> > >::value' evaluates to false
make: *** [Makefile:47: ../build/main.o] Error 1

390
src/run.cpp
View File

@ -1,390 +0,0 @@
#include <memory>
#include <type_traits>
#ifdef EMSCRIPTEN
#include <cmath>
#include <print>
// #include <array>
#include "mps.hpp"
#include "transverse_field_ising.hpp"
#include "bose_hubbard.hpp"
#include "tebd.hpp"
#include "util.hpp"
#include "type.hpp"
#include <functional>
#include <stdexcept>
#include <emscripten.h>
#include <emscripten/bind.h>
namespace ems = emscripten;
using TfiModel = TransverseFieldIsingModel1D<dtype>;
using BhModel = BoseHubbardModel1D<dtype>;
/**
* @brief A class for generating a range of values.
* @details
* Returns `nSteps` values with `start` being the first and `stop` being the last value.
* The same value is returned `countsPerStep` times before the next value is returned.
* After the last `stop` value has been returned, the cycle is restarted from the first `start` value.
*/
template<typename T>
struct Range {
using value_type = T;
Range(T start, T stop, unsigned nSteps, unsigned countsPerStep=1) : start(start), stop(stop), nSteps(nSteps), stepCounter(0), countsPerStep(countsPerStep), counter(0) {
if (nSteps <= 1) {
throw std::invalid_argument("nSteps must be greater than 1");
}
step = (stop - start) / (nSteps - 1);
};
T getNext() {
T result = start + stepCounter * step;
if (++counter >= countsPerStep) {
counter = 0;
if (++stepCounter >= nSteps) {
stepCounter = 0;
}
}
return result;
};
private:
T start;
T stop;
T step;
unsigned nSteps;
unsigned stepCounter;
unsigned countsPerStep;
unsigned counter;
};
/**
* This ugly mechanism lets the javascript pass
* any combination of Range and basic types to a function.
* Every time the function is called, the Range classes will return a new value while the rest stays the same.
*/
template<class T, template<class> class U>
inline constexpr bool is_instance_of_v = std::false_type{};
template<template<class> class U, class V>
inline constexpr bool is_instance_of_v<U<V>,U> = std::true_type{};
template<typename T>
concept isRange = is_instance_of_v<std::remove_reference_t<T>, Range>;
// All types except Range
template <typename T, typename Enable = void>
struct transform_type {
using type = T;
};
// Range
template <typename T>
struct transform_type<T, std::enable_if_t<isRange<T>>> {
using type = typename std::remove_reference_t<T>::value_type; // If T is a Range, return its value_type
};
static_assert(isRange<Range<int>>);
static_assert(isRange<Range<double>&>);
static_assert(!isRange<double>);
// Dispatcher
template<isRange T>
std::remove_reference_t<T>::value_type processType(T&& obj) {
return obj.getNext();
}
template<typename T>
requires(!isRange<T>)
T processType(T obj) {
return obj;
}
/// Make a shared model pointer from a arguments, which may contain Range objects
template<typename Model, typename... Args>
std::shared_ptr<Model> makeModel(Args&&... args) {
std::tuple<typename transform_type<Args>::type...> results = std::make_tuple(processType(std::forward<Args>(args))...);
return std::apply([](auto&&... args) { return std::make_shared<Model>(args...); }, results);
}
// template<typename Arg>
// void printArgs(Arg&& val) {
// std::cout << val << std::endl;
// }
// template<typename Arg, typename... Args>
// void printArgs(Arg&& val, Args&&... values) {
// std::cout << val << ", ";
// printArgs(std::forward<Args>(values)...);
// }
// template<typename... Args>
// void processDebug(Args&&... args) {
// std::tuple<typename transform_type<Args>::type...> results = std::make_tuple(processType(std::forward<Args>(args))...);
// std::apply([](auto&&... args){ printArgs(args...); }, results);
// }
/**
*/
template<typename Model>
class Results {
public:
virtual ~Results() {};
virtual postUpdateFunction_t<dtype> getPostUpdateFunction(std::shared_ptr<Model> model) = 0;
virtual void postUpdateFunction(std::shared_ptr<Model> model, const MPS<dtype>& psi) = 0;
};
class TfiResults : public Results<TfiModel> {
public:
using Model_t = TfiModel;
TfiResults(ems::val callback) : callback(callback) {
sigmaZ.setValues({{1, 0}, {0, -1}});
};
/// @todo write a second postUpdateFunction without the JS callback and return that one if callback is null
postUpdateFunction_t<dtype> getPostUpdateFunction(std::shared_ptr<TfiModel> model) override {
return [model, this](const MPS<dtype>& psi, unsigned i){ this->postUpdateFunction(model, psi); };
}
void postUpdateFunction(std::shared_ptr<TfiModel> model, const MPS<dtype>& psi) override {
auto Mval = psi.getSiteExpectationValueSum(this->sigmaZ).real();
auto Eval = model->energy(psi).real();
auto Sval = psi.getEntanglementEntropy(psi.getL()/2).real();
callback(Mval, Eval, Sval);
// hand control back to browser, otherwise site hangs and callback will not be executed
emscripten_sleep(0);
}
private:
e::TensorFixedSize<dtype, e::Sizes<2, 2>> sigmaZ;
ems::val callback;
};
/**
* @brief Measures <N>, E and S
*/
class BoseHubbardResults : public Results<BhModel> {
public:
BoseHubbardResults(ems::val callback) : callback(callback) {};
postUpdateFunction_t<dtype> getPostUpdateFunction(std::shared_ptr<BhModel> model) override {
return [this, model](const MPS<dtype>& psi, unsigned) { this->postUpdateFunction(model, psi); };
}
void postUpdateFunction(std::shared_ptr<BhModel> model, const MPS<dtype>& psi) override {
e::TensorMap<const e::Tensor<dtype, 2>> particleNumberOperator(model->getLadderOperators().getN().data(), model->localDim, model->localDim);
auto Nval = psi.getSiteExpectationValueSum(particleNumberOperator).real() / model->L;
auto Eval = model->energy(psi).real();
auto Sval = psi.getEntanglementEntropy(model->L/2).real();
callback(Nval, Eval, Sval);
// hand control back to browser, otherwise site hangs and callback will not be executed
emscripten_sleep(0);
}
private:
ems::val callback;
};
// EMSCRIPTEN_BINDINGS(my_class_example) {
// class_<MyClass>("MyClass")
// .constructor<int, std::string>()
// .function("incrementX", &MyClass::incrementX)
// .property("x", &MyClass::getX, &MyClass::setX)
// .property("x_readonly", &MyClass::getX)
// .class_function("getStringFromInstance", &MyClass::getStringFromInstance)
// ;
// };
enum Tfi_INIT {
INIT_SPIN_UP = 0,
};
// struct Tfi {
// Tfi(unsigned L, double _g, double _J, Tfi_INIT init) {
// throwIfLessThan(L, 2u, "L");
// dtype g(_g);
// dtype J(_J);
// throwIfNaN(g, "g");
// throwIfNaN(J, "J");
// this->model = TfiModel<dtype>(L, J, g);
// if (init == INIT_SPIN_UP) {
// mps = initMPSspinUp<dtype>(L);
// }
// }
// TfiModel<dtype> model;
// MPS<dtype> mps;
// };
// TODO unify
// Model1D<dtype> initModel_Tfi(unsigned L, double _g, double _J) {
// dtype g(_g);
// dtype J(_J);
// throwIfNaN(g, "g");
// throwIfNaN(J, "J");
// throwIfLessThan(L, 3u, "L");
// return TfiModel<dtype>(L, J, g);
// }
enum class TEBD_TYPE {
BRICKWALL = 0,
SECOND_ORDER = 1,
};
// EMSCRIPTEN_BINDINGS(tebd_module) {
// enum_<TEBD_TYPE>(TEBD_BRICKWALL)
// .value("TEBD_BRICKWALL", TEBD_BRICKWALL)
// ;
// class<MPS<dtype>>("MPS")
// .constructor<unsigned>()
// .function("getL", &MPS<dtype>::getL)
// .function("getN", &MPS<dtype>::getN)
// .function("getN_real",
// };
/**
* @brief Run TEBD on a model and measure energy, magnetization and half-chain entropy after each TEBD iteration
* @param callback: A javascript function taking 3 parameters, M, E and S after each iteration
*/
template<typename Model>
void runTEBD(MPS<dtype>& psi, const std::shared_ptr<Model> model, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, Results<Model>& result) {
dtype eps(_eps);
dtype dt(dtReal, dtImag);
// Sanity checks
throwIfNaN(eps, "eps");
throwIfNaN(dt, "dt");
throwIfLessThan(nSteps, 1u, "nSteps");
throwIfLessThan(chiMax, 1u, "chiMax");
if (psi.getL() != model->L) {
throw std::runtime_error("psi and model have different L");
}
std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast<int>(type), nSteps, dt, eps, chiMax);
// collect M over time in an array
auto F = result.getPostUpdateFunction(model);
// initial values now
F(psi, 0);
// std::println("Bonds:");
const auto& bondsMat = model->getH_BondsMatrices();
// std::cout << "First: \n" << *(bondsMat.front()) << std::endl;
// std::cout << "Secnd: \n" << *(bondsMat.at(1)) << std::endl;
// std::cout << "Last : \n" << *(bondsMat.back()) << std::endl;
// const Bonds<dtype> bonds = model.getH_Bonds();
const Bonds<dtype> bonds = getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim);
// std::println("U:");
// int i = 0;
// for (const auto& bond : bonds) {
// printTensor(*bond, std::format("Bond[{:02}]", i));
// i++;
// }
if (type == TEBD_TYPE::BRICKWALL) {
runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, F);
} else if (type == TEBD_TYPE::SECOND_ORDER) {
const Bonds<dtype> bondsHalf = getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt/2., model->localDim);
runTEBD_secondOrder(psi, bonds, bondsHalf, chiMax, eps, nSteps, F);
} else {
throw std::invalid_argument(std::format("Got invalid TEBD_TYPE: {}", static_cast<int>(type)));
}
std::println("Chis: {}", psi.getChi());
std::println("Collapse: {}", psi.collapse(0));
};
// @TODO
template<typename Model, typename... ModelArgs>
void runTEBD_phaseDiagram(MPS<dtype>& psiI, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, unsigned nValues, Results<Model>& result, ModelArgs&&... args) {
dtype eps(_eps);
dtype dt(dtReal, dtImag);
size_t L = psiI.getL();
// Sanity checks
throwIfNaN(eps, "eps");
throwIfNaN(dt, "dt");
throwIfLessThan(nSteps, 2u, "nSteps");
throwIfLessThan(chiMax, 1u, "chiMax");
std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast<int>(type), nSteps, dt, eps, chiMax);
postUpdateFunction_t<dtype> F = [](const MPS<dtype>&, unsigned) {
emscripten_sleep(0);
};
for (unsigned i = 0; i < nValues; i++) {
// TfiModel model(L, U, g);
std::shared_ptr<Model> model = makeModel<Model>(std::forward<ModelArgs>(args)...);
std::println("{:02} - Model: {}", i, *model);
auto psi = psiI;
const Bonds<dtype> bonds = getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim);
runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, F);
result.postUpdateFunction(model, psi);
}
}
void test(unsigned x) {
Range<double> range(0, 1.0, 3, 2);
std::println("test({})", x);
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("{}", processType(0.5));
std::println("{}", processType(range));
// processDebug(range, 4.3, -3, range);
// process2(range);
}
EMSCRIPTEN_BINDINGS(run) {
// see https://emscripten.org/docs/porting/connecting_cpp_and_javascript/embind.html#advanced-class-concepts
// Results
ems::class_<Results<TfiModel>>("ResultsBaseTfi");
ems::class_<Results<BhModel>>("ResultsBaseBh");
ems::class_<TfiResults, ems::base<Results<TfiModel>>>("ResultsTfi")
.constructor<ems::val>()
;
ems::class_<BoseHubbardResults, ems::base<Results<BhModel>>>("ResultsBh")
.constructor<ems::val>()
;
// Models
ems::class_<Model1D<dtype>>("Model1D")
.smart_ptr<std::shared_ptr<Model1D<dtype>>>("Model1D")
;
ems::class_<TfiModel, ems::base<Model1D<dtype>>>("ModelTfi")
// .constructor<unsigned, double ,double>()
.smart_ptr_constructor("ModelTfi", &std::make_shared<TfiModel, unsigned, double ,double>)
;
ems::class_<BhModel, ems::base<Model1D<dtype>>>("ModelBh")
// .constructor<unsigned, double, double, double, unsigned>()
.smart_ptr_constructor("ModelBh", &std::make_shared<BhModel, unsigned, double, double, double, unsigned>)
;
ems::class_<Range<double>>("DoubleRange")
.constructor<double, double, unsigned, unsigned>()
;
// Other
ems::class_<MPS<dtype>>("MPS");
ems::function("initMpsSpinUp", &initMPS_spinUp<dtype>);
ems::function("initMpsSpinRight", &initMPS_spinRight<dtype>);
ems::function("initMpsNParticles", &initMPS_nParticles<dtype>);
// ems::function("initMPS_spinUp", &initMPS_spinUp);
ems::register_vector<double>("stdVectorDouble");
ems::register_vector<std::vector<double>>("stdVectorDouble2D");
// ems::class_<Results>("TEBDResults")
// .property("M", &Results::M)
// .property("E", &Results::E)
// .property("S", &Results::S)
// ;
ems::enum_<TEBD_TYPE>("TEBD_TYPE")
.value("BRICKWALL", TEBD_TYPE::BRICKWALL)
.value("SECOND_ORDER", TEBD_TYPE::SECOND_ORDER)
;
ems::function("runTebdTfi", &runTEBD<TfiModel>);
ems::function("runTebdPhaseDiagramTfi", &runTEBD_phaseDiagram<TfiModel, unsigned, Range<double>, Range<double>>);
ems::function("runTebdBh", &runTEBD<BhModel>);
ems::function("runTebdPhaseDiagramBh", &runTEBD_phaseDiagram<BhModel, unsigned, double, Range<double>, Range<double>, unsigned>);
ems::function("test", &test);
};
#endif

2
src/type.hpp
View File

@ -2,3 +2,5 @@
/// the datatype for all template instantiations
using dtype = std::complex<double>;
// template<int N>
// using Tensor = std::shared_ptr<e::Tensor<dtype, N>>;

27
src/wasm/emscripten.hpp Normal file
View File

@ -0,0 +1,27 @@
#pragma once
#ifdef EMSCRIPTEN
#include <emscripten.h>
#include <emscripten/bind.h>
namespace ems = emscripten;
// #warning emscripten
#else
// #warning noEmscripten
#include <print>
/**
* @brief Allow using ems::val without emscripten
*/
namespace ems {
struct val {
template<typename... Args>
void operator()(Args...) {
std::println("Callback called");
}
};
}
inline void emscripten_sleep(int time) {
std::println("emscripten_sleep({})", time);
}
#endif

88
src/wasm/range.hpp Normal file
View File

@ -0,0 +1,88 @@
#pragma once
#include <memory>
#include <stdexcept>
/**
* @brief A class for generating a range of values.
* @details
* Returns `nSteps` values with `start` being the first and `stop` being the last value.
* The same value is returned `countsPerStep` times before the next value is returned.
* After the last `stop` value has been returned, the cycle is restarted from the first `start` value.
*/
template<typename T>
struct Range {
using value_type = T;
Range(T start, T stop, unsigned nSteps, unsigned countsPerStep=1) : start(start), stop(stop), nSteps(nSteps), stepCounter(0), countsPerStep(countsPerStep), counter(0) {
if (nSteps <= 1) {
throw std::invalid_argument("nSteps must be greater than 1");
}
step = (stop - start) / (nSteps - 1);
};
T getNext() {
T result = start + stepCounter * step;
if (++counter >= countsPerStep) {
counter = 0;
if (++stepCounter >= nSteps) {
stepCounter = 0;
}
}
return result;
};
private:
T start;
T stop;
T step;
unsigned nSteps;
unsigned stepCounter;
unsigned countsPerStep;
unsigned counter;
};
/**
* This ugly mechanism lets the javascript pass
* any combination of Range and basic types to a function.
* Every time the function is called, the Range classes will return a new value while the rest stays the same.
*/
template<class T, template<class> class U>
inline constexpr bool is_instance_of_v = std::false_type{};
template<template<class> class U, class V>
inline constexpr bool is_instance_of_v<U<V>,U> = std::true_type{};
template<typename T>
concept isRange = is_instance_of_v<std::remove_reference_t<T>, Range>;
// All types except Range
template <typename T, typename Enable = void>
struct transform_type {
using type = T;
};
// Range
template <typename T>
struct transform_type<T, std::enable_if_t<isRange<T>>> {
using type = typename std::remove_reference_t<T>::value_type; // If T is a Range, return its value_type
};
static_assert(isRange<Range<int>>);
static_assert(isRange<Range<double>&>);
static_assert(!isRange<double>);
// Dispatcher
template<isRange T>
std::remove_reference_t<T>::value_type processType(T&& obj) {
return obj.getNext();
}
template<typename T>
requires(!isRange<T>)
T processType(T obj) {
return obj;
}
/// Make a shared model pointer from a arguments, which may contain Range objects
template<typename Model, typename... Args>
std::shared_ptr<Model> makeModel(Args&&... args) {
std::tuple<typename transform_type<Args>::type...> results = std::make_tuple(processType(std::forward<Args>(args))...);
return std::apply([](auto&&... args) { return std::make_shared<Model>(args...); }, results);
}

81
src/wasm/result.hpp Normal file
View File

@ -0,0 +1,81 @@
#pragma once
#include "algorithms/tebd.hpp"
#include "mps.hpp"
#include "models/transverse_field_ising.hpp"
#include "models/bose_hubbard.hpp"
#include "emscripten.hpp"
#include <cstdlib>
#include <iterator>
#include <limits>
#include <numeric>
#include <print>
using TfiModel = TransverseFieldIsingModel1D<dtype>;
using BhModel = BoseHubbardModel1D<dtype>;
template<typename Model>
class Results {
public:
virtual ~Results() {};
virtual postUpdateFunction_t<dtype> getPostUpdateFunction(std::shared_ptr<Model> model) = 0;
virtual bool postUpdateFunction(std::shared_ptr<Model> model, const MPS<dtype>& psi) = 0;
};
class TfiResults : public Results<TfiModel> {
public:
using Model_t = TfiModel;
TfiResults(ems::val callback, double convergence) : callback(callback), convergence(convergence), previousE(std::numeric_limits<double>::max()) {
sigmaZ.setValues({{1, 0}, {0, -1}});
};
/// @todo write a second postUpdateFunction without the JS callback and return that one if callback is null
postUpdateFunction_t<dtype> getPostUpdateFunction(std::shared_ptr<TfiModel> model) override {
return [this, model](const MPS<dtype>& psi, unsigned i){ return this->postUpdateFunction(model, psi); };
}
bool postUpdateFunction(std::shared_ptr<TfiModel> model, const MPS<dtype>& psi) override {
auto Mval = psi.getSiteExpectationValueSum(this->sigmaZ).real() / psi.getL();
auto Eval = model->energy(psi).real();
auto Sval = psi.getEntanglementEntropy(psi.getL()/2).real();
auto chiVal = psi.getChiMax();
callback(Mval, Eval, Sval, chiVal);
// hand control back to browser, otherwise site hangs and callback will not be executed
emscripten_sleep(0);
if (std::abs(Eval - previousE) < convergence) return true;
previousE = Eval;
return false;
}
private:
e::TensorFixedSize<dtype, e::Sizes<2, 2>> sigmaZ;
ems::val callback;
double convergence;
double previousE;
};
/**
* @brief Measures <N>, E and S
*/
class BoseHubbardResults : public Results<BhModel> {
public:
BoseHubbardResults(ems::val callback, double convergence) : callback(callback), convergence(convergence), previousE(std::numeric_limits<double>::max()){};
postUpdateFunction_t<dtype> getPostUpdateFunction(std::shared_ptr<BhModel> model) override {
return [this, model](const MPS<dtype>& psi, unsigned) { return this->postUpdateFunction(model, psi); };
}
bool postUpdateFunction(std::shared_ptr<BhModel> model, const MPS<dtype>& psi) override {
e::TensorMap<const e::Tensor<dtype, 2>> particleNumberOperator(model->getLadderOperators().getN().data(), model->localDim, model->localDim);
auto Nval = psi.getSiteExpectationValueSum(particleNumberOperator).real() / model->L;
auto Eval = model->energy(psi).real();
auto Sval = psi.getEntanglementEntropy(model->L/2).real();
auto chiVal = psi.getChiMax();
callback(Nval, Eval, Sval, chiVal);
// hand control back to browser, otherwise site hangs and callback will not be executed
emscripten_sleep(0);
if (std::abs(Eval - previousE) < convergence) return true;
previousE = Eval;
return false;
}
private:
ems::val callback;
double convergence;
double previousE;
};

107
src/wasm/run.cpp Normal file
View File

@ -0,0 +1,107 @@
#include "wasm/run.hpp"
#ifndef BUILD_OPTIONS_STRING
#define BUILD_OPTIONS_STRING "unknown"
#endif
std::string getBuildOptions() {
return BUILD_OPTIONS_STRING;
}
void test(unsigned x) {
Range<double> range(0, 1.0, 3, 2);
std::println("test({})", x);
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("range({})", range.getNext());
std::println("{}", processType(0.5));
std::println("{}", processType(range));
// processDebug(range, 4.3, -3, range);
// process2(range);
}
/**
* @brief Function used to test performance in the browser
*/
int testPerformance() {
e::Tensor<dtype, 4> t4(4, 16, 16, 4);
auto data = t4.data();
for (unsigned i = 0; i < t4.size(); i++) {
data[i] = i;
}
auto [A, S, B] = split_truncate_theta<dtype>(t4, 20, 1e-10);
e::Tensor<dtype, 4> t4reconstructedButNormalized = A.contract(diagonal(S).contract(B, e::array<iP, 1>{iP{1, 0}}), e::array<iP, 1>{iP{2, 0}});
return t4reconstructedButNormalized.size();
}
#ifdef EMSCRIPTEN
EMSCRIPTEN_BINDINGS(run) {
// see https://emscripten.org/docs/porting/connecting_cpp_and_javascript/embind.html#advanced-class-concepts
// RESULTS
ems::class_<Results<TfiModel>>("ResultsBaseTfi");
ems::class_<Results<BhModel>>("ResultsBaseBh");
ems::class_<TfiResults, ems::base<Results<TfiModel>>>("ResultsTfi")
.constructor<ems::val, double>()
;
ems::class_<BoseHubbardResults, ems::base<Results<BhModel>>>("ResultsBh")
.constructor<ems::val, double>()
;
// MODELS
ems::class_<Model1D<dtype>>("Model1D")
.smart_ptr<std::shared_ptr<Model1D<dtype>>>("Model1D")
;
ems::class_<TfiModel, ems::base<Model1D<dtype>>>("ModelTfi")
// .constructor<unsigned, double ,double>()
.smart_ptr_constructor("ModelTfi", &std::make_shared<TfiModel, unsigned, double ,double>)
;
ems::class_<BhModel, ems::base<Model1D<dtype>>>("ModelBh")
// .constructor<unsigned, double, double, double, unsigned>()
.smart_ptr_constructor("ModelBh", &std::make_shared<BhModel, unsigned, double, double, double, unsigned>)
;
ems::class_<Range<double>>("DoubleRange")
.constructor<double, double, unsigned, unsigned>()
;
// MPS
ems::class_<MPS<dtype>>("MPS");
ems::function("initMpsSpinUp", &initMPS_spinUp<dtype>);
ems::function("initMpsSpinRight", &initMPS_spinRight<dtype>);
ems::function("initMpsNParticles", &initMPS_nParticles<dtype>);
// ems::function("initMPS_spinUp", &initMPS_spinUp);
ems::register_vector<double>("stdVectorDouble");
ems::register_vector<std::vector<double>>("stdVectorDouble2D");
// ems::class_<Results>("TEBDResults")
// .property("M", &Results::M)
// .property("E", &Results::E)
// .property("S", &Results::S)
// ;
// TEBD
ems::enum_<TEBD_TYPE>("TEBD_TYPE")
.value("BRICKWALL", TEBD_TYPE::BRICKWALL)
.value("SECOND_ORDER", TEBD_TYPE::SECOND_ORDER)
;
ems::function("runTebdTfi", &runTEBD<TfiModel>);
ems::function("runTebdPhaseDiagramTfi", &runTEBD_phaseDiagram<TfiModel, unsigned, Range<double>, Range<double>>);
ems::function("runTebdBh", &runTEBD<BhModel>);
ems::function("runTebdPhaseDiagramBh", &runTEBD_phaseDiagram<BhModel, unsigned, double, Range<double>, Range<double>, unsigned>);
// DMRG
ems::function("runDmrgTfi", &runDMRG<TfiModel>);
ems::function("runDmrgPhaseDiagramTfi", &runDMRG_phaseDiagram<TfiModel, unsigned, Range<double>, Range<double>>);
ems::function("runDmrgBh", &runDMRG<BhModel>);
ems::function("runDmrgPhaseDiagramBh", &runDMRG_phaseDiagram<BhModel, unsigned, double, Range<double>, Range<double>, unsigned>);
// OTHER
ems::function("test", &test);
ems::function("getBuildOptions", &getBuildOptions);
ems::function("testPerformance", &testPerformance);
};
#endif

166
src/wasm/run.hpp Normal file
View File

@ -0,0 +1,166 @@
#include <memory>
#include <print>
#include "wasm/result.hpp"
#include "algorithms/tebd.hpp"
#include "algorithms/dmrg.hpp"
#include "util.hpp"
#include "type.hpp"
#include "wasm/emscripten.hpp"
#include "wasm/range.hpp"
// template<typename Arg>
// void printArgs(Arg&& val) {
// std::cout << val << std::endl;
// }
// template<typename Arg, typename... Args>
// void printArgs(Arg&& val, Args&&... values) {
// std::cout << val << ", ";
// printArgs(std::forward<Args>(values)...);
// }
// template<typename... Args>
// void processDebug(Args&&... args) {
// std::tuple<typename transform_type<Args>::type...> results = std::make_tuple(processType(std::forward<Args>(args))...);
// std::apply([](auto&&... args){ printArgs(args...); }, results);
// }
/**
*/
// EMSCRIPTEN_BINDINGS(my_class_example) {
// class_<MyClass>("MyClass")
// .constructor<int, std::string>()
// .function("incrementX", &MyClass::incrementX)
// .property("x", &MyClass::getX, &MyClass::setX)
// .property("x_readonly", &MyClass::getX)
// .class_function("getStringFromInstance", &MyClass::getStringFromInstance)
// ;
// };
/**
* @brief Return a string that was set at compilation using the BUILD_OPTIONS_STRING variable
*/
std::string getBuildOptions();
enum class TEBD_TYPE {
BRICKWALL = 0,
SECOND_ORDER = 1,
};
// EMSCRIPTEN_BINDINGS(tebd_module) {
// enum_<TEBD_TYPE>(TEBD_BRICKWALL)
// .value("TEBD_BRICKWALL", TEBD_BRICKWALL)
// ;
// class<MPS<dtype>>("MPS")
// .constructor<unsigned>()
// .function("getL", &MPS<dtype>::getL)
// .function("getN", &MPS<dtype>::getN)
// .function("getN_real",
// };
// TEBD
template<typename Model>
void runTEBD(MPS<dtype>& psi, const std::shared_ptr<Model> model, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, Results<Model>& result) {
dtype eps(_eps);
dtype dt(dtReal, dtImag);
// sanity checks
throwIfNaN(eps, "eps");
throwIfNaN(dt, "dt");
throwIfLessThan(nSteps, 1u, "nSteps");
throwIfLessThan(chiMax, 1u, "chiMax");
if (psi.getL() != model->L) {
throw std::runtime_error("psi and model have different L");
}
std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast<int>(type), nSteps, dt, eps, chiMax);
auto F = result.getPostUpdateFunction(model);
// initial values now
F(psi, 0);
const Bonds<dtype> bonds = tebd::getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim);
if (type == TEBD_TYPE::BRICKWALL) {
tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, F);
} else if (type == TEBD_TYPE::SECOND_ORDER) {
const Bonds<dtype> bondsHalf = tebd::getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt/2., model->localDim);
tebd::runSecondOrder(psi, bonds, bondsHalf, chiMax, eps, nSteps, F);
} else {
throw std::invalid_argument(std::format("Got invalid TEBD_TYPE: {}", static_cast<int>(type)));
}
std::println("Chis: {}", psi.getChi());
std::println("Collapse: {}", psi.collapse(0));
};
template<typename Model, typename... ModelArgs>
void runTEBD_phaseDiagram(const MPS<dtype>& psiI, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, unsigned nValues, Results<Model>& result, ModelArgs&&... args) {
dtype eps(_eps);
dtype dt(dtReal, dtImag);
// sanity checks
throwIfNaN(eps, "eps");
throwIfNaN(dt, "dt");
throwIfLessThan(nSteps, 2u, "nSteps");
throwIfLessThan(chiMax, 1u, "chiMax");
std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast<int>(type), nSteps, dt, eps, chiMax);
postUpdateFunction_t<dtype> F = [](const MPS<dtype>&, unsigned) {
emscripten_sleep(0);
return false;
};
for (unsigned i = 0; i < nValues; i++) {
std::shared_ptr<Model> model = makeModel<Model>(std::forward<ModelArgs>(args)...);
std::println("{:02} - Model: {}", i, *model);
MPS<dtype> psi(psiI);
std::println("New psi: value at B[2](0, 1, 0)={}", psi.B(2)(0, 1, 0));
const Bonds<dtype> bonds = tebd::getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim);
tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, F);
result.postUpdateFunction(model, psi);
}
}
// DMRG
template<typename Model>
void runDMRG(MPS<dtype>& psi, const std::shared_ptr<Model> model, unsigned nSweeps, double _eps, unsigned chiMax, Results<Model>& result) {
dtype eps(_eps);
// sanity checks
throwIfNaN(eps, "eps");
throwIfLessThan(nSweeps, 1u, "nSweeps");
throwIfLessThan(chiMax, 1u, "chiMax");
if (psi.getL() != model->L) {
throw std::runtime_error("psi and model have different L");
}
std::println("Running DMRG nSweeps={}, eps={}, chiMax={}", nSweeps, eps, chiMax);
auto F = result.getPostUpdateFunction(model);
// initial values now
F(psi, 0);
dmrg::run(psi, model->getH_MPO(), chiMax, eps, nSweeps, F);
std::println("Chis: {}", psi.getChi());
std::println("Collapse: {}", psi.collapse(0));
};
template<typename Model, typename... ModelArgs>
void runDMRG_phaseDiagram(const MPS<dtype>& psiI, unsigned nSweeps, double _eps, unsigned chiMax, unsigned nValues, Results<Model>& result, ModelArgs&&... args) {
dtype eps(_eps);
// sanity checks
throwIfNaN(eps, "eps");
throwIfLessThan(nSweeps, 1u, "nSweeps");
throwIfLessThan(chiMax, 1u, "chiMax");
std::println("Running DMRG nSweeps={}, eps={}, chiMax={}", nSweeps, eps, chiMax);
postUpdateFunction_t<dtype> F = [](const MPS<dtype>&, unsigned) {
emscripten_sleep(0);
return false;
};
for (unsigned i = 0; i < nValues; i++) {
std::shared_ptr<Model> model = makeModel<Model>(std::forward<ModelArgs>(args)...);
std::println("{:02} - Model: {}", i, *model);
MPS<dtype> psi(psiI);
dmrg::run(psi, model->getH_MPO(), chiMax, eps, nSweeps, F);
result.postUpdateFunction(model, psi);
}
}

0
src/wasm/util.hpp Normal file
View File