added RSD for symmetric tensors with and without min_blockdim

Author: manuschneider

include/linalg.hpp

@@ -853,8 +853,8 @@ namespace cytnx {
     // Rsvd:
     //==================================================
     /**
-     * @brief Perform a truncated Singular-Value decomposition of a UniTensor.
-    @details This function will perform a truncated Singular-Value decomposition
+     * @brief Perform a randomized truncated Singular-Value decomposition of a UniTensor.
+    @details This function will perform a randomized truncated Singular-Value decomposition
     of a UniTensor. It uses the ?gesvd method for the SVD. It will perform the randomized
     SVD first, and then truncate the singular values to the given cutoff \p err. That means, given a
     UniTensor \p Tin as \f$ M \f$, then the result will be: \f[ M = U S V^\dagger, \f] where \f$ S
@@ -908,6 +908,41 @@ namespace cytnx {
                                        cytnx_uint64 power_iteration = 0,
                                        unsigned int seed = random::__static_random_device());
 
+    /**
+     * @brief Perform a truncated Singular-Value decomposition of a UniTensor and keep at most
+     * \p min_blockdim singular values in each block.
+     * @details For each block, the minimum dimension can be chosen for the truncation. This can be
+     * helpful to avoid loosing symmetry sectors in the truncated SVD. For more details, please
+     * refer to the documentation of the function \ref Rsvd(const cytnx::UniTensor &Tin,
+     * cytnx_uint64 keepdim, double err = 0., bool is_U = true, bool is_vT = true, unsigned int
+     * return_err = 0, cytnx_uint64 mindim = 1, cytnx_uint64 oversampling_summand = 10, double
+     * oversampling_factor = 1., cytnx_uint64 power_iteration = 0, unsigned int seed =
+     * random::__static_random_device()).
+     *
+     * The truncation order is as following (later constraints might be violated by previous
+     * ones):<br> 1) Keep the largest \p min_blockdim singular values in each block; reduce \p
+     * keepdim and \p mindim by the number of already kept singular values<br> 2) Keep at most \p
+     * keepdim singular values; there might be an exception in case of exact degeneracies where more
+     * singular values are kept<br> 3) Keep at least \p mindim singular values;<br> 4) Drop all
+     * singular values smaller than \p err (no normalization applied to the singular values)
+     *
+     * @param[in] min_blockdim a vector containing the minimum dimension of each block;
+     * alternatively, a vector with only one element can be given to have the same min_blockdim for
+     * each block
+     * @see Rsvd(const cytnx::UniTensor &Tin, cytnx_uint64 keepdim, double err = 0., bool is_U =
+     * true, bool is_vT = true, unsigned int return_err = 0, cytnx_uint64 mindim = 1, cytnx_uint64
+     * oversampling_summand = 10, double oversampling_factor = 1., cytnx_uint64 power_iteration = 0,
+     * unsigned int seed = random::__static_random_device())
+     */
+    std::vector<cytnx::UniTensor> Rsvd(const cytnx::UniTensor &Tin, cytnx_uint64 keepdim,
+                                       const std::vector<cytnx_uint64> min_blockdim,
+                                       double err = 0., bool is_U = true, bool is_vT = true,
+                                       unsigned int return_err = 0, cytnx_uint64 mindim = 1,
+                                       cytnx_uint64 oversampling_summand = 10,
+                                       double oversampling_factor = 1.,
+                                       cytnx_uint64 power_iteration = 0,
+                                       unsigned int seed = random::__static_random_device());
+
     std::vector<cytnx::UniTensor> Hosvd(
       const cytnx::UniTensor &Tin, const std::vector<cytnx_uint64> &mode,
       const bool &is_core = true, const bool &is_Ls = false,

pybind/linalg_py.cpp

@@ -87,7 +87,7 @@ void linalg_binding(py::module &m) {
     py::arg("power_iteration") = 2, py::arg("seed") = -1);
 
   m_linalg.def(
-    "Rsvd",
+    "Rsvd",  // for Tensor
     [](const Tensor &Tin, cytnx_uint64 keepdim, double err, bool is_U, bool is_vT,
        unsigned int return_err, cytnx_uint64 mindim, cytnx_uint64 oversampling_summand,
        double oversampling_factor, cytnx_uint64 power_iteration, int64_t seed) {
@@ -103,7 +103,7 @@ void linalg_binding(py::module &m) {
     py::arg("oversampling_summand") = 20, py::arg("oversampling_factor") = 1.,
     py::arg("power_iteration") = 0, py::arg("seed") = -1);
   m_linalg.def(
-    "Rsvd",
+    "Rsvd",  // for UniTensor, without min_blockdim
     [](const cytnx::UniTensor &Tin, cytnx_uint64 keepdim, double err, bool is_U, bool is_vT,
        unsigned int return_err, cytnx_uint64 mindim, cytnx_uint64 oversampling_summand,
        double oversampling_factor, cytnx_uint64 power_iteration, int64_t seed) {
@@ -118,6 +118,23 @@ void linalg_binding(py::module &m) {
     py::arg("is_vT") = true, py::arg("return_err") = (unsigned int)(0), py::arg("mindim") = 1,
     py::arg("oversampling_summand") = 20, py::arg("oversampling_factor") = 1.,
     py::arg("power_iteration") = 0, py::arg("seed") = -1);
+  m_linalg.def(
+    "Rsvd",  // for UniTensor, with min_blockdim
+    [](const cytnx::UniTensor &Tin, cytnx_uint64 keepdim,
+       const std::vector<cytnx_uint64> min_blockdim, double err, bool is_U, bool is_vT,
+       unsigned int return_err, cytnx_uint64 mindim, cytnx_uint64 oversampling_summand,
+       double oversampling_factor, cytnx_uint64 power_iteration, int64_t seed) {
+      if (seed == -1) {
+        // If user doesn't specify seed argument
+        seed = cytnx::random::__static_random_device();
+      }
+      return cytnx::linalg::Rsvd(Tin, keepdim, min_blockdim, err, is_U, is_vT, return_err, mindim,
+                                 oversampling_summand, oversampling_factor, power_iteration, seed);
+    },
+    py::arg("Tin"), py::arg("keepdim"), py::arg("min_blockdim"), py::arg("err") = double(0),
+    py::arg("is_U") = true, py::arg("is_vT") = true, py::arg("return_err") = (unsigned int)(0),
+    py::arg("mindim") = 1, py::arg("oversampling_summand") = 20,
+    py::arg("oversampling_factor") = 1., py::arg("power_iteration") = 0, py::arg("seed") = -1);
 
   m_linalg.def(
     "Gesvd_truncate",

src/linalg/Gesvd_truncate.cpp

@@ -211,7 +211,7 @@ namespace cytnx {
         }  // if tag
 
         if (return_err) outCyT.back().Init(outT.back(), false, 0);
-      }  // Gesvd_truncate_Dense_UT_internal no minblockdim
+      }  // Gesvd_truncate_Dense_UT_internal
 
       void Gesvd_truncate_Block_UT_internal(std::vector<UniTensor> &outCyT,
                                             const cytnx::UniTensor &Tin,