added Rsvd_notruncate for symmetric and fermionic UniTensors

Author: manuschneider

include/linalg.hpp

@@ -708,19 +708,30 @@ namespace cytnx {
     subspace, and the smallest singular values are less reliable. Use truncation with oversampling.
     @param[in] is_U if \em true, the left-unitary UniTensor U (isometry) is returned.
     @param[in] is_vT if \em true, the right-unitary UniTensor vT (isometry) is returned.
+    @param[in] mindim at least this amount of singular values are kept in each block.
+    @param[in] oversampling_summand the randomized SVD computes [(1 + oversampling_factor) *
+    keepdim*d/D + oversampling_summand] singular values in each block before further truncating,
+    where d is the block dimension and D the full tensor dimension (each being the minimum of row
+    and column dimension of the block or tensor respectively).
+    @param[in] oversampling_factor see oversampling_summand
     @param[in] power_iteration number of iterations for the power method: Y = (A *
     Adag)^power_iteration * A * Tin
     @param[in] seed the seed for the random generator. [Default] Using device entropy.
     @see Rand_isometry(const Tensor &Tin, const cytnx_uint64 &keepdim, const cytnx_uint64
     &power_iteration, const unsigned int &seed)
     @see Rsvd()
     @see Gesvd()
+    @note At least one singular value per symmetry sector is always kept.
+    @note More than keepdim singular values might be returned (depending on mindim,
+    oversampling_summand, oversampling_factor, and symmetry sector sizes).
     @warning No truncation of the singular values is performed, and the smaller ones will be
     inaccurate. Use Rsvd() for a truncated version that drops small singular values.
     */
     std::vector<cytnx::UniTensor> Rsvd_notruncate(
       const cytnx::UniTensor &Tin, cytnx_uint64 keepdim, bool is_U = true, bool is_vT = true,
-      cytnx_uint64 power_iteration = 2, unsigned int seed = random::__static_random_device());
+      cytnx_uint64 mindim = 1, cytnx_uint64 oversampling_summand = 0,
+      double oversampling_factor = 0., cytnx_uint64 power_iteration = 2,
+      unsigned int seed = random::__static_random_device());
 
     /**
     @brief Perform Singular-Value decomposition on a UniTensor using ?gesvd method.
@@ -865,10 +876,12 @@ namespace cytnx {
     the largest error will be pushed back to the vector (The smallest singular value in the return
     singular values matrix \f$ S \f$.) If \p return_err is > 1, then the full list of truncated
     singular values will be returned.
-    @param[in] oversampling_summand the randomized SVD computes (1 + oversampling_fact) * keepdim +
-    oversampling_summand before further truncating to keepdim
-    @param[in] oversampling_fact the randomized SVD computes (1 + oversampling_fact) * keepdim +
-    oversampling_summand before further truncating to keepdim
+    @param[in] mindim at least this amount of singular values are kept in each block.
+    @param[in] oversampling_summand the randomized SVD computes [(1 + oversampling_factor) *
+    keepdim*d/D + oversampling_summand] singular values in each block before further truncating,
+    where d is the block dimension and D the full tensor dimension (each being the minimum of row
+    and column dimension of the block or tensor respectively).
+    @param[in] oversampling_factor see oversampling_summand
     @param[in] power_iteration number of iterations for the power method: Y = (A *
     Adag)^power_iteration * A * Tin
     @param[in] seed the seed for the random generator. [Default] Using device entropy.
@@ -1887,10 +1900,12 @@ namespace cytnx {
     the largest error will be pushed back to the vector (The smallest singular value in the return
     singular values matrix \f$ S \f$.) If \p return_err is a \em positive int, then the
     full list of truncated singular values will be returned.
-    @param[in] oversampling_summand the randomized SVD computes (1 + oversampling_fact) * keepdim +
-    oversampling_summand before further truncating to keepdim
-    @param[in] oversampling_fact the randomized SVD computes (1 + oversampling_fact) * keepdim +
-    oversampling_summand before further truncating to keepdim
+    @param[in] mindim at least this amount of singular values are kept in each block.
+    @param[in] oversampling_summand the randomized SVD computes [(1 + oversampling_factor) *
+    keepdim*d/D + oversampling_summand] singular values in each block before further truncating,
+    where d is the block dimension and D the full tensor dimension (each being the minimum of row
+    and column dimension of the block or tensor respectively).
+    @param[in] oversampling_factor see oversampling_summand
     @param[in] power_iteration number of iterations for the power method: Y = (A *
     Adag)^power_iteration * A * Tin
     @param[in] seed the seed for the random generator. [Default] Using device entropy.

pybind/linalg_py.cpp

@@ -73,14 +73,17 @@ void linalg_binding(py::module &m) {
   m_linalg.def(
     "Rsvd_notruncate",
     [](const cytnx::UniTensor &Tin, cytnx_uint64 keepdim, bool is_U, bool is_vT,
+       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_notruncate(Tin, keepdim, is_U, is_vT, power_iteration, seed);
+      return cytnx::linalg::Rsvd_notruncate(Tin, keepdim, is_U, is_vT, mindim, oversampling_summand,
+                                            oversampling_factor, power_iteration, seed);
     },
     py::arg("Tin"), py::arg("keepdim"), py::arg("is_U") = true, py::arg("is_vT") = true,
+    py::arg("mindim") = 1, py::arg("oversampling_summand") = 0, py::arg("oversampling_factor") = 0.,
     py::arg("power_iteration") = 2, py::arg("seed") = -1);
 
   m_linalg.def(

src/algo/Vstack.cpp

@@ -16,8 +16,6 @@ namespace cytnx {
     typedef Accessor ac;
 
     Tensor Vstack(const std::vector<Tensor> &In_tensors) {
-      Tensor out;
-
       std::vector<Tensor> _Ins;
 
       // check:
@@ -82,7 +80,7 @@ namespace cytnx {
       }
 
       // allocate out!
-      out = zeros({Dcomb, Dshare}, dtype_id, device_id);
+      Tensor out = zeros({Dcomb, Dshare}, dtype_id, device_id);
 
       std::vector<void *> rawPtr(In_tensors.size());
       for (int i = 0; i < _Ins.size(); i++) {