added Rsvd_notruncate for symmetric and fermionic UniTensors

manuschneider · manuschneider · commit a8ee89cd09e8 · 2026-02-25T20:06:39.000+08:00
diff --git a/include/linalg.hpp b/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.
diff --git a/pybind/linalg_py.cpp b/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(
diff --git a/src/algo/Vstack.cpp b/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++) {
diff --git a/src/linalg/Rsvd.cpp b/src/linalg/Rsvd.cpp
@@ -240,6 +240,151 @@ namespace cytnx {
 
         if (return_err) outCyT.back().Init(outT.back(), false, 0);
       }  // Rsvd_Dense_UT_internal
+
+      void Rsvd_Block_UT_internal(std::vector<UniTensor> &outCyT, const cytnx::UniTensor &Tin,
+                                  const cytnx_uint64 &keepdim, const double &err, const bool &is_U,
+                                  const bool &is_vT, const unsigned int &return_err,
+                                  const cytnx_uint64 &mindim) {
+        cytnx_uint64 keep_dim = keepdim;
+
+        outCyT = linalg::Gesvd(Tin, is_U, is_vT);
+
+        // process truncation:
+        // 1) concate all S vals from all blk
+        Tensor Sall = outCyT[0].get_block_(0);
+        for (int i = 1; i < outCyT[0].Nblocks(); i++) {
+          Sall = algo::Concatenate(Sall, outCyT[0].get_block_(i));
+        }
+        Sall = algo::Sort(Sall);  // all singular values, starting from the smallest
+
+        // 2) get the minimum S value based on the args input.
+        Scalar Smin;
+        cytnx_uint64 smidx;
+        cytnx_uint64 Sshape = Sall.shape()[0];
+        if (keep_dim < Sshape) {
+          smidx = Sshape - keep_dim;
+          Smin = Sall.storage()(smidx);
+        } else {
+          keep_dim = Sshape;
+          smidx = 0;
+          Smin = Sall.storage()(0);
+        }
+        while ((Smin < err) and (keep_dim > (mindim < 1 ? 1 : mindim))) {
+          // at least one singular value is always kept!
+          keep_dim--;
+          // if (keep_dim == 0) break;
+          smidx++;
+          Smin = Sall.storage()(smidx);
+        }
+
+        // traversal each block and truncate!
+        UniTensor &S = outCyT[0];
+        std::vector<cytnx_uint64> new_dims;  // keep_dims for each block!
+        std::vector<cytnx_int64> keep_dims;
+        keep_dims.reserve(S.Nblocks());
+        std::vector<cytnx_int64> new_qid;
+        new_qid.reserve(S.Nblocks());
+
+        std::vector<std::vector<cytnx_uint64>>
+          new_itoi;  // assume S block is in same order as qnum:
+        std::vector<cytnx_uint64> to_be_removed;
+
+        cytnx_uint64 tot_dim = 0;
+        cytnx_uint64 cnt = 0;
+        for (int b = 0; b < S.Nblocks(); b++) {
+          Storage stmp = S.get_block_(b).storage();
+          cytnx_int64 kdim = 0;
+          for (int i = stmp.size(); i > 0; i--) {
+            if (stmp(i - 1) >= Smin) {
+              kdim = i;
+              break;
+            }
+          }
+          keep_dims.push_back(kdim);
+          if (kdim == 0) {
+            to_be_removed.push_back(b);
+            new_qid.push_back(-1);
+
+          } else {
+            new_qid.push_back(new_dims.size());
+            new_itoi.push_back({new_dims.size(), new_dims.size()});
+            new_dims.push_back(kdim);
+            tot_dim += kdim;
+            if (kdim != S.get_blocks_()[b].shape()[0])
+              S.get_blocks_()[b] = S.get_blocks_()[b].get({Accessor::range(0, kdim)});
+          }
+        }
+
+        // remove:
+        // vec_erase_(S.get_itoi(),to_be_removed);
+        S.get_itoi() = new_itoi;
+        vec_erase_(S.get_blocks_(), to_be_removed);
+        vec_erase_(S.bonds()[0].qnums(), to_be_removed);
+        S.bonds()[0]._impl->_degs = new_dims;
+        S.bonds()[0]._impl->_dim = tot_dim;
+        S.bonds()[1] = S.bonds()[0].redirect();
+
+        int t = 1;
+        if (is_U) {
+          UniTensor &U = outCyT[t];
+          to_be_removed.clear();
+          U.bonds().back() = S.bonds()[1].clone();
+          std::vector<Accessor> acs(U.rank());
+          for (int i = 0; i < U.rowrank(); i++) acs[i] = Accessor::all();
+
+          for (int b = 0; b < U.Nblocks(); b++) {
+            if (keep_dims[U.get_qindices(b).back()] == 0)
+              to_be_removed.push_back(b);
+            else {
+              /// process blocks:
+              if (keep_dims[U.get_qindices(b).back()] != U.get_blocks_()[b].shape().back()) {
+                acs.back() = Accessor::range(0, keep_dims[U.get_qindices(b).back()]);
+                U.get_blocks_()[b] = U.get_blocks_()[b].get(acs);
+              }
+
+              // change to new qindices:
+              U.get_qindices(b).back() = new_qid[U.get_qindices(b).back()];
+            }
+          }
+          vec_erase_(U.get_itoi(), to_be_removed);
+          vec_erase_(U.get_blocks_(), to_be_removed);
+
+          t++;
+        }
+
+        if (is_vT) {
+          UniTensor &vT = outCyT[t];
+          to_be_removed.clear();
+          vT.bonds().front() = S.bonds()[0].clone();
+          std::vector<Accessor> acs(vT.rank());
+          for (int i = 1; i < vT.rank(); i++) acs[i] = Accessor::all();
+
+          for (int b = 0; b < vT.Nblocks(); b++) {
+            if (keep_dims[vT.get_qindices(b)[0]] == 0)
+              to_be_removed.push_back(b);
+            else {
+              /// process blocks:
+              if (keep_dims[vT.get_qindices(b)[0]] != vT.get_blocks_()[b].shape()[0]) {
+                acs[0] = Accessor::range(0, keep_dims[vT.get_qindices(b)[0]]);
+                vT.get_blocks_()[b] = vT.get_blocks_()[b].get(acs);
+              }
+              // change to new qindices:
+              vT.get_qindices(b)[0] = new_qid[vT.get_qindices(b)[0]];
+            }
+          }
+          vec_erase_(vT.get_itoi(), to_be_removed);
+          vec_erase_(vT.get_blocks_(), to_be_removed);
+          t++;
+        }
+
+        // handle return_err!
+        if (return_err == 1) {
+          outCyT.push_back(UniTensor(Tensor({1}, Smin.dtype())));
+          outCyT.back().get_block_().storage().at(0) = Smin;
+        } else if (return_err) {
+          outCyT.push_back(UniTensor(Sall.get({Accessor::tilend(smidx)})));
+        }
+      }  // Rsvd_Block_UT_internal
     }  // unnamed namespace
 
     std::vector<cytnx::UniTensor> Rsvd(const cytnx::UniTensor &Tin, cytnx_uint64 keepdim,
@@ -258,18 +403,16 @@ namespace cytnx {
                       "\n");
 
       // check input arguments
-      // cytnx_error_msg(mindim < 0, "[ERROR][Rsvd] mindim must be >=1%s", "\n");
+      cytnx_error_msg(mindim < 0, "[ERROR][Rsvd] mindim must be >=1%s", "\n");
       cytnx_error_msg(keepdim < 1, "[ERROR][Rsvd] keepdim must be >=1%s", "\n");
-      // cytnx_error_msg(return_err < 0, "[ERROR][Rsvd] return_err cannot be negative%s",
-      //                 "\n");
+      cytnx_error_msg(return_err < 0, "[ERROR][Rsvd] return_err cannot be negative%s", "\n");
 
       std::vector<UniTensor> outCyT;
       if (Tin.uten_type() == UTenType.Dense) {
         Rsvd_Dense_UT_internal(outCyT, Tin, keepdim, err, is_U, is_vT, return_err, mindim,
                                oversampling_summand, oversampling_factor, power_iteration, seed);
-        // } else if (Tin.uten_type() == UTenType.Block) {
-        //   _Rsvd_Block_UT(outCyT, Tin, keepdim, err, is_U, is_vT,
-        //   return_err, mindim);
+      } else if (Tin.uten_type() == UTenType.Block) {
+        Rsvd_Block_UT_internal(outCyT, Tin, keepdim, err, is_U, is_vT, return_err, mindim);
       } else {
         cytnx_error_msg(true, "[ERROR][Rsvd] only Dense UniTensors are supported.%s", "\n");
       }
diff --git a/src/linalg/Rsvd_notruncate.cpp b/src/linalg/Rsvd_notruncate.cpp
diff --git a/tests/linalg_test/Rsvd_notruncate_test.cpp b/tests/linalg_test/Rsvd_notruncate_test.cpp
