compression of a MPS by site-local SVDs and truncation, and isometry …

…test for matrices

src/mps/mps.h

@@ -61,21 +61,30 @@ double mps_norm(const struct mps* psi);
 
 // orthonormalization and canonical forms
 
-void mps_local_orthonormalize_qr(struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_next);
-
-void mps_local_orthonormalize_rq(struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_prev);
-
-
 /// \brief MPS orthonormalization mode.
 enum mps_orthonormalization_mode
 {
 	MPS_ORTHONORMAL_LEFT  = 0,  //!< left-orthonormal
 	MPS_ORTHONORMAL_RIGHT = 1,  //!< right-orthonormal
 };
 
+void mps_local_orthonormalize_qr(struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_next);
+
+void mps_local_orthonormalize_rq(struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_prev);
+
 double mps_orthonormalize_qr(struct mps* mps, const enum mps_orthonormalization_mode mode);
 
 
+int mps_local_orthonormalize_left_svd(const double tol, const long max_vdim, const bool renormalize,
+	struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_next, struct trunc_info* info);
+
+int mps_local_orthonormalize_right_svd(const double tol, const long max_vdim, const bool renormalize,
+	struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_prev, struct trunc_info* info);
+
+int mps_compress(const double tol, const long max_vdim, const enum mps_orthonormalization_mode mode,
+	struct mps* mps, double* restrict norm, double* restrict scaling, struct trunc_info* info);
+
+
 //________________________________________________________________________________________________________________________
 //
 

src/tensor/block_sparse_tensor.c

@@ -2312,6 +2312,42 @@ bool block_sparse_tensor_is_identity(const struct block_sparse_tensor* t, const
 }
 
 
+//________________________________________________________________________________________________________________________
+///
+/// \brief Test whether a block-sparse tensors is an isometry within tolerance 'tol'.
+///
+bool block_sparse_tensor_is_isometry(const struct block_sparse_tensor* t, const double tol, const bool transpose)
+{
+	// must be a matrix
+	if (t->ndim != 2) {
+		return false;
+	}
+
+	// entry-wise conjugation
+	struct block_sparse_tensor tc;
+	copy_block_sparse_tensor(t, &tc);
+	conjugate_block_sparse_tensor(&tc);
+	// revert tensor axes directions for multiplication
+	tc.axis_dir[0] = -tc.axis_dir[0];
+	tc.axis_dir[1] = -tc.axis_dir[1];
+
+	struct block_sparse_tensor t2;
+	if (!transpose) {
+		block_sparse_tensor_dot(&tc, TENSOR_AXIS_RANGE_LEADING, t, TENSOR_AXIS_RANGE_LEADING, 1, &t2);
+	}
+	else {
+		block_sparse_tensor_dot(t, TENSOR_AXIS_RANGE_TRAILING, &tc, TENSOR_AXIS_RANGE_TRAILING, 1, &t2);
+	}
+
+	const bool is_isometry = block_sparse_tensor_is_identity(&t2, tol);
+
+	delete_block_sparse_tensor(&t2);
+	delete_block_sparse_tensor(&tc);
+
+	return is_isometry;
+}
+
+
 //________________________________________________________________________________________________________________________
 ///
 /// \brief Overall number of entries in the blocks of the tensor.

src/tensor/block_sparse_tensor.h

@@ -156,6 +156,8 @@ bool block_sparse_tensor_allclose(const struct block_sparse_tensor* restrict s,
 
 bool block_sparse_tensor_is_identity(const struct block_sparse_tensor* t, const double tol);
 
+bool block_sparse_tensor_is_isometry(const struct block_sparse_tensor* t, const double tol, const bool transpose);
+
 
 //________________________________________________________________________________________________________________________
 //

src/tensor/dense_tensor.c

@@ -2745,3 +2745,52 @@ bool dense_tensor_is_identity(const struct dense_tensor* t, const double tol)
 
 	return true;
 }
+
+//________________________________________________________________________________________________________________________
+///
+/// \brief Test whether a dense tensors is an isometry within tolerance 'tol'.
+///
+bool dense_tensor_is_isometry(const struct dense_tensor* t, const double tol, const bool transpose)
+{
+	// must be a matrix
+	if (t->ndim != 2) {
+		return false;
+	}
+
+	bool is_isometry;
+	if (t->dtype == CT_SINGLE_REAL || t->dtype == CT_DOUBLE_REAL)
+	{
+		struct dense_tensor t2;
+		if (!transpose) {
+			dense_tensor_dot(t, TENSOR_AXIS_RANGE_LEADING, t, TENSOR_AXIS_RANGE_LEADING, 1, &t2);
+		}
+		else {
+			dense_tensor_dot(t, TENSOR_AXIS_RANGE_TRAILING, t, TENSOR_AXIS_RANGE_TRAILING, 1, &t2);
+		}
+
+		is_isometry = dense_tensor_is_identity(&t2, tol);
+
+		delete_dense_tensor(&t2);
+	}
+	else
+	{
+		struct dense_tensor tc;
+		copy_dense_tensor(t, &tc);
+		conjugate_dense_tensor(&tc);
+
+		struct dense_tensor t2;
+		if (!transpose) {
+			dense_tensor_dot(&tc, TENSOR_AXIS_RANGE_LEADING, t, TENSOR_AXIS_RANGE_LEADING, 1, &t2);
+		}
+		else {
+			dense_tensor_dot(t, TENSOR_AXIS_RANGE_TRAILING, &tc, TENSOR_AXIS_RANGE_TRAILING, 1, &t2);
+		}
+
+		is_isometry = dense_tensor_is_identity(&t2, tol);
+
+		delete_dense_tensor(&t2);
+		delete_dense_tensor(&tc);
+	}
+
+	return is_isometry;
+}

src/tensor/dense_tensor.h

@@ -216,3 +216,5 @@ void dense_tensor_block(const struct dense_tensor* restrict t, const long* restr
 bool dense_tensor_allclose(const struct dense_tensor* restrict s, const struct dense_tensor* restrict t, const double tol);
 
 bool dense_tensor_is_identity(const struct dense_tensor* t, const double tol);
+
+bool dense_tensor_is_isometry(const struct dense_tensor* t, const double tol, const bool transpose);

test/algorithm/test_bond_ops.c

@@ -165,36 +165,16 @@ char* test_split_block_sparse_matrix_svd()
 			if (d == 0)
 			{
 				// a1 must be an isometry
-				struct block_sparse_tensor a1c;
-				copy_block_sparse_tensor(&a1, &a1c);
-				conjugate_block_sparse_tensor(&a1c);
-				// revert tensor axes directions for multiplication
-				a1c.axis_dir[0] = -a1c.axis_dir[0];
-				a1c.axis_dir[1] = -a1c.axis_dir[1];
-				struct block_sparse_tensor p;
-				block_sparse_tensor_dot(&a1, TENSOR_AXIS_RANGE_TRAILING, &a1c, TENSOR_AXIS_RANGE_TRAILING, 1, &p);
-				if (!block_sparse_tensor_is_identity(&p, 5e-6)) {
+				if (!block_sparse_tensor_is_isometry(&a1, 5e-6, true)) {
 					return "'a1' matrix is not an isometry";
 				}
-				delete_block_sparse_tensor(&p);
-				delete_block_sparse_tensor(&a1c);
 			}
 			else
 			{
 				// a0 must be an isometry
-				struct block_sparse_tensor a0c;
-				copy_block_sparse_tensor(&a0, &a0c);
-				conjugate_block_sparse_tensor(&a0c);
-				// revert tensor axes directions for multiplication
-				a0c.axis_dir[0] = -a0c.axis_dir[0];
-				a0c.axis_dir[1] = -a0c.axis_dir[1];
-				struct block_sparse_tensor p;
-				block_sparse_tensor_dot(&a0c, TENSOR_AXIS_RANGE_LEADING, &a0, TENSOR_AXIS_RANGE_LEADING, 1, &p);
-				if (!block_sparse_tensor_is_identity(&p, 5e-6)) {
-					return "'a1' matrix is not an isometry";
+				if (!block_sparse_tensor_is_isometry(&a0, 5e-6, false)) {
+					return "'a0' matrix is not an isometry";
 				}
-				delete_block_sparse_tensor(&p);
-				delete_block_sparse_tensor(&a0c);
 			}
 
 			// reassemble matrix after splitting, for comparison with reference
@@ -277,19 +257,9 @@ char* test_split_block_sparse_matrix_svd_zero()
 		if (d == 1)
 		{
 			// a0 must be an isometry
-			struct block_sparse_tensor a0c;
-			copy_block_sparse_tensor(&a0, &a0c);
-			conjugate_block_sparse_tensor(&a0c);
-			// revert tensor axes directions for multiplication
-			a0c.axis_dir[0] = -a0c.axis_dir[0];
-			a0c.axis_dir[1] = -a0c.axis_dir[1];
-			struct block_sparse_tensor p;
-			block_sparse_tensor_dot(&a0c, TENSOR_AXIS_RANGE_LEADING, &a0, TENSOR_AXIS_RANGE_LEADING, 1, &p);
-			if (!block_sparse_tensor_is_identity(&p, 5e-6)) {
-				return "'a1' matrix is not an isometry";
+			if (!block_sparse_tensor_is_isometry(&a0, 5e-6, false)) {
+				return "'a0' matrix is not an isometry";
 			}
-			delete_block_sparse_tensor(&p);
-			delete_block_sparse_tensor(&a0c);
 		}
 
 		// reassemble matrix after splitting, to ensure dimension compatibilty