From 2f77b542bdb968908e909e1eefc47a0637cbd796 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Sat, 12 Oct 2024 11:26:39 +0200
Subject: [PATCH 01/16] Add add_mps for L=1

---
 src/mps/mps.c | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++
 src/mps/mps.h |  4 ++-
 2 files changed, 72 insertions(+), 1 deletion(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index b5ad46e..d38294b 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -316,6 +316,75 @@ void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret)
 }
 
 
+//________________________________________________________________________________________________________________________
+///
+/// \brief Compute the addition of two MPS chi and psi.
+///
+void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
+{
+	// Number of lattice sites must agree.
+	assert(chi->nsites == psi->nsites);
+
+	// Number of lattice sites must be larger than 0.
+	assert(chi->nsites > 0);
+
+	// Physical quantum numbers must agree.
+	assert(chi->d == psi->d); 
+	assert(qnumber_all_equal(chi->d, chi->qsite, psi->qsite)); 
+
+	long d = chi->d;
+	int L = chi->nsites;
+	enum numeric_type dtype = chi->a[0].dtype;
+
+	if (L == 1) {
+		const struct block_sparse_tensor* chi_a = &(chi->a[0]);
+		const struct block_sparse_tensor* psi_a = &(psi->a[0]);
+		const int ndim = chi_a->ndim;
+
+        // Dummy bond quantum numbers must agree.
+		assert(
+			qnumber_all_equal(
+				ndim, 
+				chi_a->qnums_logical[0], 
+				psi_a->qnums_logical[0]
+			)
+		); 
+		assert(
+			qnumber_all_equal(
+				ndim,
+				chi_a->qnums_logical[1],
+				psi_a->qnums_logical[1]
+			)
+		); 
+
+		// Add MPS tensors. 
+		const struct block_sparse_tensor* new_a;
+		copy_block_sparse_tensor(chi_a, new_a);
+
+		// Add individual dense_tensors.
+		for (int i = 0; i < new_a->ndim; i++) {
+			dense_tensor_scalar_multiply_add(
+				1,
+				psi_a->blocks[i],
+				new_a->blocks[i]
+			);
+		}
+
+		// Allocate memory and set return value 'ret'.
+		struct mps ret_val = {
+			.a=new_a,
+			.d=d,
+			.nsites=L,
+			.qsite=chi->qsite
+		};
+		ret = ct_malloc(sizeof(struct mps));
+		memcpy((void*) &ret_val, (void*) ret, sizeof(ret));
+	} else {
+
+	}
+}
+
+
 //________________________________________________________________________________________________________________________
 ///
 /// \brief Compute the Euclidean norm of the MPS.
diff --git a/src/mps/mps.h b/src/mps/mps.h
index ba5ef8d..e34f6a2 100644
--- a/src/mps/mps.h
+++ b/src/mps/mps.h
@@ -49,10 +49,12 @@ static inline long mps_bond_dim(const struct mps* mps, const int i)
 //________________________________________________________________________________________________________________________
 //
 
-// inner product and norm
+// inner product, addition, and norm
 
 void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret);
 
+void mps_add(const struct mps* chi, const struct mps* psi, void* ret);
+
 double mps_norm(const struct mps* psi);
 
 

From 894346f5e0d0127de43afe8bd67d2b704d171c16 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Sat, 12 Oct 2024 14:28:45 +0200
Subject: [PATCH 02/16] Add simple unit test

---
 src/mps/mps.c       |  40 ++++++------
 test/mps/test_mps.c |  47 ++++++++++++++
 test/run_tests.c    | 152 ++++++++++++++++++++++----------------------
 3 files changed, 142 insertions(+), 97 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index d38294b..5fad6a5 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -335,6 +335,10 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 	long d = chi->d;
 	int L = chi->nsites;
 	enum numeric_type dtype = chi->a[0].dtype;
+	qnumber *qsite = ct_malloc(d * sizeof(qnumber));
+	for (int i = 0; i < d; i++) {
+		qsite[i] = chi->qsite[i];
+	}
 
 	if (L == 1) {
 		const struct block_sparse_tensor* chi_a = &(chi->a[0]);
@@ -346,39 +350,31 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			qnumber_all_equal(
 				ndim, 
 				chi_a->qnums_logical[0], 
-				psi_a->qnums_logical[0]
-			)
-		); 
+				psi_a->qnums_logical[0])); 
 		assert(
 			qnumber_all_equal(
 				ndim,
 				chi_a->qnums_logical[1],
-				psi_a->qnums_logical[1]
-			)
-		); 
+				psi_a->qnums_logical[1])); 
 
 		// Add MPS tensors. 
-		const struct block_sparse_tensor* new_a;
-		copy_block_sparse_tensor(chi_a, new_a);
+		struct block_sparse_tensor *ret_a = ct_malloc(sizeof(struct block_sparse_tensor));
+		copy_block_sparse_tensor(chi_a, ret_a);
 
 		// Add individual dense_tensors.
-		for (int i = 0; i < new_a->ndim; i++) {
-			dense_tensor_scalar_multiply_add(
-				1,
-				psi_a->blocks[i],
-				new_a->blocks[i]
-			);
+		for (int i = 0; i < ret_a->ndim; i++) {
+			if (psi_a->blocks[i] != NULL && ret_a->blocks[i] != NULL) {
+				dense_tensor_scalar_multiply_add(
+					numeric_one(psi_a->dtype),
+					psi_a->blocks[i],
+					ret_a->blocks[i]);
+			}
 		}
 
 		// Allocate memory and set return value 'ret'.
-		struct mps ret_val = {
-			.a=new_a,
-			.d=d,
-			.nsites=L,
-			.qsite=chi->qsite
-		};
-		ret = ct_malloc(sizeof(struct mps));
-		memcpy((void*) &ret_val, (void*) ret, sizeof(ret));
+		// ret = ct_malloc(sizeof(struct mps));
+		struct mps VALUE = {.nsites=L, .d=d, .qsite=qsite, .a=ret_a};
+		memcpy((void*)ret, (void*)&VALUE, sizeof(struct mps));
 	} else {
 
 	}
diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index 516a123..aa6f4ff 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -642,3 +642,50 @@ char* test_mps_to_statevector()
 
 	return 0;
 }
+
+char* test_mps_add()
+{
+	struct mps chi, psi, res;
+	struct block_sparse_tensor chi_vec, psi_vec, res_vec;
+	
+	const int nsites = 1;
+	const long d = 2, max_vdim = 16;
+	
+	struct rng_state rng_state;
+	seed_rng_state(42, &rng_state);
+
+	const qnumber q_pnum = 7;
+	const qnumber q_spin = 1;
+	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
+
+	qnumber *qsite = ct_malloc(d * sizeof(qnumber));
+	for (int i = 0; i < d; i++) {
+		qsite[i] = qnum_sector;
+	}
+	
+	construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &chi);
+	construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &psi);
+
+	mps_add(&chi, &psi, &res);
+
+	mps_to_statevector(&chi, &chi_vec);
+	mps_to_statevector(&psi, &psi_vec);
+	mps_to_statevector(&res, &res_vec);
+
+	double *chi_data = (double*)(chi_vec.blocks[0]->data);
+	double *psi_data = (double*)(psi_vec.blocks[0]->data);
+	double *res_data = (double*)(res_vec.blocks[0]->data);
+	for (int i = 0; i < chi_vec.blocks[0]->ndim; i++){
+		assert(chi_data[i] + psi_data[i] == res_data[i]);
+	}
+
+	delete_block_sparse_tensor(&res_vec);
+	delete_block_sparse_tensor(&psi_vec);
+	delete_block_sparse_tensor(&chi_vec);
+	ct_free(qsite);
+	delete_mps(&res);
+	delete_mps(&psi);
+	delete_mps(&chi);
+
+	return 0;
+}
\ No newline at end of file
diff --git a/test/run_tests.c b/test/run_tests.c
index 62a6117..410b0c1 100644
--- a/test/run_tests.c
+++ b/test/run_tests.c
@@ -56,6 +56,7 @@ char* test_mps_orthonormalize_qr();
 char* test_mps_compress();
 char* test_mps_split_tensor_svd();
 char* test_mps_to_statevector();
+char* test_mps_add();
 char* test_queue();
 char* test_linked_list();
 char* test_hash_table();
@@ -94,81 +95,82 @@ char* test_operator_average_coefficient_gradient();
 int main()
 {
 	struct test tests[] = {
-		TEST_FUNCTION_ENTRY(test_dense_tensor_trace),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_cyclic_partial_trace),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_transpose),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_slice),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_pad_zeros),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_pointwise),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_axis),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_dot),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_dot_update),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product_degree_zero),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_concatenate),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_block_diag),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_qr),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_rq),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_svd),
-		TEST_FUNCTION_ENTRY(test_dense_tensor_block),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_copy),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_block),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_cyclic_partial_trace),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_norm2),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_transpose),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_reshape),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_slice),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_pointwise_vector),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_axis),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_dot),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_concatenate),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_block_diag),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_qr),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_rq),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_svd),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_serialize),
-		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_entry),
-		TEST_FUNCTION_ENTRY(test_clebsch_gordan_coefficients),
-		TEST_FUNCTION_ENTRY(test_su2_tree_enumerate_charge_sectors),
-		TEST_FUNCTION_ENTRY(test_su2_fuse_split_tree_enumerate_charge_sectors),
-		TEST_FUNCTION_ENTRY(test_su2_tensor_fmove),
-		TEST_FUNCTION_ENTRY(test_su2_tensor_contract_simple),
-		TEST_FUNCTION_ENTRY(test_su2_to_dense_tensor),
-		TEST_FUNCTION_ENTRY(test_mps_vdot),
-		TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr),
-		TEST_FUNCTION_ENTRY(test_mps_compress),
-		TEST_FUNCTION_ENTRY(test_mps_split_tensor_svd),
-		TEST_FUNCTION_ENTRY(test_mps_to_statevector),
-		TEST_FUNCTION_ENTRY(test_queue),
-		TEST_FUNCTION_ENTRY(test_linked_list),
-		TEST_FUNCTION_ENTRY(test_hash_table),
-		TEST_FUNCTION_ENTRY(test_bipartite_graph_maximum_cardinality_matching),
-		TEST_FUNCTION_ENTRY(test_bipartite_graph_minimum_vertex_cover),
-		TEST_FUNCTION_ENTRY(test_integer_hermite_normal_form),
-		TEST_FUNCTION_ENTRY(test_integer_backsubstitute),
-		TEST_FUNCTION_ENTRY(test_lanczos_iteration_d),
-		TEST_FUNCTION_ENTRY(test_lanczos_iteration_z),
-		TEST_FUNCTION_ENTRY(test_eigensystem_krylov_symmetric),
-		TEST_FUNCTION_ENTRY(test_eigensystem_krylov_hermitian),
-		TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_basic),
-		TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_advanced),
-		TEST_FUNCTION_ENTRY(test_mpo_from_assembly),
-		TEST_FUNCTION_ENTRY(test_ttno_graph_from_opchains),
-		TEST_FUNCTION_ENTRY(test_ttno_from_assembly),
-		TEST_FUNCTION_ENTRY(test_ising_1d_mpo),
-		TEST_FUNCTION_ENTRY(test_heisenberg_xxz_1d_mpo),
-		TEST_FUNCTION_ENTRY(test_bose_hubbard_1d_mpo),
-		TEST_FUNCTION_ENTRY(test_fermi_hubbard_1d_mpo),
-		TEST_FUNCTION_ENTRY(test_molecular_hamiltonian_mpo),
-		TEST_FUNCTION_ENTRY(test_spin_molecular_hamiltonian_mpo),
-		TEST_FUNCTION_ENTRY(test_retained_bond_indices),
-		TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd),
-		TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd_zero),
-		TEST_FUNCTION_ENTRY(test_operator_inner_product),
-		TEST_FUNCTION_ENTRY(test_apply_operator),
-		TEST_FUNCTION_ENTRY(test_dmrg_singlesite),
-		TEST_FUNCTION_ENTRY(test_dmrg_twosite),
-		TEST_FUNCTION_ENTRY(test_operator_average_coefficient_gradient),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_trace),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_cyclic_partial_trace),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_transpose),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_slice),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_pad_zeros),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_pointwise),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_axis),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_dot),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_dot_update),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product_degree_zero),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_concatenate),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_block_diag),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_qr),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_rq),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_svd),
+		// TEST_FUNCTION_ENTRY(test_dense_tensor_block),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_copy),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_block),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_cyclic_partial_trace),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_norm2),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_transpose),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_reshape),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_slice),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_pointwise_vector),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_axis),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_dot),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_concatenate),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_block_diag),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_qr),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_rq),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_svd),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_serialize),
+		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_entry),
+		// TEST_FUNCTION_ENTRY(test_clebsch_gordan_coefficients),
+		// TEST_FUNCTION_ENTRY(test_su2_tree_enumerate_charge_sectors),
+		// TEST_FUNCTION_ENTRY(test_su2_fuse_split_tree_enumerate_charge_sectors),
+		// TEST_FUNCTION_ENTRY(test_su2_tensor_fmove),
+		// TEST_FUNCTION_ENTRY(test_su2_tensor_contract_simple),
+		// TEST_FUNCTION_ENTRY(test_su2_to_dense_tensor),
+		// TEST_FUNCTION_ENTRY(test_mps_vdot),
+		// TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr),
+		// TEST_FUNCTION_ENTRY(test_mps_compress),
+		// TEST_FUNCTION_ENTRY(test_mps_split_tensor_svd),
+		// TEST_FUNCTION_ENTRY(test_mps_to_statevector),
+		TEST_FUNCTION_ENTRY(test_mps_add),
+		// TEST_FUNCTION_ENTRY(test_queue),
+		// TEST_FUNCTION_ENTRY(test_linked_list),
+		// TEST_FUNCTION_ENTRY(test_hash_table),
+		// TEST_FUNCTION_ENTRY(test_bipartite_graph_maximum_cardinality_matching),
+		// TEST_FUNCTION_ENTRY(test_bipartite_graph_minimum_vertex_cover),
+		// TEST_FUNCTION_ENTRY(test_integer_hermite_normal_form),
+		// TEST_FUNCTION_ENTRY(test_integer_backsubstitute),
+		// TEST_FUNCTION_ENTRY(test_lanczos_iteration_d),
+		// TEST_FUNCTION_ENTRY(test_lanczos_iteration_z),
+		// TEST_FUNCTION_ENTRY(test_eigensystem_krylov_symmetric),
+		// TEST_FUNCTION_ENTRY(test_eigensystem_krylov_hermitian),
+		// TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_basic),
+		// TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_advanced),
+		// TEST_FUNCTION_ENTRY(test_mpo_from_assembly),
+		// TEST_FUNCTION_ENTRY(test_ttno_graph_from_opchains),
+		// TEST_FUNCTION_ENTRY(test_ttno_from_assembly),
+		// TEST_FUNCTION_ENTRY(test_ising_1d_mpo),
+		// TEST_FUNCTION_ENTRY(test_heisenberg_xxz_1d_mpo),
+		// TEST_FUNCTION_ENTRY(test_bose_hubbard_1d_mpo),
+		// TEST_FUNCTION_ENTRY(test_fermi_hubbard_1d_mpo),
+		// TEST_FUNCTION_ENTRY(test_molecular_hamiltonian_mpo),
+		// TEST_FUNCTION_ENTRY(test_spin_molecular_hamiltonian_mpo),
+		// TEST_FUNCTION_ENTRY(test_retained_bond_indices),
+		// TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd),
+		// TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd_zero),
+		// TEST_FUNCTION_ENTRY(test_operator_inner_product),
+		// TEST_FUNCTION_ENTRY(test_apply_operator),
+		// TEST_FUNCTION_ENTRY(test_dmrg_singlesite),
+		// TEST_FUNCTION_ENTRY(test_dmrg_twosite),
+		// TEST_FUNCTION_ENTRY(test_operator_average_coefficient_gradient),
 	};
 	int num_tests = sizeof(tests) / sizeof(struct test);
 

From d61a7af37e2324fc412c870a1df159f7d329921e Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Sun, 13 Oct 2024 07:58:15 +0200
Subject: [PATCH 03/16] Use CT_SINGLE_REAL for test

---
 test/mps/test_mps.c | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index aa6f4ff..f8afc25 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -663,8 +663,8 @@ char* test_mps_add()
 		qsite[i] = qnum_sector;
 	}
 	
-	construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &chi);
-	construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &psi);
+	construct_random_mps(CT_SINGLE_REAL, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &chi);
+	construct_random_mps(CT_SINGLE_REAL, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &psi);
 
 	mps_add(&chi, &psi, &res);
 
@@ -672,11 +672,12 @@ char* test_mps_add()
 	mps_to_statevector(&psi, &psi_vec);
 	mps_to_statevector(&res, &res_vec);
 
-	double *chi_data = (double*)(chi_vec.blocks[0]->data);
-	double *psi_data = (double*)(psi_vec.blocks[0]->data);
-	double *res_data = (double*)(res_vec.blocks[0]->data);
+	float *chi_data = (float*)(chi_vec.blocks[0]->data);
+	float *psi_data = (float*)(psi_vec.blocks[0]->data);
+	float *res_data = (float*)(res_vec.blocks[0]->data);
+	
 	for (int i = 0; i < chi_vec.blocks[0]->ndim; i++){
-		assert(chi_data[i] + psi_data[i] == res_data[i]);
+		assert(fabsf(chi_data[i] + psi_data[i] - res_data[i]) < 1e-13);
 	}
 
 	delete_block_sparse_tensor(&res_vec);

From 65ffc1dda03ac94bbaa5be2c25585fd30b54aa1d Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Wed, 16 Oct 2024 09:37:55 +0200
Subject: [PATCH 04/16] Re-work L=1 case and unit test function

---
 src/mps/mps.c       | 116 ++++++++++++++++++++++++++++++++++----------
 test/mps/test_mps.c |  74 ++++++++++++++++------------
 test/run_tests.c    |   2 +-
 3 files changed, 134 insertions(+), 58 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index 5fad6a5..830af5d 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -315,6 +315,15 @@ void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret)
 	delete_block_sparse_tensor(&r);
 }
 
+// TODO: Delete.
+void print_sparse(struct block_sparse_tensor r)
+{
+	printf("r: ndim=%d, lshape=(%ld, %ld, %ld), bshape=(%ld, %ld, %ld)\n", 
+		r.ndim, 
+		r.dim_logical[0], r.dim_logical[1], r.dim_logical[2],
+		r.dim_blocks[0], r.dim_blocks[1], r.dim_blocks[2]);
+}
+
 
 //________________________________________________________________________________________________________________________
 ///
@@ -335,48 +344,103 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 	long d = chi->d;
 	int L = chi->nsites;
 	enum numeric_type dtype = chi->a[0].dtype;
+	
+	// Copy physical quantum number array.
 	qnumber *qsite = ct_malloc(d * sizeof(qnumber));
 	for (int i = 0; i < d; i++) {
 		qsite[i] = chi->qsite[i];
 	}
 
+	// Initialize return mps.
+	struct mps VALUE;
+	VALUE.nsites = L;
+	VALUE.d = d;
+	VALUE.qsite = qsite;
+
 	if (L == 1) {
-		const struct block_sparse_tensor* chi_a = &(chi->a[0]);
-		const struct block_sparse_tensor* psi_a = &(psi->a[0]);
-		const int ndim = chi_a->ndim;
+		struct block_sparse_tensor chi_a, psi_a;
+
+		copy_block_sparse_tensor(&(chi->a[0]), &chi_a);
+		copy_block_sparse_tensor(&(psi->a[0]), &psi_a);
+
+		assert(chi_a.ndim == psi_a.ndim);
+
+		const int ndim = chi_a.ndim;
 
         // Dummy bond quantum numbers must agree.
-		assert(
-			qnumber_all_equal(
+		assert(qnumber_all_equal(
 				ndim, 
-				chi_a->qnums_logical[0], 
-				psi_a->qnums_logical[0])); 
-		assert(
-			qnumber_all_equal(
+				chi_a.qnums_blocks[0], 
+				psi_a.qnums_blocks[0])); 
+		assert(qnumber_all_equal(
 				ndim,
-				chi_a->qnums_logical[1],
-				psi_a->qnums_logical[1])); 
-
-		// Add MPS tensors. 
-		struct block_sparse_tensor *ret_a = ct_malloc(sizeof(struct block_sparse_tensor));
-		copy_block_sparse_tensor(chi_a, ret_a);
-
-		// Add individual dense_tensors.
-		for (int i = 0; i < ret_a->ndim; i++) {
-			if (psi_a->blocks[i] != NULL && ret_a->blocks[i] != NULL) {
-				dense_tensor_scalar_multiply_add(
-					numeric_one(psi_a->dtype),
-					psi_a->blocks[i],
-					ret_a->blocks[i]);
+				chi_a.qnums_blocks[2],
+				psi_a.qnums_blocks[2]));
+
+		// Add individual dense_tensors (chi = psi + chi).
+		const long nblocks = integer_product(chi_a.dim_blocks, ndim);
+		for (long k = 0; k < nblocks; k++)
+		{
+			struct dense_tensor* a = psi_a.blocks[k];
+			struct dense_tensor* b = chi_a.blocks[k];
+			if (a != NULL && b != NULL) {
+				dense_tensor_scalar_multiply_add(numeric_one(dtype), a, b);
 			}
 		}
 
-		// Allocate memory and set return value 'ret'.
-		// ret = ct_malloc(sizeof(struct mps));
-		struct mps VALUE = {.nsites=L, .d=d, .qsite=qsite, .a=ret_a};
+		delete_block_sparse_tensor(&psi_a);
+
+		// Construct "return value".		
+		VALUE.a = ct_malloc(sizeof(struct block_sparse_tensor));
+		copy_block_sparse_tensor(&chi_a, &VALUE.a[0]);
 		memcpy((void*)ret, (void*)&VALUE, sizeof(struct mps));
 	} else {
+		// TODO: Leading and trailing (dummy) bond quantum numbers must agree.
+		// ...
+		
+		struct block_sparse_tensor *ret_a = ct_malloc(L * sizeof(struct block_sparse_tensor));
+		struct block_sparse_tensor tlist[2];
+		
+		// Left-most tensor.
+		{
+			const int i_ax[1] = { 2 };
+			tlist[0] = chi->a[0];
+			tlist[1] = psi->a[0];
+
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret_a[0]);
+		}
+
+		// Intermediate tensors.
+		for (int i = 1; i < L - 1; i++) {
+			const int i_ax[1] = { 2 };
+			const int i_ax2[1] = { 0 };
+			struct block_sparse_tensor rt, rb;
+			
+			tlist[0] = chi->a[i];
+			allocate_block_sparse_tensor_like(&chi->a[i], &tlist[1]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &rt);
+
+			allocate_block_sparse_tensor_like(&chi->a[i], &tlist[0]);
+			tlist[1] = psi->a[i];
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &rb);
+
+			tlist[0] = rt;
+			tlist[1] = rb;
+			block_sparse_tensor_block_diag(tlist, 2, i_ax2, 1, &ret_a[i]);
+		}
 
+		// Right-most tensor.
+		{
+			const int i_ax[1] = { 0 };
+			tlist[0] = chi->a[L - 1];
+			tlist[1] = psi->a[L - 1];
+			
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret_a[L - 1]);
+		}
+
+		// Construct "return value".		
+		VALUE.a = ret_a;
+		memcpy((void*)ret, (void*)&VALUE, sizeof(struct mps));
 	}
 }
 
diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index f8afc25..42a7ba0 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -645,48 +645,60 @@ char* test_mps_to_statevector()
 
 char* test_mps_add()
 {
+	const long D = 2, MAX_VDIM = 128, SITE_ARRAY_LEN = 6;
+	const int SITE_ARRAY[SITE_ARRAY_LEN] = {1, 3, 5, 8, 12, 16};
+	
+	const qnumber q_pnum = 7;
+	const qnumber q_spin = 1;
+	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
+
+	const qnumber QSITE_CHI[D] = {1, 2};
+	const qnumber QSITE_PSI[D] = {1, 2};
+
 	struct mps chi, psi, res;
 	struct block_sparse_tensor chi_vec, psi_vec, res_vec;
-	
-	const int nsites = 1;
-	const long d = 2, max_vdim = 16;
+	struct dense_tensor chi_vec_dns, psi_vec_dns, res_vec_dns;
 	
 	struct rng_state rng_state;
 	seed_rng_state(42, &rng_state);
+	
+	int nsites;
+	for (int i = 0; i < SITE_ARRAY_LEN; i++) {
+		nsites = SITE_ARRAY[i];
 
-	const qnumber q_pnum = 7;
-	const qnumber q_spin = 1;
-	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
+		construct_random_mps(CT_SINGLE_REAL, nsites, D, QSITE_CHI, qnum_sector, MAX_VDIM, &rng_state, &chi);
+		construct_random_mps(CT_SINGLE_REAL, nsites, D, QSITE_PSI, qnum_sector, MAX_VDIM, &rng_state, &psi);
 
-	qnumber *qsite = ct_malloc(d * sizeof(qnumber));
-	for (int i = 0; i < d; i++) {
-		qsite[i] = qnum_sector;
-	}
-	
-	construct_random_mps(CT_SINGLE_REAL, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &chi);
-	construct_random_mps(CT_SINGLE_REAL, nsites, d, qsite, qnum_sector, max_vdim, &rng_state, &psi);
+		mps_add(&chi, &psi, &res);
 
-	mps_add(&chi, &psi, &res);
+		mps_to_statevector(&chi, &chi_vec);
+		mps_to_statevector(&psi, &psi_vec);
+		mps_to_statevector(&res, &res_vec);
 
-	mps_to_statevector(&chi, &chi_vec);
-	mps_to_statevector(&psi, &psi_vec);
-	mps_to_statevector(&res, &res_vec);
+		block_sparse_to_dense_tensor(&chi_vec, &chi_vec_dns);
+		block_sparse_to_dense_tensor(&psi_vec, &psi_vec_dns);
+		block_sparse_to_dense_tensor(&res_vec, &res_vec_dns);
 
-	float *chi_data = (float*)(chi_vec.blocks[0]->data);
-	float *psi_data = (float*)(psi_vec.blocks[0]->data);
-	float *res_data = (float*)(res_vec.blocks[0]->data);
-	
-	for (int i = 0; i < chi_vec.blocks[0]->ndim; i++){
-		assert(fabsf(chi_data[i] + psi_data[i] - res_data[i]) < 1e-13);
-	}
+		dense_tensor_scalar_multiply_add(numeric_one(chi_vec.dtype), &psi_vec_dns, &chi_vec_dns);
 
-	delete_block_sparse_tensor(&res_vec);
-	delete_block_sparse_tensor(&psi_vec);
-	delete_block_sparse_tensor(&chi_vec);
-	ct_free(qsite);
-	delete_mps(&res);
-	delete_mps(&psi);
-	delete_mps(&chi);
+		// Compare dense tensors.
+		if (!dense_tensor_allclose(&res_vec_dns, &chi_vec_dns, 1e-13)) {
+			return "addition of mps does not match reference";
+		}
+
+		// Free memory.
+		delete_dense_tensor(&res_vec_dns);
+		delete_dense_tensor(&psi_vec_dns);
+		delete_dense_tensor(&chi_vec_dns);
+
+		delete_block_sparse_tensor(&res_vec);
+		delete_block_sparse_tensor(&psi_vec);
+		delete_block_sparse_tensor(&chi_vec);
+		
+		delete_mps(&res);
+		delete_mps(&psi);
+		delete_mps(&chi);
+	}
 
 	return 0;
 }
\ No newline at end of file
diff --git a/test/run_tests.c b/test/run_tests.c
index 410b0c1..29f422d 100644
--- a/test/run_tests.c
+++ b/test/run_tests.c
@@ -139,7 +139,7 @@ int main()
 		// TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr),
 		// TEST_FUNCTION_ENTRY(test_mps_compress),
 		// TEST_FUNCTION_ENTRY(test_mps_split_tensor_svd),
-		// TEST_FUNCTION_ENTRY(test_mps_to_statevector),
+		TEST_FUNCTION_ENTRY(test_mps_to_statevector),
 		TEST_FUNCTION_ENTRY(test_mps_add),
 		// TEST_FUNCTION_ENTRY(test_queue),
 		// TEST_FUNCTION_ENTRY(test_linked_list),

From a1c6e69b5433f728990c6cee01a5727c61d851be Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Wed, 16 Oct 2024 09:38:56 +0200
Subject: [PATCH 05/16] Add .vscode/ to .gitignore

---
 .gitignore | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index e945163..85b2be8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -39,8 +39,9 @@
 *.idb
 *.pdb
 
-# Visual Studio folders
+# Visual Studio and Visual Code folders
 .vs/
+.vscode/
 
 # Python
 __pycache__

From 4e32fdd820b71e54ebf8914c5220e11decef14bb Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Wed, 16 Oct 2024 10:37:10 +0200
Subject: [PATCH 06/16] Cleanup code

---
 src/mps/mps.c | 92 ++++++++++++++++++++-------------------------------
 1 file changed, 36 insertions(+), 56 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index 830af5d..02c5df4 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -315,15 +315,6 @@ void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret)
 	delete_block_sparse_tensor(&r);
 }
 
-// TODO: Delete.
-void print_sparse(struct block_sparse_tensor r)
-{
-	printf("r: ndim=%d, lshape=(%ld, %ld, %ld), bshape=(%ld, %ld, %ld)\n", 
-		r.ndim, 
-		r.dim_logical[0], r.dim_logical[1], r.dim_logical[2],
-		r.dim_blocks[0], r.dim_blocks[1], r.dim_blocks[2]);
-}
-
 
 //________________________________________________________________________________________________________________________
 ///
@@ -341,32 +332,29 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 	assert(chi->d == psi->d); 
 	assert(qnumber_all_equal(chi->d, chi->qsite, psi->qsite)); 
 
-	long d = chi->d;
-	int L = chi->nsites;
-	enum numeric_type dtype = chi->a[0].dtype;
+	const long d = chi->d;
+	const int L = chi->nsites;
+	const enum numeric_type dtype = chi->a[0].dtype;
 	
 	// Copy physical quantum number array.
-	qnumber *qsite = ct_malloc(d * sizeof(qnumber));
+	qnumber* qsite = ct_malloc(d * sizeof(qnumber));
 	for (int i = 0; i < d; i++) {
 		qsite[i] = chi->qsite[i];
 	}
 
 	// Initialize return mps.
-	struct mps VALUE;
-	VALUE.nsites = L;
-	VALUE.d = d;
-	VALUE.qsite = qsite;
-
-	if (L == 1) {
-		struct block_sparse_tensor chi_a, psi_a;
-
-		copy_block_sparse_tensor(&(chi->a[0]), &chi_a);
-		copy_block_sparse_tensor(&(psi->a[0]), &psi_a);
+	struct mps res;
+	res.nsites = L;
+	res.d = d;
+	res.qsite = qsite;
+
+	if (L == 1) {		
+		const struct block_sparse_tensor chi_a = chi->a[0];
+		const struct block_sparse_tensor psi_a = psi->a[0];
+		const int ndim = chi_a.ndim;
 
 		assert(chi_a.ndim == psi_a.ndim);
 
-		const int ndim = chi_a.ndim;
-
         // Dummy bond quantum numbers must agree.
 		assert(qnumber_all_equal(
 				ndim, 
@@ -377,28 +365,32 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 				chi_a.qnums_blocks[2],
 				psi_a.qnums_blocks[2]));
 
-		// Add individual dense_tensors (chi = psi + chi).
-		const long nblocks = integer_product(chi_a.dim_blocks, ndim);
+		// Allocate memory and copy sparse tensor into resulting tensor.
+		res.a = ct_malloc(sizeof(struct block_sparse_tensor));
+		copy_block_sparse_tensor(&chi_a, &res.a[0]);
+
+		// Add individual dense_tensors.
+		const long nblocks = integer_product(res.a->dim_blocks, ndim);
 		for (long k = 0; k < nblocks; k++)
 		{
 			struct dense_tensor* a = psi_a.blocks[k];
-			struct dense_tensor* b = chi_a.blocks[k];
+			struct dense_tensor* b = res.a->blocks[k];
 			if (a != NULL && b != NULL) {
 				dense_tensor_scalar_multiply_add(numeric_one(dtype), a, b);
 			}
 		}
-
-		delete_block_sparse_tensor(&psi_a);
-
-		// Construct "return value".		
-		VALUE.a = ct_malloc(sizeof(struct block_sparse_tensor));
-		copy_block_sparse_tensor(&chi_a, &VALUE.a[0]);
-		memcpy((void*)ret, (void*)&VALUE, sizeof(struct mps));
 	} else {
-		// TODO: Leading and trailing (dummy) bond quantum numbers must agree.
-		// ...
+		// Leading and trailing (dummy) bond quantum numbers must agree.
+		assert(qnumber_all_equal(
+				chi->a[0].ndim, 
+				chi->a[0].qnums_blocks[0], 
+				psi->a[0].qnums_blocks[0])); 
+		assert(qnumber_all_equal(
+				chi->a[chi->nsites - 1].ndim, 
+				chi->a[chi->nsites - 1].qnums_blocks[2], 
+				psi->a[chi->nsites - 1].qnums_blocks[2])); 
 		
-		struct block_sparse_tensor *ret_a = ct_malloc(L * sizeof(struct block_sparse_tensor));
+		res.a = ct_malloc(L * sizeof(struct block_sparse_tensor));
 		struct block_sparse_tensor tlist[2];
 		
 		// Left-most tensor.
@@ -407,26 +399,16 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			tlist[0] = chi->a[0];
 			tlist[1] = psi->a[0];
 
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret_a[0]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &res.a[0]);
 		}
 
 		// Intermediate tensors.
 		for (int i = 1; i < L - 1; i++) {
-			const int i_ax[1] = { 2 };
-			const int i_ax2[1] = { 0 };
-			struct block_sparse_tensor rt, rb;
-			
+			const int i_ax[2] = { 0, 2 };
 			tlist[0] = chi->a[i];
-			allocate_block_sparse_tensor_like(&chi->a[i], &tlist[1]);
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &rt);
-
-			allocate_block_sparse_tensor_like(&chi->a[i], &tlist[0]);
 			tlist[1] = psi->a[i];
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &rb);
 
-			tlist[0] = rt;
-			tlist[1] = rb;
-			block_sparse_tensor_block_diag(tlist, 2, i_ax2, 1, &ret_a[i]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 2, &res.a[i]);
 		}
 
 		// Right-most tensor.
@@ -435,13 +417,11 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			tlist[0] = chi->a[L - 1];
 			tlist[1] = psi->a[L - 1];
 			
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret_a[L - 1]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &res.a[L - 1]);
 		}
-
-		// Construct "return value".		
-		VALUE.a = ret_a;
-		memcpy((void*)ret, (void*)&VALUE, sizeof(struct mps));
 	}
+
+	memcpy((void*)ret, (void*)&res, sizeof(struct mps));
 }
 
 

From 7f168ae6867a1702c0e319ae34c6085ecb3d9c76 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Wed, 16 Oct 2024 10:39:33 +0200
Subject: [PATCH 07/16] Undo uncomment in run_tests.c

---
 test/mps/test_mps.c |  16 ++---
 test/run_tests.c    | 148 ++++++++++++++++++++++----------------------
 2 files changed, 82 insertions(+), 82 deletions(-)

diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index 42a7ba0..e27ea41 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -645,15 +645,15 @@ char* test_mps_to_statevector()
 
 char* test_mps_add()
 {
-	const long D = 2, MAX_VDIM = 128, SITE_ARRAY_LEN = 6;
-	const int SITE_ARRAY[SITE_ARRAY_LEN] = {1, 3, 5, 8, 12, 16};
+	const long d = 2, max_vdim = 16, site_array_len = 6;
+	const int site_array[site_array_len] = {1, 3, 5, 8, 12, 16};
 	
 	const qnumber q_pnum = 7;
 	const qnumber q_spin = 1;
 	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
 
-	const qnumber QSITE_CHI[D] = {1, 2};
-	const qnumber QSITE_PSI[D] = {1, 2};
+	const qnumber qsite_chi[d] = {1, 2};
+	const qnumber qsite_psi[d] = {1, 2};
 
 	struct mps chi, psi, res;
 	struct block_sparse_tensor chi_vec, psi_vec, res_vec;
@@ -663,11 +663,11 @@ char* test_mps_add()
 	seed_rng_state(42, &rng_state);
 	
 	int nsites;
-	for (int i = 0; i < SITE_ARRAY_LEN; i++) {
-		nsites = SITE_ARRAY[i];
+	for (int i = 0; i < site_array_len; i++) {
+		nsites = site_array[i];
 
-		construct_random_mps(CT_SINGLE_REAL, nsites, D, QSITE_CHI, qnum_sector, MAX_VDIM, &rng_state, &chi);
-		construct_random_mps(CT_SINGLE_REAL, nsites, D, QSITE_PSI, qnum_sector, MAX_VDIM, &rng_state, &psi);
+		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_chi, qnum_sector, max_vdim, &rng_state, &chi);
+		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_psi, qnum_sector, max_vdim, &rng_state, &psi);
 
 		mps_add(&chi, &psi, &res);
 
diff --git a/test/run_tests.c b/test/run_tests.c
index 29f422d..0083350 100644
--- a/test/run_tests.c
+++ b/test/run_tests.c
@@ -95,82 +95,82 @@ char* test_operator_average_coefficient_gradient();
 int main()
 {
 	struct test tests[] = {
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_trace),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_cyclic_partial_trace),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_transpose),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_slice),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_pad_zeros),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_pointwise),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_axis),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_dot),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_dot_update),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product_degree_zero),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_concatenate),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_block_diag),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_qr),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_rq),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_svd),
-		// TEST_FUNCTION_ENTRY(test_dense_tensor_block),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_copy),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_block),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_cyclic_partial_trace),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_norm2),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_transpose),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_reshape),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_slice),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_pointwise_vector),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_axis),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_dot),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_concatenate),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_block_diag),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_qr),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_rq),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_svd),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_serialize),
-		// TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_entry),
-		// TEST_FUNCTION_ENTRY(test_clebsch_gordan_coefficients),
-		// TEST_FUNCTION_ENTRY(test_su2_tree_enumerate_charge_sectors),
-		// TEST_FUNCTION_ENTRY(test_su2_fuse_split_tree_enumerate_charge_sectors),
-		// TEST_FUNCTION_ENTRY(test_su2_tensor_fmove),
-		// TEST_FUNCTION_ENTRY(test_su2_tensor_contract_simple),
-		// TEST_FUNCTION_ENTRY(test_su2_to_dense_tensor),
-		// TEST_FUNCTION_ENTRY(test_mps_vdot),
-		// TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr),
-		// TEST_FUNCTION_ENTRY(test_mps_compress),
-		// TEST_FUNCTION_ENTRY(test_mps_split_tensor_svd),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_trace),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_cyclic_partial_trace),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_transpose),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_slice),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_pad_zeros),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_pointwise),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_multiply_axis),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_dot),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_dot_update),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_kronecker_product_degree_zero),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_concatenate),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_block_diag),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_qr),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_rq),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_svd),
+		TEST_FUNCTION_ENTRY(test_dense_tensor_block),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_copy),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_block),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_cyclic_partial_trace),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_norm2),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_transpose),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_reshape),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_slice),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_pointwise_vector),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_multiply_axis),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_dot),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_concatenate),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_block_diag),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_qr),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_rq),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_svd),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_serialize),
+		TEST_FUNCTION_ENTRY(test_block_sparse_tensor_get_entry),
+		TEST_FUNCTION_ENTRY(test_clebsch_gordan_coefficients),
+		TEST_FUNCTION_ENTRY(test_su2_tree_enumerate_charge_sectors),
+		TEST_FUNCTION_ENTRY(test_su2_fuse_split_tree_enumerate_charge_sectors),
+		TEST_FUNCTION_ENTRY(test_su2_tensor_fmove),
+		TEST_FUNCTION_ENTRY(test_su2_tensor_contract_simple),
+		TEST_FUNCTION_ENTRY(test_su2_to_dense_tensor),
+		TEST_FUNCTION_ENTRY(test_mps_vdot),
+		TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr),
+		TEST_FUNCTION_ENTRY(test_mps_compress),
+		TEST_FUNCTION_ENTRY(test_mps_split_tensor_svd),
 		TEST_FUNCTION_ENTRY(test_mps_to_statevector),
 		TEST_FUNCTION_ENTRY(test_mps_add),
-		// TEST_FUNCTION_ENTRY(test_queue),
-		// TEST_FUNCTION_ENTRY(test_linked_list),
-		// TEST_FUNCTION_ENTRY(test_hash_table),
-		// TEST_FUNCTION_ENTRY(test_bipartite_graph_maximum_cardinality_matching),
-		// TEST_FUNCTION_ENTRY(test_bipartite_graph_minimum_vertex_cover),
-		// TEST_FUNCTION_ENTRY(test_integer_hermite_normal_form),
-		// TEST_FUNCTION_ENTRY(test_integer_backsubstitute),
-		// TEST_FUNCTION_ENTRY(test_lanczos_iteration_d),
-		// TEST_FUNCTION_ENTRY(test_lanczos_iteration_z),
-		// TEST_FUNCTION_ENTRY(test_eigensystem_krylov_symmetric),
-		// TEST_FUNCTION_ENTRY(test_eigensystem_krylov_hermitian),
-		// TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_basic),
-		// TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_advanced),
-		// TEST_FUNCTION_ENTRY(test_mpo_from_assembly),
-		// TEST_FUNCTION_ENTRY(test_ttno_graph_from_opchains),
-		// TEST_FUNCTION_ENTRY(test_ttno_from_assembly),
-		// TEST_FUNCTION_ENTRY(test_ising_1d_mpo),
-		// TEST_FUNCTION_ENTRY(test_heisenberg_xxz_1d_mpo),
-		// TEST_FUNCTION_ENTRY(test_bose_hubbard_1d_mpo),
-		// TEST_FUNCTION_ENTRY(test_fermi_hubbard_1d_mpo),
-		// TEST_FUNCTION_ENTRY(test_molecular_hamiltonian_mpo),
-		// TEST_FUNCTION_ENTRY(test_spin_molecular_hamiltonian_mpo),
-		// TEST_FUNCTION_ENTRY(test_retained_bond_indices),
-		// TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd),
-		// TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd_zero),
-		// TEST_FUNCTION_ENTRY(test_operator_inner_product),
-		// TEST_FUNCTION_ENTRY(test_apply_operator),
-		// TEST_FUNCTION_ENTRY(test_dmrg_singlesite),
-		// TEST_FUNCTION_ENTRY(test_dmrg_twosite),
-		// TEST_FUNCTION_ENTRY(test_operator_average_coefficient_gradient),
+		TEST_FUNCTION_ENTRY(test_queue),
+		TEST_FUNCTION_ENTRY(test_linked_list),
+		TEST_FUNCTION_ENTRY(test_hash_table),
+		TEST_FUNCTION_ENTRY(test_bipartite_graph_maximum_cardinality_matching),
+		TEST_FUNCTION_ENTRY(test_bipartite_graph_minimum_vertex_cover),
+		TEST_FUNCTION_ENTRY(test_integer_hermite_normal_form),
+		TEST_FUNCTION_ENTRY(test_integer_backsubstitute),
+		TEST_FUNCTION_ENTRY(test_lanczos_iteration_d),
+		TEST_FUNCTION_ENTRY(test_lanczos_iteration_z),
+		TEST_FUNCTION_ENTRY(test_eigensystem_krylov_symmetric),
+		TEST_FUNCTION_ENTRY(test_eigensystem_krylov_hermitian),
+		TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_basic),
+		TEST_FUNCTION_ENTRY(test_mpo_graph_from_opchains_advanced),
+		TEST_FUNCTION_ENTRY(test_mpo_from_assembly),
+		TEST_FUNCTION_ENTRY(test_ttno_graph_from_opchains),
+		TEST_FUNCTION_ENTRY(test_ttno_from_assembly),
+		TEST_FUNCTION_ENTRY(test_ising_1d_mpo),
+		TEST_FUNCTION_ENTRY(test_heisenberg_xxz_1d_mpo),
+		TEST_FUNCTION_ENTRY(test_bose_hubbard_1d_mpo),
+		TEST_FUNCTION_ENTRY(test_fermi_hubbard_1d_mpo),
+		TEST_FUNCTION_ENTRY(test_molecular_hamiltonian_mpo),
+		TEST_FUNCTION_ENTRY(test_spin_molecular_hamiltonian_mpo),
+		TEST_FUNCTION_ENTRY(test_retained_bond_indices),
+		TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd),
+		TEST_FUNCTION_ENTRY(test_split_block_sparse_matrix_svd_zero),
+		TEST_FUNCTION_ENTRY(test_operator_inner_product),
+		TEST_FUNCTION_ENTRY(test_apply_operator),
+		TEST_FUNCTION_ENTRY(test_dmrg_singlesite),
+		TEST_FUNCTION_ENTRY(test_dmrg_twosite),
+		TEST_FUNCTION_ENTRY(test_operator_average_coefficient_gradient),
 	};
 	int num_tests = sizeof(tests) / sizeof(struct test);
 

From aee96aab6e89ef8e201ed7e0f31194249f15b115 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Wed, 16 Oct 2024 10:40:17 +0200
Subject: [PATCH 08/16] Add dependency installation instructions to README.md

---
 README.md | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/README.md b/README.md
index 4f9279b..5bab8ea 100644
--- a/README.md
+++ b/README.md
@@ -9,14 +9,16 @@ Tensor network algorithms for chemical systems.
 
 Building
 --------
-The code requires the BLAS, LAPACKE, HDF5 and Python 3 development libraries with NumPy. These can be installed via `sudo apt install libblas-dev liblapacke-dev libhdf5-dev python3-dev python3-numpy` (on Ubuntu Linux) or similar.
+The code requires the BLAS, LAPACKE, HDF5 and Python 3 development libraries with NumPy. These can be installed via 
+- `sudo apt install libblas-dev liblapacke-dev libhdf5-dev python3-dev python3-numpy` (on Ubuntu Linux)
+- `brew install openblas lapack hdf5 python3 numpy` (on MacOS)
 
 From the project directory, use `cmake` to build the project:
-```
+```bash
 mkdir build && cd build
 cmake ../
 cmake --build .
-````
+```
 
 Currently, this will compile the unit tests, which you can run via `./chemtensor_test`, as well as the demo examples and Python module library.
 

From 898dea7be6c0af1e6b5cd45e27c0b5c1f0c2a783 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Thu, 17 Oct 2024 07:00:29 +0200
Subject: [PATCH 09/16] lowercase first letter in comments and remove dots

---
 src/mps/mps.c       | 24 ++++++++++++------------
 test/mps/test_mps.c |  5 +++--
 2 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index 02c5df4..94fa2c2 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -322,13 +322,13 @@ void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret)
 ///
 void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 {
-	// Number of lattice sites must agree.
+	// number of lattice sites must agree
 	assert(chi->nsites == psi->nsites);
 
-	// Number of lattice sites must be larger than 0.
+	// number of lattice sites must be larger than 0
 	assert(chi->nsites > 0);
 
-	// Physical quantum numbers must agree.
+	// physical quantum numbers must agree
 	assert(chi->d == psi->d); 
 	assert(qnumber_all_equal(chi->d, chi->qsite, psi->qsite)); 
 
@@ -336,13 +336,13 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 	const int L = chi->nsites;
 	const enum numeric_type dtype = chi->a[0].dtype;
 	
-	// Copy physical quantum number array.
+	// copy physical quantum number array
 	qnumber* qsite = ct_malloc(d * sizeof(qnumber));
 	for (int i = 0; i < d; i++) {
 		qsite[i] = chi->qsite[i];
 	}
 
-	// Initialize return mps.
+	// initialize return mps
 	struct mps res;
 	res.nsites = L;
 	res.d = d;
@@ -355,7 +355,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 
 		assert(chi_a.ndim == psi_a.ndim);
 
-        // Dummy bond quantum numbers must agree.
+		// dummy bond quantum numbers must agree
 		assert(qnumber_all_equal(
 				ndim, 
 				chi_a.qnums_blocks[0], 
@@ -365,11 +365,11 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 				chi_a.qnums_blocks[2],
 				psi_a.qnums_blocks[2]));
 
-		// Allocate memory and copy sparse tensor into resulting tensor.
+		// allocate memory and copy sparse tensor into resulting tensor
 		res.a = ct_malloc(sizeof(struct block_sparse_tensor));
 		copy_block_sparse_tensor(&chi_a, &res.a[0]);
 
-		// Add individual dense_tensors.
+		// add individual dense_tensors
 		const long nblocks = integer_product(res.a->dim_blocks, ndim);
 		for (long k = 0; k < nblocks; k++)
 		{
@@ -380,7 +380,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			}
 		}
 	} else {
-		// Leading and trailing (dummy) bond quantum numbers must agree.
+		// leading and trailing (dummy) bond quantum numbers must agree
 		assert(qnumber_all_equal(
 				chi->a[0].ndim, 
 				chi->a[0].qnums_blocks[0], 
@@ -393,7 +393,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 		res.a = ct_malloc(L * sizeof(struct block_sparse_tensor));
 		struct block_sparse_tensor tlist[2];
 		
-		// Left-most tensor.
+		// left-most tensor
 		{
 			const int i_ax[1] = { 2 };
 			tlist[0] = chi->a[0];
@@ -402,7 +402,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &res.a[0]);
 		}
 
-		// Intermediate tensors.
+		// intermediate tensors
 		for (int i = 1; i < L - 1; i++) {
 			const int i_ax[2] = { 0, 2 };
 			tlist[0] = chi->a[i];
@@ -411,7 +411,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			block_sparse_tensor_block_diag(tlist, 2, i_ax, 2, &res.a[i]);
 		}
 
-		// Right-most tensor.
+		// right-most tensor
 		{
 			const int i_ax[1] = { 0 };
 			tlist[0] = chi->a[L - 1];
diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index e27ea41..b76bbdd 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -643,6 +643,7 @@ char* test_mps_to_statevector()
 	return 0;
 }
 
+
 char* test_mps_add()
 {
 	const long d = 2, max_vdim = 16, site_array_len = 6;
@@ -681,12 +682,12 @@ char* test_mps_add()
 
 		dense_tensor_scalar_multiply_add(numeric_one(chi_vec.dtype), &psi_vec_dns, &chi_vec_dns);
 
-		// Compare dense tensors.
+		// compare dense tensors
 		if (!dense_tensor_allclose(&res_vec_dns, &chi_vec_dns, 1e-13)) {
 			return "addition of mps does not match reference";
 		}
 
-		// Free memory.
+		// free memory
 		delete_dense_tensor(&res_vec_dns);
 		delete_dense_tensor(&psi_vec_dns);
 		delete_dense_tensor(&chi_vec_dns);

From 66b72bcd2ff2f94690c406862c7c8d614084edbe Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Thu, 17 Oct 2024 07:02:34 +0200
Subject: [PATCH 10/16] Use literal number instead of site_array_len

---
 test/mps/test_mps.c | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index b76bbdd..b5ceaa6 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -646,8 +646,8 @@ char* test_mps_to_statevector()
 
 char* test_mps_add()
 {
-	const long d = 2, max_vdim = 16, site_array_len = 6;
-	const int site_array[site_array_len] = {1, 3, 5, 8, 12, 16};
+	const long d = 2, max_vdim = 16;
+	const int site_array[6] = {1, 3, 5, 8, 12, 16};
 	
 	const qnumber q_pnum = 7;
 	const qnumber q_spin = 1;
@@ -664,7 +664,7 @@ char* test_mps_add()
 	seed_rng_state(42, &rng_state);
 	
 	int nsites;
-	for (int i = 0; i < site_array_len; i++) {
+	for (int i = 0; i < 6; i++) {
 		nsites = site_array[i];
 
 		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_chi, qnum_sector, max_vdim, &rng_state, &chi);

From d219e63b63a49bc94674d794f107ab405019651e Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Thu, 17 Oct 2024 07:04:11 +0200
Subject: [PATCH 11/16] Move test_mps_add below test_mps_vdot

---
 test/mps/test_mps.c | 122 ++++++++++++++++++++++----------------------
 test/run_tests.c    |   4 +-
 2 files changed, 63 insertions(+), 63 deletions(-)

diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index b5ceaa6..4da4d9e 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -127,6 +127,67 @@ char* test_mps_vdot()
 }
 
 
+char* test_mps_add()
+{
+	const long d = 2, max_vdim = 16;
+	const int site_array[6] = {1, 3, 5, 8, 12, 16};
+	
+	const qnumber q_pnum = 7;
+	const qnumber q_spin = 1;
+	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
+
+	const qnumber qsite_chi[d] = {1, 2};
+	const qnumber qsite_psi[d] = {1, 2};
+
+	struct mps chi, psi, res;
+	struct block_sparse_tensor chi_vec, psi_vec, res_vec;
+	struct dense_tensor chi_vec_dns, psi_vec_dns, res_vec_dns;
+	
+	struct rng_state rng_state;
+	seed_rng_state(42, &rng_state);
+	
+	int nsites;
+	for (int i = 0; i < 6; i++) {
+		nsites = site_array[i];
+
+		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_chi, qnum_sector, max_vdim, &rng_state, &chi);
+		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_psi, qnum_sector, max_vdim, &rng_state, &psi);
+
+		mps_add(&chi, &psi, &res);
+
+		mps_to_statevector(&chi, &chi_vec);
+		mps_to_statevector(&psi, &psi_vec);
+		mps_to_statevector(&res, &res_vec);
+
+		block_sparse_to_dense_tensor(&chi_vec, &chi_vec_dns);
+		block_sparse_to_dense_tensor(&psi_vec, &psi_vec_dns);
+		block_sparse_to_dense_tensor(&res_vec, &res_vec_dns);
+
+		dense_tensor_scalar_multiply_add(numeric_one(chi_vec.dtype), &psi_vec_dns, &chi_vec_dns);
+
+		// compare dense tensors
+		if (!dense_tensor_allclose(&res_vec_dns, &chi_vec_dns, 1e-13)) {
+			return "addition of mps does not match reference";
+		}
+
+		// free memory
+		delete_dense_tensor(&res_vec_dns);
+		delete_dense_tensor(&psi_vec_dns);
+		delete_dense_tensor(&chi_vec_dns);
+
+		delete_block_sparse_tensor(&res_vec);
+		delete_block_sparse_tensor(&psi_vec);
+		delete_block_sparse_tensor(&chi_vec);
+		
+		delete_mps(&res);
+		delete_mps(&psi);
+		delete_mps(&chi);
+	}
+
+	return 0;
+}
+
+
 char* test_mps_orthonormalize_qr()
 {
 	hid_t file = H5Fopen("../test/mps/data/test_mps_orthonormalize_qr.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT);
@@ -640,66 +701,5 @@ char* test_mps_to_statevector()
 
 	H5Fclose(file);
 
-	return 0;
-}
-
-
-char* test_mps_add()
-{
-	const long d = 2, max_vdim = 16;
-	const int site_array[6] = {1, 3, 5, 8, 12, 16};
-	
-	const qnumber q_pnum = 7;
-	const qnumber q_spin = 1;
-	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
-
-	const qnumber qsite_chi[d] = {1, 2};
-	const qnumber qsite_psi[d] = {1, 2};
-
-	struct mps chi, psi, res;
-	struct block_sparse_tensor chi_vec, psi_vec, res_vec;
-	struct dense_tensor chi_vec_dns, psi_vec_dns, res_vec_dns;
-	
-	struct rng_state rng_state;
-	seed_rng_state(42, &rng_state);
-	
-	int nsites;
-	for (int i = 0; i < 6; i++) {
-		nsites = site_array[i];
-
-		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_chi, qnum_sector, max_vdim, &rng_state, &chi);
-		construct_random_mps(CT_DOUBLE_COMPLEX, nsites, d, qsite_psi, qnum_sector, max_vdim, &rng_state, &psi);
-
-		mps_add(&chi, &psi, &res);
-
-		mps_to_statevector(&chi, &chi_vec);
-		mps_to_statevector(&psi, &psi_vec);
-		mps_to_statevector(&res, &res_vec);
-
-		block_sparse_to_dense_tensor(&chi_vec, &chi_vec_dns);
-		block_sparse_to_dense_tensor(&psi_vec, &psi_vec_dns);
-		block_sparse_to_dense_tensor(&res_vec, &res_vec_dns);
-
-		dense_tensor_scalar_multiply_add(numeric_one(chi_vec.dtype), &psi_vec_dns, &chi_vec_dns);
-
-		// compare dense tensors
-		if (!dense_tensor_allclose(&res_vec_dns, &chi_vec_dns, 1e-13)) {
-			return "addition of mps does not match reference";
-		}
-
-		// free memory
-		delete_dense_tensor(&res_vec_dns);
-		delete_dense_tensor(&psi_vec_dns);
-		delete_dense_tensor(&chi_vec_dns);
-
-		delete_block_sparse_tensor(&res_vec);
-		delete_block_sparse_tensor(&psi_vec);
-		delete_block_sparse_tensor(&chi_vec);
-		
-		delete_mps(&res);
-		delete_mps(&psi);
-		delete_mps(&chi);
-	}
-
 	return 0;
 }
\ No newline at end of file
diff --git a/test/run_tests.c b/test/run_tests.c
index 0083350..81cd5a9 100644
--- a/test/run_tests.c
+++ b/test/run_tests.c
@@ -52,11 +52,11 @@ char* test_su2_tensor_fmove();
 char* test_su2_tensor_contract_simple();
 char* test_su2_to_dense_tensor();
 char* test_mps_vdot();
+char* test_mps_add();
 char* test_mps_orthonormalize_qr();
 char* test_mps_compress();
 char* test_mps_split_tensor_svd();
 char* test_mps_to_statevector();
-char* test_mps_add();
 char* test_queue();
 char* test_linked_list();
 char* test_hash_table();
@@ -136,11 +136,11 @@ int main()
 		TEST_FUNCTION_ENTRY(test_su2_tensor_contract_simple),
 		TEST_FUNCTION_ENTRY(test_su2_to_dense_tensor),
 		TEST_FUNCTION_ENTRY(test_mps_vdot),
+		TEST_FUNCTION_ENTRY(test_mps_add),
 		TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr),
 		TEST_FUNCTION_ENTRY(test_mps_compress),
 		TEST_FUNCTION_ENTRY(test_mps_split_tensor_svd),
 		TEST_FUNCTION_ENTRY(test_mps_to_statevector),
-		TEST_FUNCTION_ENTRY(test_mps_add),
 		TEST_FUNCTION_ENTRY(test_queue),
 		TEST_FUNCTION_ENTRY(test_linked_list),
 		TEST_FUNCTION_ENTRY(test_hash_table),

From 24c82c88f6c98b40e10c7d1dc308b11cfb885de3 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Thu, 17 Oct 2024 07:10:49 +0200
Subject: [PATCH 12/16] Refactor void* to struct mps* and call
 allocate_empty_mps

---
 src/mps/mps.c | 25 +++++++++----------------
 src/mps/mps.h |  2 +-
 2 files changed, 10 insertions(+), 17 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index 94fa2c2..d092812 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -320,7 +320,7 @@ void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret)
 ///
 /// \brief Compute the addition of two MPS chi and psi.
 ///
-void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
+void mps_add(const struct mps* chi, const struct mps* psi, struct mps* ret)
 {
 	// number of lattice sites must agree
 	assert(chi->nsites == psi->nsites);
@@ -343,10 +343,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 	}
 
 	// initialize return mps
-	struct mps res;
-	res.nsites = L;
-	res.d = d;
-	res.qsite = qsite;
+	allocate_empty_mps(L, d, qsite, ret);
 
 	if (L == 1) {		
 		const struct block_sparse_tensor chi_a = chi->a[0];
@@ -365,16 +362,15 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 				chi_a.qnums_blocks[2],
 				psi_a.qnums_blocks[2]));
 
-		// allocate memory and copy sparse tensor into resulting tensor
-		res.a = ct_malloc(sizeof(struct block_sparse_tensor));
-		copy_block_sparse_tensor(&chi_a, &res.a[0]);
+		// copy sparse tensor into resulting tensor
+		copy_block_sparse_tensor(&chi_a, &ret->a[0]);
 
 		// add individual dense_tensors
-		const long nblocks = integer_product(res.a->dim_blocks, ndim);
+		const long nblocks = integer_product(ret->a->dim_blocks, ndim);
 		for (long k = 0; k < nblocks; k++)
 		{
 			struct dense_tensor* a = psi_a.blocks[k];
-			struct dense_tensor* b = res.a->blocks[k];
+			struct dense_tensor* b = ret->a->blocks[k];
 			if (a != NULL && b != NULL) {
 				dense_tensor_scalar_multiply_add(numeric_one(dtype), a, b);
 			}
@@ -390,7 +386,6 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 				chi->a[chi->nsites - 1].qnums_blocks[2], 
 				psi->a[chi->nsites - 1].qnums_blocks[2])); 
 		
-		res.a = ct_malloc(L * sizeof(struct block_sparse_tensor));
 		struct block_sparse_tensor tlist[2];
 		
 		// left-most tensor
@@ -399,7 +394,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			tlist[0] = chi->a[0];
 			tlist[1] = psi->a[0];
 
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &res.a[0]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret->a[0]);
 		}
 
 		// intermediate tensors
@@ -408,7 +403,7 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			tlist[0] = chi->a[i];
 			tlist[1] = psi->a[i];
 
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 2, &res.a[i]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 2, &ret->a[i]);
 		}
 
 		// right-most tensor
@@ -417,11 +412,9 @@ void mps_add(const struct mps* chi, const struct mps* psi, void* ret)
 			tlist[0] = chi->a[L - 1];
 			tlist[1] = psi->a[L - 1];
 			
-			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &res.a[L - 1]);
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 1, &ret->a[L - 1]);
 		}
 	}
-
-	memcpy((void*)ret, (void*)&res, sizeof(struct mps));
 }
 
 
diff --git a/src/mps/mps.h b/src/mps/mps.h
index e34f6a2..c7f35e4 100644
--- a/src/mps/mps.h
+++ b/src/mps/mps.h
@@ -53,7 +53,7 @@ static inline long mps_bond_dim(const struct mps* mps, const int i)
 
 void mps_vdot(const struct mps* chi, const struct mps* psi, void* ret);
 
-void mps_add(const struct mps* chi, const struct mps* psi, void* ret);
+void mps_add(const struct mps* chi, const struct mps* psi, struct mps* ret);
 
 double mps_norm(const struct mps* psi);
 

From 28270b9971a4e514e7cd6c95d2b986f995b5738d Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Fri, 18 Oct 2024 12:44:07 +0200
Subject: [PATCH 13/16] Remove copying qsite array in add_mps

---
 src/mps/mps.c | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index d092812..09c32b3 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -336,14 +336,8 @@ void mps_add(const struct mps* chi, const struct mps* psi, struct mps* ret)
 	const int L = chi->nsites;
 	const enum numeric_type dtype = chi->a[0].dtype;
 	
-	// copy physical quantum number array
-	qnumber* qsite = ct_malloc(d * sizeof(qnumber));
-	for (int i = 0; i < d; i++) {
-		qsite[i] = chi->qsite[i];
-	}
-
 	// initialize return mps
-	allocate_empty_mps(L, d, qsite, ret);
+	allocate_empty_mps(L, d, chi->qsite, ret);
 
 	if (L == 1) {		
 		const struct block_sparse_tensor chi_a = chi->a[0];

From f4670663e02c9f8314c78b3548ab930a9b2f8978 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Fri, 18 Oct 2024 15:37:03 +0200
Subject: [PATCH 14/16] Use explicit integer for qsite_chi and qsite_psi

---
 test/mps/test_mps.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index 4da4d9e..3845e3e 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -136,8 +136,8 @@ char* test_mps_add()
 	const qnumber q_spin = 1;
 	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
 
-	const qnumber qsite_chi[d] = {1, 2};
-	const qnumber qsite_psi[d] = {1, 2};
+	const qnumber qsite_chi[2] = {1, 2};
+	const qnumber qsite_psi[2] = {1, 2};
 
 	struct mps chi, psi, res;
 	struct block_sparse_tensor chi_vec, psi_vec, res_vec;

From d9048cf78b2ac83da0b8c9407f4b51c0dad14d43 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Sat, 19 Oct 2024 07:30:03 +0200
Subject: [PATCH 15/16] Use dim_logical instead of ndim

---
 src/mps/mps.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/mps/mps.c b/src/mps/mps.c
index 09c32b3..3f3a8f1 100644
--- a/src/mps/mps.c
+++ b/src/mps/mps.c
@@ -348,11 +348,11 @@ void mps_add(const struct mps* chi, const struct mps* psi, struct mps* ret)
 
 		// dummy bond quantum numbers must agree
 		assert(qnumber_all_equal(
-				ndim, 
+				chi_a.dim_logical[0], 
 				chi_a.qnums_blocks[0], 
 				psi_a.qnums_blocks[0])); 
 		assert(qnumber_all_equal(
-				ndim,
+				chi_a.dim_logical[2],
 				chi_a.qnums_blocks[2],
 				psi_a.qnums_blocks[2]));
 
@@ -372,11 +372,11 @@ void mps_add(const struct mps* chi, const struct mps* psi, struct mps* ret)
 	} else {
 		// leading and trailing (dummy) bond quantum numbers must agree
 		assert(qnumber_all_equal(
-				chi->a[0].ndim, 
+				chi->a[0].dim_logical[0], 
 				chi->a[0].qnums_blocks[0], 
-				psi->a[0].qnums_blocks[0])); 
+				psi->a[0].qnums_blocks[0]));
 		assert(qnumber_all_equal(
-				chi->a[chi->nsites - 1].ndim, 
+				chi->a[chi->nsites - 1].dim_logical[2], 
 				chi->a[chi->nsites - 1].qnums_blocks[2], 
 				psi->a[chi->nsites - 1].qnums_blocks[2])); 
 		

From fe9300c59881f870f4af32996f430610473810e4 Mon Sep 17 00:00:00 2001
From: Matthias Reumann <matthias.reumann@tum.de>
Date: Sat, 19 Oct 2024 07:32:09 +0200
Subject: [PATCH 16/16] Set qnum_sector to a small integer

---
 test/mps/test_mps.c | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/test/mps/test_mps.c b/test/mps/test_mps.c
index 3845e3e..be8cd7c 100644
--- a/test/mps/test_mps.c
+++ b/test/mps/test_mps.c
@@ -132,10 +132,7 @@ char* test_mps_add()
 	const long d = 2, max_vdim = 16;
 	const int site_array[6] = {1, 3, 5, 8, 12, 16};
 	
-	const qnumber q_pnum = 7;
-	const qnumber q_spin = 1;
-	const qnumber qnum_sector = encode_quantum_number_pair(q_pnum, q_spin);
-
+	const qnumber qnum_sector = 8;
 	const qnumber qsite_chi[2] = {1, 2};
 	const qnumber qsite_psi[2] = {1, 2};