Add mps_add function

.gitignore

@@ -39,8 +39,9 @@
 *.idb
 *.pdb
 
-# Visual Studio folders
+# Visual Studio and Visual Code folders
 .vs/
+.vscode/
 
 # Python
 __pycache__

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.
 

src/mps/mps.c

@@ -316,6 +316,102 @@ 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, struct mps* 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)); 
+
+	const long d = chi->d;
+	const int L = chi->nsites;
+	const enum numeric_type dtype = chi->a[0].dtype;
+
+	// initialize return mps
+	allocate_empty_mps(L, d, chi->qsite, ret);
+
+	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);
+
+		// dummy bond quantum numbers must agree
+		assert(qnumber_all_equal(
+				chi_a.dim_logical[0], 
+				chi_a.qnums_blocks[0], 
+				psi_a.qnums_blocks[0])); 
+		assert(qnumber_all_equal(
+				chi_a.dim_logical[2],
+				chi_a.qnums_blocks[2],
+				psi_a.qnums_blocks[2]));
+
+		// copy sparse tensor into resulting tensor
+		copy_block_sparse_tensor(&chi_a, &ret->a[0]);
+
+		// add individual dense_tensors
+		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 = ret->a->blocks[k];
+			if (a != NULL && b != NULL) {
+				dense_tensor_scalar_multiply_add(numeric_one(dtype), a, b);
+			}
+		}
+	} else {
+		// leading and trailing (dummy) bond quantum numbers must agree
+		assert(qnumber_all_equal(
+				chi->a[0].dim_logical[0], 
+				chi->a[0].qnums_blocks[0], 
+				psi->a[0].qnums_blocks[0]));
+		assert(qnumber_all_equal(
+				chi->a[chi->nsites - 1].dim_logical[2], 
+				chi->a[chi->nsites - 1].qnums_blocks[2], 
+				psi->a[chi->nsites - 1].qnums_blocks[2])); 
+
+		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[2] = { 0, 2 };
+			tlist[0] = chi->a[i];
+			tlist[1] = psi->a[i];
+
+			block_sparse_tensor_block_diag(tlist, 2, i_ax, 2, &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]);
+		}
+	}
+}
+
+
 //________________________________________________________________________________________________________________________
 ///
 /// \brief Compute the Euclidean norm of the MPS.

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, struct mps* ret);
+
 double mps_norm(const struct mps* psi);
 
 

test/mps/test_mps.c

@@ -127,6 +127,64 @@ 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 qnum_sector = 8;
+	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;
+	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);
@@ -641,4 +699,4 @@ char* test_mps_to_statevector()
 	H5Fclose(file);
 
 	return 0;
-}
+}

test/run_tests.c

@@ -52,6 +52,7 @@ 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();
@@ -135,6 +136,7 @@ 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),