UniTensor.hpp

Go to the documentation of this file.

#ifndef CYTNX_UNITENSOR_H_
#define CYTNX_UNITENSOR_H_
 
#include "Type.hpp"
#include "cytnx_error.hpp"
#include "Device.hpp"
#include "Tensor.hpp"
#include "utils/utils.hpp"
#include "intrusive_ptr_base.hpp"
#include <iostream>
#include <vector>
#include <map>
#include <utility>
#include <initializer_list>
#include <fstream>
#include <algorithm>
#include "Symmetry.hpp"
#include "Bond.hpp"
#include "Generator.hpp"
#include <random>
 
#ifdef BACKEND_TORCH
#else
  #include "backend/Scalar.hpp"
 
// namespace cytnx{
namespace cytnx {
  namespace random {
    extern std::random_device __static_random_device;
  }
 
  class UniTensorType_class {
   public:
    enum : int { Void = -99, Dense = 0, Sparse = 1, Block = 2, BlockFermionic = 3 };
    std::string getname(const int &ut_type) const;
  };
 
  extern UniTensorType_class UTenType;
 
  // class DenseUniTensor;
  class UniTensor_base : public intrusive_ptr_base<UniTensor_base> {
   public:
    int uten_type_id;  // the unitensor type id.
    bool _is_braket_form;
    bool _is_tag;
    bool _is_diag;
    cytnx_int64 _rowrank;
    std::string _name;
    std::vector<std::string> _labels;
    std::vector<Bond> _bonds;
    std::vector<Symmetry> syms_cache;
 
    bool _update_braket() {
      if (_bonds.size() == 0) return false;
 
      if (this->_bonds[0].type() != bondType::BD_REG) {
        // check:
        for (unsigned int i = 0; i < this->_bonds.size(); i++) {
          if (i < this->_rowrank) {
            if (this->_bonds[i].type() != bondType::BD_KET) return false;
          } else {
            if (this->_bonds[i].type() != bondType::BD_BRA) return false;
          }
        }
        return true;
      } else {
        return false;
      }
    }
 
    friend class UniTensor;  // allow wrapper to access the private elems
    friend class DenseUniTensor;
    friend class BlockUniTensor;
    friend class BlockFermionicUniTensor;
 
    UniTensor_base()
        : _is_tag(false),
          _name(std::string("")),
          _is_braket_form(false),
          _rowrank(0),
          _is_diag(false),
          uten_type_id(UTenType.Void){};
 
    // copy&assignment constr., use intrusive_ptr's !!
    UniTensor_base(const UniTensor_base &rhs);
    UniTensor_base &operator=(UniTensor_base &rhs);
 
    cytnx_uint64 rowrank() const { return this->_rowrank; }
    bool is_diag() const { return this->_is_diag; }
    const bool &is_braket_form() const { return this->_is_braket_form; }
    const bool &is_tag() const { return this->_is_tag; }
    const std::vector<std::string> &labels() const { return this->_labels; }
    cytnx_int64 get_index(std::string label) const {
      std::vector<std::string> labels = this->_labels;
      for (cytnx_uint64 i = 0; i < labels.size(); i++) {
        if (labels[i] == label) return i;
      }
      return -1;
    }
    const std::vector<Bond> &bonds() const { return this->_bonds; }
 
    Bond &bond_(const cytnx_uint64 &idx) {
      cytnx_error_msg(idx >= this->rank(), "[ERROR][bond] index %d out of bound, total %d\n", idx,
                      this->rank());
      return this->_bonds[idx];
    }
 
    Bond &bond_(const std::string &label) {
      auto res = std::find(this->_labels.begin(), this->_labels.end(), label);
      cytnx_error_msg(res == this->_labels.end(), "[ERROR] label %s not exists.\n", label.c_str());
      cytnx_uint64 idx = std::distance(this->_labels.begin(), res);
 
      return this->bond_(idx);
    }
 
    const std::string &name() const { return this->_name; }
    cytnx_uint64 rank() const {
      cytnx_error_msg(this->_labels.size() != this->_bonds.size(),
                      "[ERROR][UniTensor_base][rank] inconsistent metadata: labels.size()=%zu, "
                      "bonds.size()=%zu.%s",
                      this->_labels.size(), this->_bonds.size(), "\n");
      return this->_labels.size();
    }
    void set_name_(const std::string &in) { this->_name = in; }
    [[deprecated("Please use set_name_(const std::string &in) instead.")]] void set_name(
      const std::string &in) {
      this->set_name_(in);
    }
 
    void set_label_(const std::string &oldlabel, const std::string &new_label) {
      cytnx_int64 idx;
      auto res = std::find(this->_labels.begin(), this->_labels.end(), oldlabel);
      cytnx_error_msg(res == this->_labels.end(), "[ERROR] label %s not exists.\n",
                      oldlabel.c_str());
      idx = std::distance(this->_labels.begin(), res);
 
      cytnx_error_msg(idx >= this->_labels.size(), "[ERROR] index exceed the rank of UniTensor%s",
                      "\n");
      // check in:
      bool is_dup = false;
      for (cytnx_uint64 i = 0; i < this->_labels.size(); i++) {
        if (i == idx) continue;
        if (new_label == this->_labels[i]) {
          is_dup = true;
          break;
        }
      }
      cytnx_error_msg(is_dup, "[ERROR] alreay has a label that is the same as the input label%s",
                      "\n");
      this->_labels[idx] = new_label;
    }
    void set_label_(const cytnx_int64 &inx, const std::string &new_label) {
      cytnx_error_msg(inx < 0, "[ERROR] index is negative%s", "\n");
      cytnx_error_msg(inx >= this->_labels.size(), "[ERROR] index exceed the rank of UniTensor%s",
                      "\n");
      // check in:
      bool is_dup = false;
      for (cytnx_uint64 i = 0; i < this->_labels.size(); i++) {
        if (i == inx) continue;
        if (new_label == this->_labels[i]) {
          is_dup = true;
          break;
        }
      }
      cytnx_error_msg(is_dup, "[ERROR] alreay has a label that is the same as the input label%s",
                      "\n");
      this->_labels[inx] = new_label;
    }
    [[deprecated(
      "Please use set_label_(const std::string &oldlabel, const std::string "
      "&new_label) instead.")]] void
      set_label(const std::string &oldlabel, const std::string &new_label) {
      this->set_label_(oldlabel, new_label);
    }
    [[deprecated(
      "Please use set_label_(const cytnx_int64 &inx, const std::string &new_label) "
      "instead.")]] void
      set_label(const cytnx_int64 &inx, const std::string &new_label) {
      this->set_label_(inx, new_label);
    }
 
    [[deprecated("Please use relabel_(const std::vector<std::string> &new_labels) instead.")]] void
      set_labels(const std::vector<std::string> &new_labels);
    void relabel_(const std::vector<std::string> &new_labels);  // implemented
    [[deprecated("Please use relabel_(const std::vector<std::string> &new_labels) instead.")]] void
      relabels_(const std::vector<std::string> &new_labels);  // implemented
    void relabel_(const std::vector<std::string> &old_labels,
                  const std::vector<std::string> &new_labels);  // implemented
    [[deprecated(
      "Please use relabel_(const std::vector<std::string> &old_labels, const "
      "std::vector<std::string> &new_labels) instead.")]] void
      relabels_(const std::vector<std::string> &old_labels,
                const std::vector<std::string> &new_labels);  // implemented
    void relabel_(const std::string &old_label, const std::string &new_label) {
      this->set_label_(old_label, new_label);
    }
    void relabel_(const cytnx_int64 &inx, const std::string &new_label) {
      this->set_label_(inx, new_label);
    }
 
    int uten_type() { return this->uten_type_id; }
    std::string uten_type_str() const { return UTenType.getname(this->uten_type_id); }
 
 
    // string labels!
    virtual void Init(const std::vector<Bond> &bonds,
                      const std::vector<std::string> &in_labels = {},
                      const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
                      const int &device = Device.cpu, const bool &is_diag = false,
                      const bool &no_alloc = false, const std::string &name = "");
 
    virtual void Init_by_Tensor(const Tensor &in, const bool &is_diag = false,
                                const cytnx_int64 &rowrank = -1, const std::string &name = "");
    virtual std::vector<cytnx_uint64> shape() const;
    virtual std::vector<bool> signflip() const;
    virtual bool is_blockform() const;
    virtual bool is_contiguous() const;
    virtual void to_(const int &device);
    virtual boost::intrusive_ptr<UniTensor_base> to(const int &device);
    virtual boost::intrusive_ptr<UniTensor_base> clone() const;
    virtual unsigned int dtype() const;
    virtual int device() const;
    virtual std::string dtype_str() const;
    virtual std::string device_str() const;
    virtual void set_rowrank_(const cytnx_uint64 &new_rowrank);
    virtual boost::intrusive_ptr<UniTensor_base> set_rowrank(const cytnx_uint64 &new_rowrank) const;
 
    virtual boost::intrusive_ptr<UniTensor_base> permute(const std::vector<cytnx_int64> &mapper,
                                                         const cytnx_int64 &rowrank = -1);
    virtual boost::intrusive_ptr<UniTensor_base> permute(const std::vector<std::string> &mapper,
                                                         const cytnx_int64 &rowrank = -1);
    // virtual boost::intrusive_ptr<UniTensor_base> permute(const std::vector<cytnx_int64> &mapper,
    //                                                      const cytnx_int64 &rowrank = -1);
 
    virtual void permute_(const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1);
    virtual void permute_(const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1);
 
    virtual boost::intrusive_ptr<UniTensor_base> permute_nosignflip(
      const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1);
    virtual boost::intrusive_ptr<UniTensor_base> permute_nosignflip(
      const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1);
    virtual void permute_nosignflip_(const std::vector<cytnx_int64> &mapper,
                                     const cytnx_int64 &rowrank = -1);
    virtual void permute_nosignflip_(const std::vector<std::string> &mapper,
                                     const cytnx_int64 &rowrank = -1);
    // virtual void permute_(const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank =
    // -1);
 
    virtual void twist_(const cytnx_int64 &idx);
    virtual void twist_(const std::string &label);
    virtual void fermion_twists_();
 
    virtual boost::intrusive_ptr<UniTensor_base> contiguous_();
    virtual boost::intrusive_ptr<UniTensor_base> contiguous();
    virtual boost::intrusive_ptr<UniTensor_base> apply_();
    virtual boost::intrusive_ptr<UniTensor_base> apply();
    virtual void print_diagram(const bool &bond_info = false) const;
    virtual void print_blocks(const bool &full_info = true) const;
    virtual void print_block(const cytnx_int64 &idx, const bool &full_info = true) const;
 
    virtual boost::intrusive_ptr<UniTensor_base> astype(const unsigned int &dtype) const;
 
    virtual cytnx_uint64 Nblocks() const { return 0; };
    virtual Tensor get_block(const cytnx_uint64 &idx = 0) const;  // return a copy of block
    virtual Tensor get_block(const std::vector<cytnx_int64> &qnum,
                             const bool &force) const;  // return a copy of block
 
    virtual const Tensor &get_block_(const cytnx_uint64 &idx = 0)
      const;  // return a share view of block, this only work for non-symm tensor.
    virtual const Tensor &get_block_(const std::vector<cytnx_int64> &qnum,
                                     const bool &force) const;  // return a copy of block
    virtual Tensor &get_block_(const cytnx_uint64 &idx = 0);  // return a share view of block, this
                                                              // only work for non-symm tensor.
    virtual Tensor &get_block_(const std::vector<cytnx_int64> &qnum,
                               const bool &force);  // return a copy of block
    virtual bool same_data(const boost::intrusive_ptr<UniTensor_base> &rhs) const;
 
    virtual std::vector<Tensor> get_blocks() const;
    virtual const std::vector<Tensor> &get_blocks_(const bool &) const;
    virtual std::vector<Tensor> &get_blocks_(const bool &);
 
    virtual void put_block(const Tensor &in, const cytnx_uint64 &idx = 0);
    virtual void put_block_(Tensor &in, const cytnx_uint64 &idx = 0);
    virtual void put_block(const Tensor &in, const std::vector<cytnx_int64> &qnum);
    virtual void put_block_(Tensor &in, const std::vector<cytnx_int64> &qnum);
 
    // this will only work on non-symm tensor (DenseUniTensor)
    virtual boost::intrusive_ptr<UniTensor_base> get(const std::vector<Accessor> &accessors);
 
    // this will only work on non-symm tensor (DenseUniTensor)
    virtual void set(const std::vector<Accessor> &accessors, const Tensor &rhs);
 
    virtual void reshape_(const std::vector<cytnx_int64> &new_shape,
                          const cytnx_uint64 &rowrank = 0);
    virtual boost::intrusive_ptr<UniTensor_base> reshape(const std::vector<cytnx_int64> &new_shape,
                                                         const cytnx_uint64 &rowrank = 0);
    virtual boost::intrusive_ptr<UniTensor_base> to_dense();
    virtual void to_dense_();
    virtual void combineBond(const std::vector<std::string> &indicators, const bool &force = false);
    virtual void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force,
                              const bool &by_label);
    virtual void combineBonds(const std::vector<std::string> &indicators,
                              const bool &force = false);
    virtual void combineBonds(const std::vector<cytnx_int64> &indicators,
                              const bool &force = false);
    virtual boost::intrusive_ptr<UniTensor_base> contract(
      const boost::intrusive_ptr<UniTensor_base> &rhs, const bool &mv_elem_self = false,
      const bool &mv_elem_rhs = false);
    virtual std::vector<Bond> getTotalQnums(const bool &physical = false);
    virtual std::vector<std::vector<cytnx_int64>> get_blocks_qnums() const;
    virtual void Trace_(const std::string &a, const std::string &b);
    virtual void Trace_(const cytnx_int64 &a, const cytnx_int64 &b);
 
    virtual boost::intrusive_ptr<UniTensor_base> Trace(const std::string &a, const std::string &b);
    virtual boost::intrusive_ptr<UniTensor_base> Trace(const cytnx_int64 &a, const cytnx_int64 &b);
 
    virtual boost::intrusive_ptr<UniTensor_base> relabel(
      const std::vector<std::string> &new_labels);
    virtual boost::intrusive_ptr<UniTensor_base> relabels(
      const std::vector<std::string> &new_labels);
 
    virtual boost::intrusive_ptr<UniTensor_base> relabel(
      const std::vector<std::string> &old_labels, const std::vector<std::string> &new_labels);
    virtual boost::intrusive_ptr<UniTensor_base> relabels(
      const std::vector<std::string> &old_labels, const std::vector<std::string> &new_labels);
 
    virtual boost::intrusive_ptr<UniTensor_base> relabel(const std::string &old_label,
                                                         const std::string &new_label);
 
    virtual boost::intrusive_ptr<UniTensor_base> relabel(const cytnx_int64 &inx,
                                                         const std::string &new_label);
 
    virtual std::vector<Symmetry> syms() const;
 
    // arithmetic
    // internal/advanced: the UniTensor-vs-UniTensor overloads below (Add_/Mul_/Sub_/Div_
    // taking a boost::intrusive_ptr<UniTensor_base>) are not TN operations; python surface
    // removed per #934. See operator+(const UniTensor&, const UniTensor&) in
    // include/linalg.hpp for the full rationale. The Scalar overloads (scalar scaling)
    // remain fully public, including from python.
    virtual void Add_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    virtual void Add_(const Scalar &rhs);
 
    virtual void Mul_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    virtual void Mul_(const Scalar &rhs);
 
    virtual void Sub_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    virtual void Sub_(const Scalar &rhs);
    virtual void lSub_(const Scalar &lhs);
 
    virtual void Div_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    virtual void Div_(const Scalar &rhs);
    virtual void lDiv_(const Scalar &lhs);
 
    virtual Tensor Norm() const;
    virtual boost::intrusive_ptr<UniTensor_base> normalize();
    virtual void normalize_();
 
    virtual boost::intrusive_ptr<UniTensor_base> Conj();
    virtual void Conj_();
 
    virtual boost::intrusive_ptr<UniTensor_base> Transpose();
    virtual void Transpose_();
 
    virtual boost::intrusive_ptr<UniTensor_base> Dagger();
    virtual void Dagger_();
 
    virtual void tag_();
 
    virtual void truncate_(const std::string &label, const cytnx_uint64 &dim);
    virtual void truncate_(const cytnx_int64 &bond_idx, const cytnx_uint64 &dim);
 
    virtual bool elem_exists(const std::vector<cytnx_uint64> &locator) const;
 
    // this a workaround, as virtual function cannot template.
    virtual Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator);
    virtual const Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator) const;
 
    virtual cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                            const cytnx_complex128 &aux);
    virtual cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                           const cytnx_complex64 &aux);
    virtual cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                        const cytnx_double &aux);
    virtual cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                       const cytnx_float &aux);
    virtual cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                        const cytnx_uint64 &aux);
    virtual cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                       const cytnx_int64 &aux);
    virtual cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                        const cytnx_uint32 &aux);
    virtual cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                       const cytnx_int32 &aux);
    virtual cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                        const cytnx_uint16 &aux);
    virtual cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                       const cytnx_int16 &aux);
 
    virtual const cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                                  const cytnx_complex128 &aux) const;
    virtual const cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                                 const cytnx_complex64 &aux) const;
    virtual const cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                              const cytnx_double &aux) const;
    virtual const cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                             const cytnx_float &aux) const;
    virtual const cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                              const cytnx_uint64 &aux) const;
    virtual const cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                             const cytnx_int64 &aux) const;
    virtual const cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                              const cytnx_uint32 &aux) const;
    virtual const cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                             const cytnx_int32 &aux) const;
    virtual const cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                              const cytnx_uint16 &aux) const;
    virtual const cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                             const cytnx_int16 &aux) const;
 
    virtual void from_(const boost::intrusive_ptr<UniTensor_base> &rhs, bool force,
                       cytnx_double tol = 0.);
 
    virtual void group_basis_();
    virtual const std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) const;
    virtual std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx);
    virtual const vec2d<cytnx_uint64> &get_itoi() const;
    virtual vec2d<cytnx_uint64> &get_itoi();
 
    virtual void _save_dispatch(std::fstream &f) const;
    virtual void _load_dispatch(std::fstream &f, unsigned int version);
 
    virtual ~UniTensor_base(){};
  };
 
  //======================================================================
  class DenseUniTensor : public UniTensor_base {
   protected:
   public:
    Tensor _block;
    std::vector<Tensor> _interface_block;  // this is serves as interface for get_blocks_();
    boost::intrusive_ptr<DenseUniTensor> clone_meta() const {
      boost::intrusive_ptr<DenseUniTensor> tmp(new DenseUniTensor());
      tmp->_bonds = vec_clone(this->_bonds);
      tmp->_labels = this->_labels;
      tmp->syms_cache = vec_clone(this->syms_cache);
      tmp->_is_braket_form = this->_is_braket_form;
      tmp->_rowrank = this->_rowrank;
      tmp->_is_diag = this->_is_diag;
      tmp->_name = this->_name;
      tmp->_is_tag = this->_is_tag;
      return tmp;
    }
    //------------------------------------------
 
    DenseUniTensor() { this->uten_type_id = UTenType.Dense; };
    friend class UniTensor;  // allow wrapper to access the private elems
    // virtual functions
 
    // void Init(const std::vector<Bond> &bonds, const std::vector<cytnx_int64> &in_labels = {},
    //           const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
    //           const int &device = Device.cpu, const bool &is_diag = false,
    //           const bool &no_alloc = false);
 
    void Init(const std::vector<Bond> &bonds, const std::vector<std::string> &in_labels = {},
              const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
              const int &device = Device.cpu, const bool &is_diag = false,
              const bool &no_alloc = false, const std::string &name = "");
    // this only work for non-symm tensor
    void Init_by_Tensor(const Tensor &in_tensor, const bool &is_diag = false,
                        const cytnx_int64 &rowrank = -1, const std::string &name = "");
    std::vector<cytnx_uint64> shape() const {
      if (this->_is_diag) {
        std::vector<cytnx_uint64> shape = this->_block.shape();
        shape.push_back(shape[0]);
        return shape;
      } else {
        return this->_block.shape();
      }
    }
    bool is_blockform() const { return false; }
    void to_(const int &device) { this->_block.to_(device); }
    boost::intrusive_ptr<UniTensor_base> to(const int &device) {
      if (this->device() == device) {
        std::vector<Tensor> _interface_block;  // this is serves as interface for get_blocks_();
        return this;
      } else {
        boost::intrusive_ptr<UniTensor_base> out = this->clone();
        out->to_(device);
        return out;
      }
    }
    void set_rowrank_(const cytnx_uint64 &new_rowrank) {
      cytnx_error_msg(new_rowrank > this->_labels.size(),
                      "[ERROR] rowrank cannot exceed the rank of UniTensor.%s", "\n");
      if (this->is_diag()) {
        cytnx_error_msg(new_rowrank != 1, "[ERROR] rowrank should be [==1] when is_diag =true!.%s",
                        "\n");
      }
 
      this->_rowrank = new_rowrank;
    }
 
    boost::intrusive_ptr<UniTensor_base> set_rowrank(const cytnx_uint64 &new_rowrank) const {
      boost::intrusive_ptr<DenseUniTensor> out_raw = this->clone_meta();
      out_raw->_block = this->_block;
      out_raw->set_rowrank_(new_rowrank);
      return out_raw;
    }
 
    boost::intrusive_ptr<UniTensor_base> clone() const {
      boost::intrusive_ptr<DenseUniTensor> tmp = this->clone_meta();
      tmp->_block = this->_block.clone();
      return tmp;
    };
    bool is_contiguous() const { return this->_block.is_contiguous(); }
    unsigned int dtype() const { return this->_block.dtype(); }
    int device() const { return this->_block.device(); }
    std::string dtype_str() const { return Type.getname(this->_block.dtype()); }
    std::string device_str() const { return Device.getname(this->_block.device()); }
    boost::intrusive_ptr<UniTensor_base> permute(const std::vector<cytnx_int64> &mapper,
                                                 const cytnx_int64 &rowrank = -1);
    boost::intrusive_ptr<UniTensor_base> permute(const std::vector<std::string> &mapper,
                                                 const cytnx_int64 &rowrank = -1);
 
    void permute_(const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1);
    void permute_(const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1);
 
    void twist_(const cytnx_int64 &idx) override {
      // do nothing for bosonic UniTensor
      return;
    }
    void twist_(const std::string &label) override {
      // do nothing for bosonic UniTensor
      return;
    }
 
    void fermion_twists_() override {
      // do nothing for bosonic UniTensor
      return;
    }
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::vector<std::string> &new_labels);
    boost::intrusive_ptr<UniTensor_base> relabels(const std::vector<std::string> &new_labels);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::vector<std::string> &old_labels,
                                                 const std::vector<std::string> &new_labels);
    boost::intrusive_ptr<UniTensor_base> relabels(const std::vector<std::string> &old_labels,
                                                  const std::vector<std::string> &new_labels);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::string &old_label,
                                                 const std::string &new_label);
    boost::intrusive_ptr<UniTensor_base> relabel(const cytnx_int64 &inx,
                                                 const std::string &new_label);
 
    boost::intrusive_ptr<UniTensor_base> astype(const unsigned int &dtype) const {
      boost::intrusive_ptr<DenseUniTensor> tmp = this->clone_meta();
      tmp->_block = this->_block.astype(dtype);
      return tmp;
    }
 
    std::vector<Symmetry> syms() const {
      cytnx_error_msg(true, "[ERROR][DenseUniTensor] dense unitensor does not have symmetry.%s",
                      "\n");
      return std::vector<Symmetry>();
    }
 
    boost::intrusive_ptr<UniTensor_base> contiguous_() {
      // Rebind _block to the non-mutating Tensor::contiguous() result instead of calling
      // contiguous_() in place, so a block Tensor shared with another UniTensor (e.g. via
      // relabel()) is not corrupted (#724). Tensor::contiguous() is cheap when the block is
      // already contiguous (no data copy; same impl returned).
      this->_block = this->_block.contiguous();
      return boost::intrusive_ptr<UniTensor_base>(this);
    }
    boost::intrusive_ptr<UniTensor_base> contiguous() {
      // if contiguous then return self!
      if (this->is_contiguous()) {
        boost::intrusive_ptr<UniTensor_base> out(this);
        return out;
      } else {
        boost::intrusive_ptr<DenseUniTensor> tmp = this->clone_meta();
        tmp->_block = this->_block.contiguous();
        return tmp;
      }
    }
 
    boost::intrusive_ptr<UniTensor_base> apply_() {
      return boost::intrusive_ptr<UniTensor_base>(this);
    }
    boost::intrusive_ptr<UniTensor_base> apply() {
      // just return self
      boost::intrusive_ptr<UniTensor_base> out(this);
      return out;
    }
 
    void print_diagram(const bool &bond_info = false) const;
    void print_blocks(const bool &full_info = true) const;
    void print_block(const cytnx_int64 &idx, const bool &full_info = true) const;
    Tensor get_block() const { return this->_block.clone(); }
    Tensor get_block(const cytnx_uint64 &idx) const {
      cytnx_error_msg(idx != 0,
                      "[ERROR][DenseUniTensor] Dense tensor has only one block, block number %llu "
                      "invalid. Use get_block(0).\n",
                      (unsigned long long)idx);
      return this->_block.clone();
    }
 
    Tensor get_block(const std::vector<cytnx_int64> &qnum, const bool &force) const {
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] try to get_block() using qnum on a non-symmetry UniTensor%s",
        "\n");
      return Tensor();
    }
    // return a share view of block, this only work for non-symm tensor.
    const Tensor &get_block_(const std::vector<cytnx_int64> &qnum, const bool &force) const {
      cytnx_error_msg(
        true,
        "[ERROR][DenseUniTensor] try to get_block_() using qnum on a non-symmetry UniTensor%s",
        "\n");
      return this->_block;
    }
    Tensor &get_block_(const std::vector<cytnx_int64> &qnum, const bool &force) {
      cytnx_error_msg(
        true,
        "[ERROR][DenseUniTensor] try to get_block_() using qnum on a non-symmetry UniTensor%s",
        "\n");
      return this->_block;
    }
 
    // return a share view of block, this only work for non-symm tensor.
    Tensor &get_block_() { return this->_block; }
    Tensor &get_block_(const cytnx_uint64 &idx) {
      cytnx_error_msg(idx != 0,
                      "[ERROR][DenseUniTensor] Dense tensor has only one block, block number %llu "
                      "invalid. Use get_block_(0).\n",
                      (unsigned long long)idx);
      return this->_block;
    }
    // return a share view of block, this only work for non-symm tensor.
    const Tensor &get_block_() const { return this->_block; }
    const Tensor &get_block_(const cytnx_uint64 &idx) const {
      cytnx_error_msg(idx != 0,
                      "[ERROR][DenseUniTensor] Dense tensor has only one block, block number %llu "
                      "invalid. Use get_block_(0).\n",
                      (unsigned long long)idx);
      return this->_block;
    }
 
    cytnx_uint64 Nblocks() const { return 1; };
    std::vector<Tensor> get_blocks() const {
      std::vector<Tensor> out;
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] Cannot use get_blocks(), use get_block() instead!%s", "\n");
      return out;  // this will not share memory!!
    }
    const std::vector<Tensor> &get_blocks_(const bool &silent = false) const {
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] Cannot use get_blocks_(), use get_block_() instead!%s",
        "\n");
      return this->_interface_block;  // this will not share memory!!
    }
    std::vector<Tensor> &get_blocks_(const bool &silent = false) {
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] Cannot use get_blocks_(), use get_block_() instead!%s",
        "\n");
      return this->_interface_block;  // this will not share memory!!
    }
 
    void put_block(const Tensor &in) {
      // We don't check the dtype for DenseUniTensor, since it'll be more convenient to change
      // DenseUniTensor's dtype
 
      // cytnx_error_msg(in.dtype() != this->dtype(),
      //                 "[ERROR][DenseUniTensor][put_block] The input tensor dtype does not
      //                 match.%s",
      //                 "\n");
      cytnx_error_msg(in.device() != this->device(),
                      "[ERROR][DenseUniTensor][put_block] The input tensor device does not "
                      "match.%s",
                      "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      if (this->is_diag()) {
        cytnx_error_msg(
          in.shape() != this->_block.shape(),
          "[ERROR][DenseUniTensor][put_block] the input tensor shape does not match.%s", "\n");
        this->_block = in.clone();
      } else {
        cytnx_error_msg(
          in.shape() != this->shape(),
          "[ERROR][DenseUniTensor][put_block] the input tensor shape does not match.%s", "\n");
        this->_block = in.clone();
      }
    }
    void put_block(const Tensor &in, const cytnx_uint64 &idx) {
      cytnx_error_msg(idx != 0,
                      "[ERROR][DenseUniTensor] Dense tensor has only one block, block number %llu "
                      "invalid. Use put_block(0).\n",
                      (unsigned long long)idx);
      put_block(in);
    }
    // share view of the block
    void put_block_(Tensor &in) {
      // We don't check the dtype for DenseUniTensor, since it'll be more convenient to change
      // DenseUniTensor's dtype
 
      // cytnx_error_msg(in.dtype() != this->dtype(),
      //                 "[ERROR][DenseUniTensor][put_block] The input tensor dtype does not
      //                 match.%s",
      //                 "\n");
      cytnx_error_msg(in.device() != this->device(),
                      "[ERROR][DenseUniTensor][put_block] The input tensor device does not "
                      "match.%s",
                      "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      if (this->is_diag()) {
        cytnx_error_msg(
          in.shape() != this->_block.shape(),
          "[ERROR][DenseUniTensor][put_block] the input tensor shape does not match.%s", "\n");
        this->_block = in;
      } else {
        cytnx_error_msg(
          in.shape() != this->shape(),
          "[ERROR][DenseUniTensor][put_block] the input tensor shape does not match.%s", "\n");
        this->_block = in;
      }
    }
    void put_block_(Tensor &in, const cytnx_uint64 &idx) {
      cytnx_error_msg(idx != 0,
                      "[ERROR][DenseUniTensor] Dense tensor has only one block, block number %llu "
                      "invalid. Use put_block_(0).\n",
                      (unsigned long long)idx);
      put_block_(in);
    }
 
    void put_block(const Tensor &in, const std::vector<cytnx_int64> &qnum) {
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] try to put_block using qnum on a non-symmetry UniTensor%s",
        "\n");
    }
    void put_block_(Tensor &in, const std::vector<cytnx_int64> &qnum) {
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] try to put_block using qnum on a non-symmetry UniTensor%s",
        "\n");
    }
    // these two methods only work on non-symm tensor (DenseUniTensor)
    boost::intrusive_ptr<UniTensor_base> get(const std::vector<Accessor> &accessors);
    void set(const std::vector<Accessor> &accessors, const Tensor &rhs);
 
    void reshape_(const std::vector<cytnx_int64> &new_shape, const cytnx_uint64 &rowrank = 0);
    boost::intrusive_ptr<UniTensor_base> reshape(const std::vector<cytnx_int64> &new_shape,
                                                 const cytnx_uint64 &rowrank = 0);
    boost::intrusive_ptr<UniTensor_base> to_dense();
    void to_dense_();
    void combineBond(const std::vector<std::string> &indicators, const bool &force = true);
    void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force,
                      const bool &by_label);
    void combineBonds(const std::vector<std::string> &indicators, const bool &force = true);
    void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force = true);
    boost::intrusive_ptr<UniTensor_base> contract(const boost::intrusive_ptr<UniTensor_base> &rhs,
                                                  const bool &mv_elem_self = false,
                                                  const bool &mv_elem_rhs = false);
    std::vector<Bond> getTotalQnums(const bool &physical = false) {
      cytnx_error_msg(true, "[ERROR][DenseUniTensor] %s",
                      "getTotalQnums can only operate on UniTensor with symmetry.\n");
      return std::vector<Bond>();
    }
 
    std::vector<std::vector<cytnx_int64>> get_blocks_qnums() const {
      cytnx_error_msg(true, "[ERROR][DenseUniTensor] %s",
                      "get_blocks_qnums can only operate on UniTensor with symmetry.\n");
      return std::vector<std::vector<cytnx_int64>>();
    }
 
    bool same_data(const boost::intrusive_ptr<UniTensor_base> &rhs) const {
      if (rhs->uten_type() != UTenType.Dense) return false;
 
      return this->get_block_().same_data(rhs->get_block_());
    }
 
    ~DenseUniTensor(){};
 
    // arithmetic
    void Add_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Add_(const Scalar &rhs);
 
    void Mul_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Mul_(const Scalar &rhs);
 
    void Sub_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Sub_(const Scalar &rhs);
    void lSub_(const Scalar &lhs);
 
    void Div_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Div_(const Scalar &rhs);
    void lDiv_(const Scalar &lhs);
 
    void Conj_() { this->_block.Conj_(); };
 
    boost::intrusive_ptr<UniTensor_base> Conj() {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Conj_();
      return out;
    }
 
    boost::intrusive_ptr<UniTensor_base> Transpose() {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Transpose_();
      return out;
    }
    void Transpose_();
 
    boost::intrusive_ptr<UniTensor_base> normalize() {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->normalize_();
      return out;
    }
    void normalize_();
 
    boost::intrusive_ptr<UniTensor_base> Dagger() {
      boost::intrusive_ptr<UniTensor_base> out = this->Conj();
      out->Transpose_();
      return out;
    }
    void Dagger_() {
      this->Conj_();
      this->Transpose_();
    }
    void Trace_(const cytnx_int64 &a, const cytnx_int64 &b);
    void Trace_(const std::string &a, const std::string &b);
    boost::intrusive_ptr<UniTensor_base> Trace(const std::string &a, const std::string &b) {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Trace_(a, b);
      return out;
    }
    boost::intrusive_ptr<UniTensor_base> Trace(const cytnx_int64 &a, const cytnx_int64 &b) {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Trace_(a, b);
      return out;
    }
 
    Tensor Norm() const;
 
    const Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return Scalar::Sproxy();
    }
    const cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                          const cytnx_complex128 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return cytnx_complex128(0, 0);
    }
    const cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                         const cytnx_complex64 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return cytnx_complex64(0, 0);
    }
    const cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_double &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_float &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint64 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int64 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint32 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int32 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint16 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
    const cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int16 &aux) const {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return 0;
    }
 
    Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return Scalar::Sproxy();
    }
    cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                    const cytnx_complex128 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_complex128 *)nullptr;
    }
    cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                   const cytnx_complex64 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_complex64 *)nullptr;
    }
    cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_double &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_double *)nullptr;
    }
    cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_float &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_float *)nullptr;
    }
    cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint64 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_uint64 *)nullptr;
    }
    cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int64 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_int64 *)nullptr;
    }
    cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint32 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_uint32 *)nullptr;
    }
    cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int32 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_int32 *)nullptr;
    }
    cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint16 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_uint16 *)nullptr;
    }
    cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int16 &aux) {
      cytnx_error_msg(
        true, "[ERROR][Internal] This shouldn't be called by DenseUniTensor, something wrong.%s",
        "\n");
      return *(cytnx_int16 *)nullptr;
    }
 
    bool elem_exists(const std::vector<cytnx_uint64> &locator) const {
      cytnx_error_msg(
        true, "[ERROR][DenseUniTensor] elem_exists can only be used on UniTensor with Symmetry.%s",
        "\n");
    }
    void tag_() {
      if (!this->is_tag()) {
        for (int i = 0; i < this->_rowrank; i++) {
          this->_bonds[i] = this->_bonds[i].retype(BD_KET);
        }
        for (int i = this->_rowrank; i < this->_bonds.size(); i++) {
          this->_bonds[i] = this->_bonds[i].retype(BD_BRA);
        }
        this->_is_tag = true;
        this->_is_braket_form = this->_update_braket();
      }
    }
    void truncate_(const cytnx_int64 &bond_idx, const cytnx_uint64 &dim);
    void truncate_(const std::string &label, const cytnx_uint64 &dim);
 
    void from_(const boost::intrusive_ptr<UniTensor_base> &rhs, bool force, cytnx_double tol = 0.);
 
    void group_basis_() {
      cytnx_warning_msg(true, "[WARNING] group basis will not have any effect on DensUniTensor.%s",
                        "\n");
    }
 
    void _save_dispatch(std::fstream &f) const;
    void _load_dispatch(std::fstream &f, unsigned int version);
 
    const std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) const {
      cytnx_error_msg(true, "[ERROR] get_qindices can only be unsed on UniTensor with Symmetry.%s",
                      "\n");
    }
    std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) {
      cytnx_error_msg(true, "[ERROR] get_qindices can only be unsed on UniTensor with Symmetry.%s",
                      "\n");
    }
 
    const vec2d<cytnx_uint64> &get_itoi() const {
      cytnx_error_msg(true, "[ERROR] get_itoi can only be unsed on UniTensor with Symmetry.%s",
                      "\n");
    }
    vec2d<cytnx_uint64> &get_itoi() {
      cytnx_error_msg(true, "[ERROR] get_itoi can only be unsed on UniTensor with Symmetry.%s",
                      "\n");
    }
 
    // end virtual function
  };
 
  //======================================================================
  class BlockUniTensor : public UniTensor_base {
   protected:
   public:
    std::vector<std::vector<cytnx_uint64>>
      _inner_to_outer_idx;  // stores the qindices for each block
    std::vector<Tensor> _blocks;
    Tensor NullRefTensor;  // this returns when access block does not exists!
 
    // given an index list [loc], get qnums from this->_bonds[loc] and return the combined qnums
    // calculated from Symm object! this assume 1. symmetry are the same for each bond!
    //             2. total_qns are feeded with size len(symmetry)
    void _fx_get_total_fluxs(std::vector<cytnx_uint64> &loc, const std::vector<Symmetry> &syms,
                             std::vector<cytnx_int64> &total_qns) {
      cytnx_error_msg(
        loc.size() != this->_bonds.size(),
        "[ERROR][BlockUniTensor] qnum location rank does not match the number of bonds.%s", "\n");
      cytnx_error_msg(total_qns.size() != syms.size(),
                      "[ERROR][BlockUniTensor] total_qns size does not match symmetry count.%s",
                      "\n");
      std::fill(total_qns.begin(), total_qns.end(), 0);
      if (this->_bonds.empty()) return;
 
      for (cytnx_int32 i = 0; i < syms.size(); i++) {
        if (this->_bonds[0].type() == BD_BRA)
          total_qns[i] = syms[i].reverse_rule(this->_bonds[0]._impl->_qnums[loc[0]][i]);
        else
          total_qns[i] = this->_bonds[0]._impl->_qnums[loc[0]][i];
 
        for (auto j = 1; j < loc.size(); j++) {
          if (this->_bonds[j].type() == BD_BRA)
            total_qns[i] = syms[i].combine_rule(
              total_qns[i], syms[i].reverse_rule(this->_bonds[j]._impl->_qnums[loc[j]][i]));
          else {
            total_qns[i] =
              syms[i].combine_rule(total_qns[i], this->_bonds[j]._impl->_qnums[loc[j]][i]);
          }
        }
      }
    }
 
    void _fx_locate_elem(cytnx_int64 &bidx, std::vector<cytnx_uint64> &loc_in_T,
                         const std::vector<cytnx_uint64> &locator) const;
 
    // internal function, grouping all duplicate qnums in all bonds
    void _fx_group_duplicates(const std::vector<cytnx_uint64> &dup_bond_idxs,
                              const std::vector<std::vector<cytnx_uint64>> &idx_mappers);
 
    void set_meta(BlockUniTensor *tmp, const bool &inner, const bool &outer) const {
      // outer meta
      if (outer) {
        tmp->_bonds = vec_clone(this->_bonds);
        tmp->_labels = this->_labels;
        tmp->syms_cache = vec_clone(this->syms_cache);
        tmp->_is_braket_form = this->_is_braket_form;
        tmp->_rowrank = this->_rowrank;
        tmp->_name = this->_name;
      }
 
      tmp->_is_diag = this->_is_diag;
 
      // inner meta
      if (inner) {
        tmp->_inner_to_outer_idx = this->_inner_to_outer_idx;
      }
    }
 
    boost::intrusive_ptr<BlockUniTensor> clone_meta(const bool &inner, const bool &outer) const {
      boost::intrusive_ptr<BlockUniTensor> tmp(new BlockUniTensor());
      this->set_meta(tmp.get(), inner, outer);
      return tmp;
    };
 
    friend class UniTensor;
    BlockUniTensor() {
      this->uten_type_id = UTenType.Block;
      this->_is_tag = true;
    }
 
    // virtual functions:
    // void Init(const std::vector<Bond> &bonds, const std::vector<cytnx_int64> &in_labels = {},
    //           const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
    //           const int &device = Device.cpu, const bool &is_diag = false,
    //           const bool &no_alloc = false);
 
    void Init(const std::vector<Bond> &bonds, const std::vector<std::string> &in_labels = {},
              const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
              const int &device = Device.cpu, const bool &is_diag = false,
              const bool &no_alloc = false, const std::string &name = "");
 
    void Init_by_Tensor(const Tensor &in_tensor, const bool &is_diag = false,
                        const cytnx_int64 &rowrank = -1, const std::string &name = "") {
      cytnx_error_msg(
        true, "[ERROR][BlockUniTensor] Cannot use Init_by_tensor() on a BlockUniTensor.%s", "\n");
    }
 
    std::vector<cytnx_uint64> shape() const {
      if (this->rank() == 0) return {};
      std::vector<cytnx_uint64> out(this->_bonds.size());
      for (cytnx_uint64 i = 0; i < out.size(); i++) {
        out[i] = this->_bonds[i].dim();
      }
      return out;
    }
 
    bool is_blockform() const { return true; }
    bool is_contiguous() const {
      bool out = true;
      for (int i = 0; i < this->_blocks.size(); i++) {
        out &= this->_blocks[i].is_contiguous();
      }
      return out;
    };
 
    cytnx_uint64 Nblocks() const { return this->_blocks.size(); };
 
    void to_(const int &device) {
      for (cytnx_uint64 i = 0; i < this->_blocks.size(); i++) {
        this->_blocks[i].to_(device);
      }
    };
 
    boost::intrusive_ptr<UniTensor_base> to(const int &device) {
      if (this->device() == device) {
        return this;
      } else {
        boost::intrusive_ptr<UniTensor_base> out = this->clone();
        out->to_(device);
        return out;
      }
    };
 
    boost::intrusive_ptr<UniTensor_base> clone() const {
      boost::intrusive_ptr<BlockUniTensor> tmp = this->clone_meta(true, true);
      tmp->_blocks = vec_clone(this->_blocks);
      return tmp;
    };
 
    unsigned int dtype() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? Type.Void : this->_blocks[0].dtype();
    };
    int device() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? -404 : this->_blocks[0].device();
    };
    std::string dtype_str() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? "Void, no valid blocks" : this->_blocks[0].dtype_str();
    };
    std::string device_str() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? "None, no valid blocks" : this->_blocks[0].device_str();
    };
 
    Tensor get_block(const cytnx_uint64 &idx = 0) const {
      cytnx_error_msg(idx >= this->_blocks.size(), "[ERROR][BlockUniTensor] index out of range%s",
                      "\n");
      return this->_blocks[idx].clone();
    };
 
    // this one for Block will return the indicies!!
    Tensor get_block(const std::vector<cytnx_int64> &indices, const bool &force_return) const {
      if (this->rank() == 0) {
        cytnx_error_msg(!indices.empty(),
                        "[ERROR][get_block][BlockUniTensor] rank-0 scalar block expects no "
                        "qnum indices.%s",
                        "\n");
        cytnx_error_msg(this->_blocks.empty(), "[ERROR][BlockUniTensor] index out of range%s",
                        "\n");
        return this->_blocks[0].clone();
      }
      cytnx_error_msg(indices.size() != this->rank(),
                      "[ERROR][get_block][BlockUniTensor] len(indices) must be the same as the "
                      "Tensor rank (number of legs).%s",
                      "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        if (force_return) {
          return NullRefTensor;
        } else {
          cytnx_error_msg(true,
                          "[ERROR][get_block][BlockUniTensor] no avaliable block exists, "
                          "force_return=false, so "
                          "error throws. \n    If you want to return an empty block without "
                          "error when block is "
                          "not avaliable, set force_return=True.%s",
                          "\n");
        }
      } else {
        return this->_blocks[b].clone();
      }
    }
 
    const Tensor &get_block_(const cytnx_uint64 &idx = 0) const {
      cytnx_error_msg(idx >= this->_blocks.size(), "[ERROR][BlockUniTensor] index out of range%s",
                      "\n");
      return this->_blocks[idx];
    };
 
    Tensor &get_block_(const cytnx_uint64 &idx = 0) {
      cytnx_error_msg(idx >= this->_blocks.size(), "[ERROR][BlockUniTensor] index out of range%s",
                      "\n");
      return this->_blocks[idx];
    };
 
    const Tensor &get_block_(const std::vector<cytnx_int64> &indices,
                             const bool &force_return) const {
      if (this->rank() == 0) {
        cytnx_error_msg(!indices.empty(),
                        "[ERROR][get_block][BlockUniTensor] rank-0 scalar block expects no "
                        "qnum indices.%s",
                        "\n");
        cytnx_error_msg(this->_blocks.empty(), "[ERROR][BlockUniTensor] index out of range%s",
                        "\n");
        return this->_blocks[0];
      }
      cytnx_error_msg(indices.size() != this->rank(),
                      "[ERROR][get_block][BlockUniTensor] len(indices) must be the same as the "
                      "Tensor rank (number of legs).%s",
                      "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        if (force_return) {
          return this->NullRefTensor;
        } else {
          cytnx_error_msg(true,
                          "[ERROR][get_block][BlockUniTensor] no avaliable block exists, "
                          "force_return=false, so "
                          "error throws. \n    If you want to return an empty block without "
                          "error when block is "
                          "not avaliable, set force_return=True.%s",
                          "\n");
        }
      } else {
        return this->_blocks[b];
      }
    }
 
    Tensor &get_block_(const std::vector<cytnx_int64> &indices, const bool &force_return) {
      if (this->rank() == 0) {
        cytnx_error_msg(!indices.empty(),
                        "[ERROR][get_block][BlockUniTensor] rank-0 scalar block expects no "
                        "qnum indices.%s",
                        "\n");
        cytnx_error_msg(this->_blocks.empty(), "[ERROR][BlockUniTensor] index out of range%s",
                        "\n");
        return this->_blocks[0];
      }
      cytnx_error_msg(indices.size() != this->rank(),
                      "[ERROR][get_block][BlockUniTensor] len(indices) must be the same as the "
                      "Tensor rank (number of legs).%s",
                      "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        if (force_return) {
          return this->NullRefTensor;
        } else {
          cytnx_error_msg(true,
                          "[ERROR][get_block][BlockUniTensor] no avaliable block exists, "
                          "force_return=false, so "
                          "error throws. \n    If you want to return an empty block without "
                          "error when block is "
                          "not avaliable, set force_return=True.%s",
                          "\n");
        }
      } else {
        return this->_blocks[b];
      }
    }
 
    std::vector<Tensor> get_blocks() const { return vec_clone(this->_blocks); }
    const std::vector<Tensor> &get_blocks_(const bool &) const { return this->_blocks; }
    std::vector<Tensor> &get_blocks_(const bool &) { return this->_blocks; }
 
    bool same_data(const boost::intrusive_ptr<UniTensor_base> &rhs) const {
      if (rhs->uten_type() != UTenType.Block) return false;
      if (rhs->get_blocks_(1).size() != this->get_blocks_(1).size()) return false;
 
      for (int i = 0; i < rhs->get_blocks_(1).size(); i++)
        if (this->get_blocks_(1)[i].same_data(rhs->get_blocks_(1)[i]) == false) return false;
 
      return true;
    }
 
    void set_rowrank_(const cytnx_uint64 &new_rowrank) {
      cytnx_error_msg(new_rowrank > this->rank(),
                      "[ERROR][BlockUniTensor] rowrank should be [>=0] and [<=UniTensor.rank].%s",
                      "\n");
      if (this->is_diag()) {
        cytnx_error_msg(new_rowrank != 1,
                        "[ERROR][BlockUniTensor] rowrank should be [==1] when is_diag =true!.%s",
                        "\n");
      }
      this->_rowrank = new_rowrank;
      this->_is_braket_form = this->_update_braket();
    }
 
    boost::intrusive_ptr<UniTensor_base> set_rowrank(const cytnx_uint64 &new_rowrank) const {
      boost::intrusive_ptr<BlockUniTensor> tmp = this->clone_meta(true, true);
      tmp->_blocks = this->_blocks;
      tmp->set_rowrank_(new_rowrank);
      return tmp;
    }
 
    boost::intrusive_ptr<UniTensor_base> permute(const std::vector<cytnx_int64> &mapper,
                                                 const cytnx_int64 &rowrank = -1);
    boost::intrusive_ptr<UniTensor_base> permute(const std::vector<std::string> &mapper,
                                                 const cytnx_int64 &rowrank = -1);
 
    void permute_(const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1);
    void permute_(const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1);
 
    void twist_(const cytnx_int64 &idx) override {
      // do nothing for bosonic UniTensor
      return;
    }
    void fermion_twists_() override {
      // do nothing for bosonic UniTensor
      return;
    }
    void twist_(const std::string &label) override {
      // do nothing for bosonic UniTensor
      return;
    }
 
    boost::intrusive_ptr<UniTensor_base> contiguous_() {
      // Rebind each block to the non-mutating Tensor::contiguous() result instead of calling
      // contiguous_() in place, so a block Tensor shared with another UniTensor (e.g. via
      // relabel()) is not corrupted (#724). Tensor::contiguous() is cheap when the block is
      // already contiguous (no data copy; same impl returned).
      for (unsigned int b = 0; b < this->_blocks.size(); b++)
        this->_blocks[b] = this->_blocks[b].contiguous();
      return boost::intrusive_ptr<UniTensor_base>(this);
    }
    boost::intrusive_ptr<UniTensor_base> contiguous();
 
    boost::intrusive_ptr<UniTensor_base> apply_() {
      return boost::intrusive_ptr<UniTensor_base>(this);
    }
    boost::intrusive_ptr<UniTensor_base> apply() {
      // just return self
      boost::intrusive_ptr<UniTensor_base> out(this);
      return out;
    }
 
    void print_diagram(const bool &bond_info = false) const;
    void print_blocks(const bool &full_info = true) const;
    void print_block(const cytnx_int64 &idx, const bool &full_info = true) const;
 
    boost::intrusive_ptr<UniTensor_base> contract(const boost::intrusive_ptr<UniTensor_base> &rhs,
                                                  const bool &mv_elem_self = false,
                                                  const bool &mv_elem_rhs = false);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::vector<std::string> &new_labels);
    boost::intrusive_ptr<UniTensor_base> relabels(const std::vector<std::string> &new_labels);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::vector<std::string> &old_labels,
                                                 const std::vector<std::string> &new_labels);
    boost::intrusive_ptr<UniTensor_base> relabels(const std::vector<std::string> &old_labels,
                                                  const std::vector<std::string> &new_labels);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::string &old_label,
                                                 const std::string &new_label);
    boost::intrusive_ptr<UniTensor_base> relabel(const cytnx_int64 &inx,
                                                 const std::string &new_label);
 
    std::vector<Symmetry> syms() const;
 
    void reshape_(const std::vector<cytnx_int64> &new_shape, const cytnx_uint64 &rowrank = 0) {
      cytnx_error_msg(true, "[ERROR] Cannot reshape a UniTensor with symmetry.%s", "\n");
    }
    boost::intrusive_ptr<UniTensor_base> reshape(const std::vector<cytnx_int64> &new_shape,
                                                 const cytnx_uint64 &rowrank = 0) {
      cytnx_error_msg(true, "[ERROR] Cannot reshape a UniTensor with symmetry.%s", "\n");
      return nullptr;
    }
 
    boost::intrusive_ptr<UniTensor_base> to_dense();
    void to_dense_();
 
    boost::intrusive_ptr<UniTensor_base> astype(const unsigned int &dtype) const {
      boost::intrusive_ptr<BlockUniTensor> tmp = this->clone_meta(true, true);
      tmp->_blocks.resize(this->_blocks.size());
      for (cytnx_int64 blk = 0; blk < this->_blocks.size(); blk++) {
        tmp->_blocks[blk] = this->_blocks[blk].astype(dtype);
      }
      return tmp;
    };
 
    // this will only work on non-symm tensor (DenseUniTensor)
    boost::intrusive_ptr<UniTensor_base> get(const std::vector<Accessor> &accessors) {
      cytnx_error_msg(
        true,
        "[ERROR][BlockUniTensor][get] Cannot use get on a UniTensor with "
        "Symmetry.\n suggestion: try get_block/get_block_/get_blocks/get_blocks_ first.%s",
        "\n");
      return nullptr;
    }
 
    // this will only work on non-symm tensor (DenseUniTensor)
    void set(const std::vector<Accessor> &accessors, const Tensor &rhs) {
      cytnx_error_msg(
        true,
        "[ERROR][BlockUniTensor][get] Cannot use get on a UniTensor with "
        "Symmetry.\n suggestion: try get_block/get_block_/get_blocks/get_blocks_ first.%s",
        "\n");
    }
 
    void put_block(const Tensor &in, const cytnx_uint64 &idx = 0) {
      cytnx_error_msg(in.dtype() != this->dtype(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor dtype does not match.%s",
                      "\n");
      cytnx_error_msg(in.device() != this->device(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor device does not "
                      "match.%s",
                      "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(idx >= this->_blocks.size(), "[ERROR][BlockUniTensor] index out of range%s",
                      "\n");
      cytnx_error_msg(in.shape() != this->_blocks[idx].shape(),
                      "[ERROR][BlockUniTensor] the shape of input tensor does not match the shape "
                      "of block @ idx=%d\n",
                      idx);
 
      this->_blocks[idx] = in.clone();
    }
    void put_block_(Tensor &in, const cytnx_uint64 &idx = 0) {
      cytnx_error_msg(in.dtype() != this->dtype(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor dtype does not match.%s",
                      "\n");
      cytnx_error_msg(in.device() != this->device(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor device does not "
                      "match.%s",
                      "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(idx >= this->_blocks.size(), "[ERROR][BlockUniTensor] index out of range%s",
                      "\n");
      cytnx_error_msg(in.shape() != this->_blocks[idx].shape(),
                      "[ERROR][BlockUniTensor] the shape of input tensor does not match the shape "
                      "of block @ idx=%d\n",
                      idx);
 
      this->_blocks[idx] = in;
    }
    void put_block(const Tensor &in, const std::vector<cytnx_int64> &indices) {
      cytnx_error_msg(in.dtype() != this->dtype(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor dtype does not match.%s",
                      "\n");
      cytnx_error_msg(in.device() != this->device(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor device does not "
                      "match.%s",
                      "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(indices.size() != this->rank(),
                      "[ERROR][put_block][BlockUniTensor] len(indices) must be the same as the "
                      "Tensor rank (number of legs).%s",
                      "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        cytnx_error_msg(true, "[ERROR][put_block][BlockUniTensor] no avaliable block exists.%s",
                        "\n");
      } else {
        cytnx_error_msg(
          in.shape() != this->_blocks[b].shape(),
          "[ERROR][BlockUniTensor] the shape of input tensor does not match the shape "
          "of block @ idx=%d\n",
          b);
 
        this->_blocks[b] = in.clone();
      }
    }
    void put_block_(Tensor &in, const std::vector<cytnx_int64> &indices) {
      cytnx_error_msg(in.dtype() != this->dtype(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor dtype does not match.%s",
                      "\n");
      cytnx_error_msg(in.device() != this->device(),
                      "[ERROR][BlockUniTensor][put_block] The input tensor device does not "
                      "match.%s",
                      "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(indices.size() != this->rank(),
                      "[ERROR][put_block][BlockUniTensor] len(indices) must be the same as the "
                      "Tensor rank (number of legs).%s",
                      "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        cytnx_error_msg(true, "[ERROR][put_block][BlockUniTensor] no avaliable block exists.%s",
                        "\n");
      } else {
        cytnx_error_msg(
          in.shape() != this->_blocks[b].shape(),
          "[ERROR][BlockUniTensor] the shape of input tensor does not match the shape "
          "of block @ idx=%d\n",
          b);
        this->_blocks[b] = in;
      }
    }
 
    void tag_() {
      // no-use!
    }
 
    boost::intrusive_ptr<UniTensor_base> Conj() {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Conj_();
      return out;
    }
 
    void Conj_() {
      for (int i = 0; i < this->_blocks.size(); i++) {
        this->_blocks[i].Conj_();
      }
    };
 
    void Transpose_();
    boost::intrusive_ptr<UniTensor_base> Transpose() {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Transpose_();
      return out;
    }
 
    void normalize_();
    boost::intrusive_ptr<UniTensor_base> normalize() {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->normalize_();
      return out;
    }
 
    boost::intrusive_ptr<UniTensor_base> Dagger() {
      boost::intrusive_ptr<UniTensor_base> out = this->Conj();
      out->Transpose_();
      return out;
    }
    void Dagger_() {
      this->Conj_();
      this->Transpose_();
    }
 
    void Trace_(const std::string &a, const std::string &b);
    void Trace_(const cytnx_int64 &a, const cytnx_int64 &b);
 
    boost::intrusive_ptr<UniTensor_base> Trace(const std::string &a, const std::string &b) {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Trace_(a, b);
      return out;
    }
    boost::intrusive_ptr<UniTensor_base> Trace(const cytnx_int64 &a, const cytnx_int64 &b) {
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Trace_(a, b);
      return out;
    }
 
    Tensor Norm() const;
 
    bool elem_exists(const std::vector<cytnx_uint64> &locator) const;
 
    const Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator) const;
    const cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                          const cytnx_complex128 &aux) const;
    const cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                         const cytnx_complex64 &aux) const;
    const cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_double &aux) const;
    const cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_float &aux) const;
    const cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint64 &aux) const;
    const cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int64 &aux) const;
    const cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint32 &aux) const;
    const cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int32 &aux) const;
    const cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint16 &aux) const;
    const cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int16 &aux) const;
 
    Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator);
    cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                    const cytnx_complex128 &aux);
    cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                   const cytnx_complex64 &aux);
    cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_double &aux);
    cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_float &aux);
    cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint64 &aux);
    cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int64 &aux);
    cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint32 &aux);
    cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int32 &aux);
    cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint16 &aux);
    cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int16 &aux);
 
    void _save_dispatch(std::fstream &f) const;
    void _load_dispatch(std::fstream &f, unsigned int version);
 
    // this will remove the [q_index]-th qnum at [bond_idx]-th Bond!
    void truncate_(const std::string &label, const cytnx_uint64 &q_index);
    void truncate_(const cytnx_int64 &bond_idx, const cytnx_uint64 &q_index);
 
    void Add_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Add_(const Scalar &rhs) {
      cytnx_error_msg(
        true,
        "[ERROR] Cannot perform elementwise arithmetic '+' between Scalar and BlockUniTensor.\n %s "
        "\n",
        "This operation would destroy the block structure. [Suggest] Avoid or use get/put_block(s) "
        "to do operation on blocks.");
    }
 
    void Mul_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Mul_(const Scalar &rhs);
 
    void Sub_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Sub_(const Scalar &rhs) {
      cytnx_error_msg(
        true,
        "[ERROR] Cannot perform elementwise arithmetic '-' between Scalar and BlockUniTensor.\n %s "
        "\n",
        "This operation would destroy the block structure. [Suggest] Avoid or use get/put_block(s) "
        "to do operation on blocks.");
    }
    void lSub_(const Scalar &lhs) {
      cytnx_error_msg(
        true,
        "[ERROR] Cannot perform elementwise arithmetic '-' between Scalar and BlockUniTensor.\n %s "
        "\n",
        "This operation would destroy the block structure. [Suggest] Avoid or use get/put_block(s) "
        "to do operation on blocks.");
    }
 
    void Div_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Div_(const Scalar &rhs);
    void lDiv_(const Scalar &lhs) {
      cytnx_error_msg(
        true,
        "[ERROR] Cannot perform elementwise arithmetic '/' between Scalar and BlockUniTensor.\n %s "
        "\n",
        "This operation would cause division by zero on non-block elements. [Suggest] Avoid or use "
        "get/put_block(s) to do operation on blocks.");
    }
    void from_(const boost::intrusive_ptr<UniTensor_base> &rhs, bool force, cytnx_double tol = 0.);
 
    void group_basis_();
 
    void combineBond(const std::vector<std::string> &indicators, const bool &force = false);
    void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force = false);
    void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force,
                      const bool &by_label);
    void combineBonds(const std::vector<std::string> &indicators, const bool &force = false);
 
    const std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) const {
      cytnx_error_msg(
        bidx >= this->Nblocks(),
        "[ERROR][BlockUniTensor] bidx out of bound! only %d blocks in current UTen.\n",
        this->Nblocks());
      return this->_inner_to_outer_idx[bidx];
    }
    std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) {
      cytnx_error_msg(
        bidx >= this->Nblocks(),
        "[ERROR][BlockUniTensor] bidx out of bound! only %d blocks in current UTen.\n",
        this->Nblocks());
      return this->_inner_to_outer_idx[bidx];
    }
 
    const vec2d<cytnx_uint64> &get_itoi() const { return this->_inner_to_outer_idx; }
    vec2d<cytnx_uint64> &get_itoi() { return this->_inner_to_outer_idx; }
  };
 
  //======================================================================
  class BlockFermionicUniTensor : public UniTensor_base {
    //[21 Aug 2024] This is a copy from BlockUniTensor; additionally sign flips are stored as
    //_signflip and taken care of in all the methods
   public:
    std::vector<std::vector<cytnx_uint64>>
      _inner_to_outer_idx;  // stores the qindices for each block
    std::vector<Tensor> _blocks;
    Tensor NullRefTensor;  // this returns when accessed block does not exists!
 
    // given an index list [loc], get qnums from this->_bonds[loc] and return the combined qnums
    // calculated from Symm object! this assume 1. symmetry are the same for each bond!
    //             2. total_qns are feeded with size len(symmetry)
    void _fx_get_total_fluxs(std::vector<cytnx_uint64> &loc, const std::vector<Symmetry> &syms,
                             std::vector<cytnx_int64> &total_qns) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        loc.size() != this->_bonds.size(),
        "[ERROR][BlockFermionicUniTensor] qnum location rank does not match the number of bonds.%s",
        "\n");
      cytnx_error_msg(
        total_qns.size() != syms.size(),
        "[ERROR][BlockFermionicUniTensor] total_qns size does not match symmetry count.%s", "\n");
      std::fill(total_qns.begin(), total_qns.end(), 0);
      if (this->_bonds.empty()) return;
 
      for (cytnx_int32 i = 0; i < syms.size(); i++) {
        if (this->_bonds[0].type() == BD_BRA)
          total_qns[i] = syms[i].reverse_rule(this->_bonds[0]._impl->_qnums[loc[0]][i]);
        else
          total_qns[i] = this->_bonds[0]._impl->_qnums[loc[0]][i];
 
        for (auto j = 1; j < loc.size(); j++) {
          if (this->_bonds[j].type() == BD_BRA)
            total_qns[i] = syms[i].combine_rule(
              total_qns[i], syms[i].reverse_rule(this->_bonds[j]._impl->_qnums[loc[j]][i]));
          else {
            total_qns[i] =
              syms[i].combine_rule(total_qns[i], this->_bonds[j]._impl->_qnums[loc[j]][i]);
          }
        }
      }
    }
 
    void _fx_locate_elem(cytnx_int64 &bidx, std::vector<cytnx_uint64> &loc_in_T,
                         const std::vector<cytnx_uint64> &locator) const;
 
    // internal function, grouping all duplicate qnums in all bonds
    void _fx_group_duplicates(const std::vector<cytnx_uint64> &dup_bond_idxs,
                              const std::vector<std::vector<cytnx_uint64>> &idx_mappers);
 
    void set_meta(BlockFermionicUniTensor *tmp, const bool &inner, const bool &outer) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor; additionally, _signflip is set
      // outer meta
      if (outer) {
        tmp->_bonds = vec_clone(this->_bonds);
        tmp->_labels = this->_labels;
        tmp->syms_cache = vec_clone(this->syms_cache);
        tmp->_is_braket_form = this->_is_braket_form;
        tmp->_rowrank = this->_rowrank;
        tmp->_name = this->_name;
        // tmp->_signflip = vec_clone(this->_signflip);
        tmp->_signflip = this->_signflip;
      }
 
      tmp->_is_diag = this->_is_diag;
 
      // inner meta
      if (inner) {
        tmp->_inner_to_outer_idx = this->_inner_to_outer_idx;
      }
    }
 
    boost::intrusive_ptr<BlockFermionicUniTensor> clone_meta(const bool &inner,
                                                             const bool &outer) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      boost::intrusive_ptr<BlockFermionicUniTensor> tmp(new BlockFermionicUniTensor());
      this->set_meta(tmp.get(), inner, outer);
      return tmp;
    };
 
    friend class UniTensor;
    BlockFermionicUniTensor() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      this->uten_type_id = UTenType.BlockFermionic;
      this->_is_tag = true;
    }
 
    void Init(const std::vector<Bond> &bonds, const std::vector<std::string> &in_labels = {},
              const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
              const int &device = Device.cpu, const bool &is_diag = false,
              const bool &no_alloc = false, const std::string &name = "");
 
    void Init_by_Tensor(const Tensor &in_tensor, const bool &is_diag = false,
                        const cytnx_int64 &rowrank = -1, const std::string &name = "") {
      cytnx_error_msg(true,
                      "[ERROR][BlockFermionicUniTensor] Cannot use Init_by_tensor() on a "
                      "BlockFermionicUniTensor.%s",
                      "\n");
    }
 
    std::vector<cytnx_uint64> shape() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      if (this->rank() == 0) return {};
      std::vector<cytnx_uint64> out(this->_bonds.size());
      for (cytnx_uint64 i = 0; i < out.size(); i++) {
        out[i] = this->_bonds[i].dim();
      }
      return out;
    }
 
    bool is_blockform() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return true;
    }
 
    bool is_contiguous() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      bool out = true;
      for (int i = 0; i < this->_blocks.size(); i++) {
        out &= this->_blocks[i].is_contiguous();
      }
      return out;
    };
 
    cytnx_uint64 Nblocks() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return this->_blocks.size();
    };
 
    std::vector<bool> signflip() const override { return this->_signflip; };
 
    void reset_signflip_() { this->_signflip.assign(this->_blocks.size(), false); }
 
    void erase_signflip_(const std::vector<cytnx_uint64> &positions);
 
    void to_(const int &device) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      for (cytnx_uint64 i = 0; i < this->_blocks.size(); i++) {
        this->_blocks[i].to_(device);
      }
    };
 
    boost::intrusive_ptr<UniTensor_base> to(const int &device) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      if (this->device() == device) {
        return this;
      } else {
        boost::intrusive_ptr<UniTensor_base> out = this->clone();
        out->to_(device);
        return out;
      }
    };
 
    boost::intrusive_ptr<UniTensor_base> clone() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor; changed to BlockFermionicUniTensor as
      // output
      boost::intrusive_ptr<BlockFermionicUniTensor> tmp = this->clone_meta(true, true);
      tmp->_blocks = vec_clone(this->_blocks);
      return tmp;
    };
 
    unsigned int dtype() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? Type.Void : this->_blocks[0].dtype();
    };
    int device() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? -404 : this->_blocks[0].device();
    };
    std::string dtype_str() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? "Void, no valid blocks" : this->_blocks[0].dtype_str();
    };
    std::string device_str() const {
        //[21 Aug 2024] This is a copy from BlockUniTensor;
  #ifdef UNI_DEBUG
      cytnx_error_msg(this->_blocks.size() == 0, "[ERROR][internal] empty blocks for blockform.%s",
                      "\n");
  #endif
      return this->_blocks.size() < 1 ? "None, no valid blocks" : this->_blocks[0].device_str();
    };
 
    Tensor get_block(const cytnx_uint64 &idx = 0) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(idx >= this->_blocks.size(),
                      "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
      return this->_blocks[idx].clone();
    };
 
    // this one for Block will return the indicies!!
    Tensor get_block(const std::vector<cytnx_int64> &indices, const bool &force_return) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      if (this->rank() == 0) {
        cytnx_error_msg(!indices.empty(),
                        "[ERROR][get_block][BlockFermionicUniTensor] rank-0 scalar block expects "
                        "no qnum indices.%s",
                        "\n");
        cytnx_error_msg(this->_blocks.empty(),
                        "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
        return this->_blocks[0].clone();
      }
      cytnx_error_msg(
        indices.size() != this->rank(),
        "[ERROR][get_block][BlockFermionicUniTensor] len(indices) must be the same as the "
        "Tensor rank (number of legs).%s",
        "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        if (force_return) {
          return NullRefTensor;
        } else {
          cytnx_error_msg(true,
                          "[ERROR][get_block][BlockFermionicUniTensor] no avaliable block exists, "
                          "force_return=false, so "
                          "error throws. \n    If you want to return an empty block without "
                          "error when block is "
                          "not avaliable, set force_return=True.%s",
                          "\n");
        }
      } else {
        return this->_blocks[b].clone();
      }
    }
 
    const Tensor &get_block_(const cytnx_uint64 &idx = 0) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(idx >= this->_blocks.size(),
                      "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
      return this->_blocks[idx];
    };
 
    Tensor &get_block_(const cytnx_uint64 &idx = 0) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(idx >= this->_blocks.size(),
                      "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
      return this->_blocks[idx];
    };
 
    const Tensor &get_block_(const std::vector<cytnx_int64> &indices,
                             const bool &force_return) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      if (this->rank() == 0) {
        cytnx_error_msg(!indices.empty(),
                        "[ERROR][get_block][BlockFermionicUniTensor] rank-0 scalar block expects "
                        "no qnum indices.%s",
                        "\n");
        cytnx_error_msg(this->_blocks.empty(),
                        "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
        return this->_blocks[0];
      }
      cytnx_error_msg(
        indices.size() != this->rank(),
        "[ERROR][get_block][BlockFermionicUniTensor] len(indices) must be the same as the "
        "Tensor rank (number of legs).%s",
        "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        if (force_return) {
          return this->NullRefTensor;
        } else {
          cytnx_error_msg(true,
                          "[ERROR][get_block][BlockFermionicUniTensor] no avaliable block exists, "
                          "force_return=false, so "
                          "error throws. \n    If you want to return an empty block without "
                          "error when block is "
                          "not avaliable, set force_return=True.%s",
                          "\n");
        }
      } else {
        return this->_blocks[b];
      }
    }
 
    Tensor &get_block_(const std::vector<cytnx_int64> &indices, const bool &force_return) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      if (this->rank() == 0) {
        cytnx_error_msg(!indices.empty(),
                        "[ERROR][get_block][BlockFermionicUniTensor] rank-0 scalar block expects "
                        "no qnum indices.%s",
                        "\n");
        cytnx_error_msg(this->_blocks.empty(),
                        "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
        return this->_blocks[0];
      }
      cytnx_error_msg(
        indices.size() != this->rank(),
        "[ERROR][get_block][BlockFermionicUniTensor] len(indices) must be the same as the "
        "Tensor rank (number of legs).%s",
        "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        if (force_return) {
          return this->NullRefTensor;
        } else {
          cytnx_error_msg(true,
                          "[ERROR][get_block][BlockFermionicUniTensor] no avaliable block exists, "
                          "force_return=false, so "
                          "error throws. \n    If you want to return an empty block without "
                          "error when block is "
                          "not avaliable, set force_return=True.%s",
                          "\n");
        }
      } else {
        return this->_blocks[b];
      }
    }
 
    std::vector<Tensor> get_blocks() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return vec_clone(this->_blocks);
    }
    const std::vector<Tensor> &get_blocks_(const bool &) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return this->_blocks;
    }
    std::vector<Tensor> &get_blocks_(const bool &) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return this->_blocks;
    }
 
    bool same_data(const boost::intrusive_ptr<UniTensor_base> &rhs) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor; changed to BlockFermionic UTenType
      if (rhs->uten_type() != UTenType.BlockFermionic) return false;
      if (rhs->get_blocks_(1).size() != this->get_blocks_(1).size()) return false;
 
      for (int i = 0; i < rhs->get_blocks_(1).size(); i++)
        if (this->get_blocks_(1)[i].same_data(rhs->get_blocks_(1)[i]) == false) return false;
 
      return true;
    }
 
    void set_rowrank_(const cytnx_uint64 &new_rowrank) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        new_rowrank > this->rank(),
        "[ERROR][BlockFermionicUniTensor] rowrank should be [>=0] and [<=UniTensor.rank].%s", "\n");
      if (this->is_diag()) {
        cytnx_error_msg(
          new_rowrank != 1,
          "[ERROR][BlockFermionicUniTensor] rowrank should be [==1] when is_diag =true!.%s", "\n");
      }
      this->_rowrank = new_rowrank;
      this->_is_braket_form = this->_update_braket();
    }
 
    boost::intrusive_ptr<UniTensor_base> set_rowrank(const cytnx_uint64 &new_rowrank) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor; output type changed to
      // BlockFermionicUniTensor
      boost::intrusive_ptr<BlockFermionicUniTensor> tmp = this->clone_meta(true, true);
      tmp->_blocks = this->_blocks;
      tmp->set_rowrank_(new_rowrank);
      return tmp;
    }
 
    boost::intrusive_ptr<UniTensor_base> permute(const std::vector<cytnx_int64> &mapper,
                                                 const cytnx_int64 &rowrank = -1);
    boost::intrusive_ptr<UniTensor_base> permute(const std::vector<std::string> &mapper,
                                                 const cytnx_int64 &rowrank = -1);
 
    void permute_(const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1) override;
    void permute_(const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1) override;
 
    boost::intrusive_ptr<UniTensor_base> permute_nosignflip(
      const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1) override;
    boost::intrusive_ptr<UniTensor_base> permute_nosignflip(
      const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1) override;
    void permute_nosignflip_(const std::vector<cytnx_int64> &mapper,
                             const cytnx_int64 &rowrank = -1) override;
    void permute_nosignflip_(const std::vector<std::string> &mapper,
                             const cytnx_int64 &rowrank = -1) override;
 
    void twist_(const cytnx_int64 &idx) override;
    void twist_(const std::string &label) override;
    void fermion_twists_() override;
 
    // Helper function; implements the sign flips when permuting indices
    std::vector<bool> _swapsigns_(const std::vector<cytnx_int64> &mapper) const;
    std::vector<bool> _lhssigns_(const std::vector<cytnx_int64> &mapper,
                                 const cytnx_int64 contrno) const;
    std::vector<bool> _swapsigns_(const std::vector<cytnx_uint64> &mapper) const;
    std::vector<bool> _lhssigns_(const std::vector<cytnx_uint64> &mapper,
                                 const cytnx_uint64 contrno) const;
 
    boost::intrusive_ptr<UniTensor_base> contiguous_() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      // Rebind each block to the non-mutating Tensor::contiguous() result instead of calling
      // contiguous_() in place, so a block Tensor shared with another UniTensor (e.g. via
      // relabel()) is not corrupted (#724).
      for (unsigned int b = 0; b < this->_blocks.size(); b++)
        this->_blocks[b] = this->_blocks[b].contiguous();
      return boost::intrusive_ptr<UniTensor_base>(this);
    }
    boost::intrusive_ptr<UniTensor_base> contiguous();
 
    boost::intrusive_ptr<UniTensor_base> apply_();
    boost::intrusive_ptr<UniTensor_base> apply();
 
    void print_diagram(const bool &bond_info = false) const;
    void print_blocks(const bool &full_info = true) const;
    void print_block(const cytnx_int64 &idx, const bool &full_info = true) const;
 
    boost::intrusive_ptr<UniTensor_base> contract(const boost::intrusive_ptr<UniTensor_base> &rhs,
                                                  const bool &mv_elem_self = false,
                                                  const bool &mv_elem_rhs = false);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::vector<std::string> &new_labels);
    boost::intrusive_ptr<UniTensor_base> relabels(const std::vector<std::string> &new_labels);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::vector<std::string> &old_labels,
                                                 const std::vector<std::string> &new_labels);
    boost::intrusive_ptr<UniTensor_base> relabels(const std::vector<std::string> &old_labels,
                                                  const std::vector<std::string> &new_labels);
 
    boost::intrusive_ptr<UniTensor_base> relabel(const std::string &old_label,
                                                 const std::string &new_label);
    boost::intrusive_ptr<UniTensor_base> relabel(const cytnx_int64 &inx,
                                                 const std::string &new_label);
    std::vector<Symmetry> syms() const;
 
    void reshape_(const std::vector<cytnx_int64> &new_shape, const cytnx_uint64 &rowrank = 0) {
      cytnx_error_msg(true, "[ERROR] Cannot reshape a UniTensor with symmetry.%s", "\n");
    }
    boost::intrusive_ptr<UniTensor_base> reshape(const std::vector<cytnx_int64> &new_shape,
                                                 const cytnx_uint64 &rowrank = 0) {
      cytnx_error_msg(true, "[ERROR] Cannot reshape a UniTensor with symmetry.%s", "\n");
      return nullptr;
    }
 
    boost::intrusive_ptr<UniTensor_base> to_dense();
    void to_dense_();
 
    boost::intrusive_ptr<UniTensor_base> astype(const unsigned int &dtype) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor; the tensor type was adapted
      boost::intrusive_ptr<BlockFermionicUniTensor> tmp = this->clone_meta(true, true);
      tmp->_blocks.resize(this->_blocks.size());
      for (cytnx_int64 blk = 0; blk < this->_blocks.size(); blk++) {
        tmp->_blocks[blk] = this->_blocks[blk].astype(dtype);
      }
      return tmp;
    };
 
    // this will only work on non-symm tensor (DenseUniTensor)
    boost::intrusive_ptr<UniTensor_base> get(const std::vector<Accessor> &accessors) {
      cytnx_error_msg(
        true,
        "[ERROR][BlockFermionicUniTensor][get] Cannot use get on a UniTensor with "
        "Symmetry.\n suggestion: try get_block/get_block_/get_blocks/get_blocks_ first.%s",
        "\n");
      return nullptr;
    }
 
    // this will only work on non-symm tensor (DenseUniTensor)
    void set(const std::vector<Accessor> &accessors, const Tensor &rhs) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        true,
        "[ERROR][BlockFermionicUniTensor][get] Cannot use get on a UniTensor with "
        "Symmetry.\n suggestion: try get_block/get_block_/get_blocks/get_blocks_ first.%s",
        "\n");
    }
 
    void put_block(const Tensor &in, const cytnx_uint64 &idx = 0) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        in.dtype() != this->dtype(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor dtype does not match.%s",
        "\n");
      cytnx_error_msg(
        in.device() != this->device(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor device does not "
        "match.%s",
        "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(idx >= this->_blocks.size(),
                      "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
      cytnx_error_msg(
        in.shape() != this->_blocks[idx].shape(),
        "[ERROR][BlockFermionicUniTensor] the shape of input tensor does not match the shape "
        "of block @ idx=%d\n",
        idx);
 
      this->_blocks[idx] = in.clone();
    }
    void put_block_(Tensor &in, const cytnx_uint64 &idx = 0) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        in.dtype() != this->dtype(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor dtype does not match.%s",
        "\n");
      cytnx_error_msg(
        in.device() != this->device(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor device does not "
        "match.%s",
        "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(idx >= this->_blocks.size(),
                      "[ERROR][BlockFermionicUniTensor] index out of range%s", "\n");
      cytnx_error_msg(
        in.shape() != this->_blocks[idx].shape(),
        "[ERROR][BlockFermionicUniTensor] the shape of input tensor does not match the shape "
        "of block @ idx=%d\n",
        idx);
 
      this->_blocks[idx] = in;
    }
    void put_block(const Tensor &in, const std::vector<cytnx_int64> &indices) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        in.dtype() != this->dtype(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor dtype does not match.%s",
        "\n");
      cytnx_error_msg(
        in.device() != this->device(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor device does not "
        "match.%s",
        "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(
        indices.size() != this->rank(),
        "[ERROR][put_block][BlockFermionicUniTensor] len(indices) must be the same as the "
        "Tensor rank (number of legs).%s",
        "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        cytnx_error_msg(
          true, "[ERROR][put_block][BlockFermionicUniTensor] no avaliable block exists.%s", "\n");
      } else {
        cytnx_error_msg(
          in.shape() != this->_blocks[b].shape(),
          "[ERROR][BlockFermionicUniTensor] the shape of input tensor does not match the shape "
          "of block @ idx=%d\n",
          b);
 
        this->_blocks[b] = in.clone();
      }
    }
    void put_block_(Tensor &in, const std::vector<cytnx_int64> &indices) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        in.dtype() != this->dtype(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor dtype does not match.%s",
        "\n");
      cytnx_error_msg(
        in.device() != this->device(),
        "[ERROR][BlockFermionicUniTensor][put_block] The input tensor device does not "
        "match.%s",
        "\n");
      // We shouldn't check the contiguous
      // cytnx_error_msg(!in.contiguous());
      cytnx_error_msg(
        indices.size() != this->rank(),
        "[ERROR][put_block][BlockFermionicUniTensor] len(indices) must be the same as the "
        "Tensor rank (number of legs).%s",
        "\n");
 
      std::vector<cytnx_uint64> inds(indices.begin(), indices.end());
 
      // find if the indices specify exists!
      cytnx_int64 b = -1;
      for (cytnx_uint64 i = 0; i < this->_inner_to_outer_idx.size(); i++) {
        if (inds == this->_inner_to_outer_idx[i]) {
          b = i;
          break;
        }
      }
 
      if (b < 0) {
        cytnx_error_msg(
          true, "[ERROR][put_block][BlockFermionicUniTensor] no avaliable block exists.%s", "\n");
      } else {
        cytnx_error_msg(
          in.shape() != this->_blocks[b].shape(),
          "[ERROR][BlockFermionicUniTensor] the shape of input tensor does not match the shape "
          "of block @ idx=%d\n",
          b);
        this->_blocks[b] = in;
      }
    }
 
    void tag_() {
      // no-use!
    }
 
    boost::intrusive_ptr<UniTensor_base> Conj() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Conj_();
      return out;
    }
 
    void Conj_() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      for (int i = 0; i < this->_blocks.size(); i++) {
        this->_blocks[i].Conj_();
      }
    };
 
    // Transpose(_) changes the index order without sign flips and reverses the bond directions for
    // fermionic tensors
    void Transpose_();
    boost::intrusive_ptr<UniTensor_base> Transpose() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Transpose_();
      return out;
    }
 
    void normalize_();
    boost::intrusive_ptr<UniTensor_base> normalize() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->normalize_();
      return out;
    }
 
    boost::intrusive_ptr<UniTensor_base> Dagger() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      boost::intrusive_ptr<UniTensor_base> out = this->Conj();
      out->Transpose_();
      return out;
    }
    void Dagger_() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      this->Conj_();
      this->Transpose_();
    }
 
    void Trace_(const std::string &a, const std::string &b);
    void Trace_(const cytnx_int64 &a, const cytnx_int64 &b);
 
    boost::intrusive_ptr<UniTensor_base> Trace(const std::string &a, const std::string &b) {
      //[21 Aug 2024] This is a copy from BlockUniTensor; the tensor type was adapted
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Trace_(a, b);
      return out;
    }
    boost::intrusive_ptr<UniTensor_base> Trace(const cytnx_int64 &a, const cytnx_int64 &b) {
      //[21 Aug 2024] This is a copy from BlockUniTensor; the tensor type was adapted
      boost::intrusive_ptr<UniTensor_base> out = this->clone();
      out->Trace_(a, b);
      return out;
    }
 
    Tensor Norm() const;
 
    bool elem_exists(const std::vector<cytnx_uint64> &locator) const;
 
    const Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator) const;
    const cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                          const cytnx_complex128 &aux) const;
    const cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                         const cytnx_complex64 &aux) const;
    const cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_double &aux) const;
    const cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_float &aux) const;
    const cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint64 &aux) const;
    const cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int64 &aux) const;
    const cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint32 &aux) const;
    const cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int32 &aux) const;
    const cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                      const cytnx_uint16 &aux) const;
    const cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                     const cytnx_int16 &aux) const;
 
    Scalar::Sproxy at_for_sparse(const std::vector<cytnx_uint64> &locator);
    cytnx_complex128 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                    const cytnx_complex128 &aux);
    cytnx_complex64 &at_for_sparse(const std::vector<cytnx_uint64> &locator,
                                   const cytnx_complex64 &aux);
    cytnx_double &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_double &aux);
    cytnx_float &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_float &aux);
    cytnx_uint64 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint64 &aux);
    cytnx_int64 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int64 &aux);
    cytnx_uint32 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint32 &aux);
    cytnx_int32 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int32 &aux);
    cytnx_uint16 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_uint16 &aux);
    cytnx_int16 &at_for_sparse(const std::vector<cytnx_uint64> &locator, const cytnx_int16 &aux);
 
    void _save_dispatch(std::fstream &f) const;
    void _load_dispatch(std::fstream &f, unsigned int version);
 
    // this will remove the [q_index]-th qnum at [bond_idx]-th Bond!
    void truncate_(const std::string &label, const cytnx_uint64 &q_index);
    void truncate_(const cytnx_int64 &bond_idx, const cytnx_uint64 &q_index);
 
    void Add_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Add_(const Scalar &rhs) {
      cytnx_error_msg(true,
                      "[ERROR] Cannot perform elementwise arithmetic '+' between Scalar and "
                      "BlockFermionicUniTensor.\n %s "
                      "\n",
                      "This operation would destroy the block structure. [Suggest] Avoid or use "
                      "get/put_block(s) to do operation on blocks.");
    }
 
    void Mul_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Mul_(const Scalar &rhs);
 
    void Sub_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Sub_(const Scalar &rhs) {
      cytnx_error_msg(true,
                      "[ERROR] Cannot perform elementwise arithmetic '-' between Scalar and "
                      "BlockFermionicUniTensor.\n %s "
                      "\n",
                      "This operation would destroy the block structure. [Suggest] Avoid or use "
                      "get/put_block(s) to do operation on blocks.");
    }
    void lSub_(const Scalar &lhs) {
      cytnx_error_msg(true,
                      "[ERROR] Cannot perform elementwise arithmetic '-' between Scalar and "
                      "BlockFermionicUniTensor.\n %s "
                      "\n",
                      "This operation would destroy the block structure. [Suggest] Avoid or use "
                      "get/put_block(s) to do operation on blocks.");
    }
 
    void Div_(const boost::intrusive_ptr<UniTensor_base> &rhs);
    void Div_(const Scalar &rhs);
    void lDiv_(const Scalar &lhs) {
      cytnx_error_msg(true,
                      "[ERROR] Cannot perform elementwise arithmetic '/' between Scalar and "
                      "BlockFermionicUniTensor.\n %s "
                      "\n",
                      "This operation would cause division by zero on non-block elements. "
                      "[Suggest] Avoid or use get/put_block(s) to do operation on blocks.");
    }
    void from_(const boost::intrusive_ptr<UniTensor_base> &rhs, bool force, cytnx_double tol = 0.);
 
    void group_basis_();
 
    void combineBond(const std::vector<std::string> &indicators, const bool &force = false);
    void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force = false);
    void combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force,
                      const bool &by_label);
    void combineBonds(const std::vector<std::string> &indicators, const bool &force = false);
 
    const std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        bidx >= this->Nblocks(),
        "[ERROR][BlockFermionicUniTensor] bidx out of bound! only %d blocks in current UTen.\n",
        this->Nblocks());
      return this->_inner_to_outer_idx[bidx];
    }
    std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      cytnx_error_msg(
        bidx >= this->Nblocks(),
        "[ERROR][BlockFermionicUniTensor] bidx out of bound! only %d blocks in current UTen.\n",
        this->Nblocks());
      return this->_inner_to_outer_idx[bidx];
    }
 
    const vec2d<cytnx_uint64> &get_itoi() const {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return this->_inner_to_outer_idx;
    }
    vec2d<cytnx_uint64> &get_itoi() {
      //[21 Aug 2024] This is a copy from BlockUniTensor;
      return this->_inner_to_outer_idx;
    }
 
    void beauty_print_block(std::ostream &os, const cytnx_uint64 &Nin, const cytnx_uint64 &Nout,
                            const std::vector<cytnx_uint64> &qn_indices,
                            const std::vector<Bond> &bonds, const Tensor &block) const;
 
   private:
    // additional information for fermions (#841): if true, the sign of the corresponding
    // block needs to be flipped. Must stay in lockstep with _blocks (one entry per block);
    // outside the class it can only be modified through reset_signflip_()/erase_signflip_()
    // above, each of which preserves that invariant.
    std::vector<cytnx_bool> _signflip;
  };
 
  //======================================================================
 
  class UniTensor_options {
   public:
    bool _is_diag;
    int _dtype;
    int _device;
    int _rowrank;
 
    UniTensor_options() {
      this->_is_diag = false;
      this->_dtype = Type.Double;
      this->_device = Device.cpu;
      this->_rowrank = -1;
    }
 
    UniTensor_options(const UniTensor_options &rhs) {
      this->_is_diag = rhs._is_diag;
      this->_dtype = rhs._dtype;
      this->_device = rhs._device;
      this->_rowrank = rhs._rowrank;
    }
 
    UniTensor_options &operator=(const UniTensor_options &rhs) {
      this->_is_diag = rhs._is_diag;
      this->_dtype = rhs._dtype;
      this->_device = rhs._device;
      this->_rowrank = rhs._rowrank;
      return *this;
    }
 
    UniTensor_options &is_diag(const bool &in) {
      this->_is_diag = in;
      return *this;
    }
    UniTensor_options &dtype(const int &in) {
      this->_dtype = in;
      return *this;
    }
    UniTensor_options &device(const int &in) {
      this->_device = in;
      return *this;
    }
    UniTensor_options &rowrank(const int &in) {
      this->_rowrank = in;
      return *this;
    }
  };
 
  class UniTensor {
   public:
    boost::intrusive_ptr<UniTensor_base> _impl;
    UniTensor() : _impl(new UniTensor_base()){};
    UniTensor(const UniTensor &rhs) { this->_impl = rhs._impl; }
    UniTensor &operator=(const UniTensor &rhs) {
      this->_impl = rhs._impl;
      return *this;
    }
 
 
    explicit UniTensor(const Tensor &in_tensor, const bool &is_diag = false,
                       const cytnx_int64 &rowrank = -1,
                       const std::vector<std::string> &in_labels = {}, const std::string &name = "")
        : _impl(new UniTensor_base()) {
      this->Init(in_tensor, is_diag, rowrank, in_labels, name);
    }
    void Init(const Tensor &in_tensor, const bool &is_diag = false, const cytnx_int64 &rowrank = -1,
              const std::vector<std::string> &in_labels = {}, const std::string &name = "") {
      boost::intrusive_ptr<UniTensor_base> out(new DenseUniTensor());
      out->Init_by_Tensor(in_tensor, is_diag, rowrank, name);
      this->_impl = out;
      if (in_labels.size() != 0) this->set_labels(in_labels);
    }
 
 
    UniTensor(const std::vector<Bond> &bonds, const std::vector<std::string> &in_labels = {},
              const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
              const int &device = Device.cpu, const bool &is_diag = false,
              const std::string &name = "")
        : _impl(new UniTensor_base()) {
  #ifdef UNI_DEBUG
      cytnx_warning_msg(
        true,
        "[DEBUG] message: entry for UniTensor(const std::vector<Bond> &bonds, const "
        "std::vector<std::string> &in_labels={}, const cytnx_int64 &rowrank=-1, const unsigned "
        "int "
        "&dtype=Type.Double, const int &device = Device.cpu, const bool &is_diag=false)%s",
        "\n");
  #endif
      this->Init(bonds, in_labels, rowrank, dtype, device, is_diag, name);
    }
 
    /* [developing]
    void Init(const std::vector<Bond> &bonds, const std::vector<std::string> &in_labels = {},
              const UniTensor_options &UToptions = UniTensor_options(), const std::string &name =
    ""){ this->Init(bonds,in_labels, UToptions._rowrank, UToptions._dtype  , UToptions._device ,
                                    UToptions._is_diag,
                                    name);
    }
    */
 
    UniTensor &Init(const std::vector<Bond> &bonds, const std::vector<std::string> &in_labels = {},
                    const cytnx_int64 &rowrank = -1, const unsigned int &dtype = Type.Double,
                    const int &device = Device.cpu, const bool &is_diag = false,
                    const std::string &name = "") {
      // checking type:
      bool is_sym = false;
      int sym_fver = -1;
      bool fermionic = false;
 
      for (cytnx_uint64 i = 0; i < bonds.size(); i++) {
        // check
        if (bonds[i].syms().size() != 0) {
          is_sym = true;
          if (sym_fver == -1)
            sym_fver = bonds[i]._impl->_degs.size();
          else {
            cytnx_error_msg((bool(sym_fver) ^ bool(bonds[i]._impl->_degs.size())),
                            "[ERROR] When initializing a UniTensor with symmetries, all Bonds must "
                            "be in the same format!%s",
                            "\n");
          }
          if (!fermionic) {
            std::vector<Symmetry> symms = bonds[i].syms();
            for (cytnx_uint64 i = 0; i < symms.size(); i++) {
              if (symms[i].is_fermionic()) fermionic = true;
            }
          }
        } else
          cytnx_error_msg(
            is_sym, "[ERROR] Cannot have bonds with mixing of symmetry and non-symmetry.%s", "\n");
      }
 
      // dynamical dispatch:
      if (is_sym) {
  #ifdef UNI_DEBUG
        cytnx_warning_msg(true, "[DEBUG] message: entry dispatch: UniTensor: symmetric%s", "\n");
  #endif
        cytnx_error_msg(
          sym_fver < 1,
          "[ERROR] Symmetric tensor, but no degeneracies given. The UniTensor seems broken.%s",
          "\n");
        if (fermionic) {
          boost::intrusive_ptr<UniTensor_base> out(new BlockFermionicUniTensor());
          this->_impl = out;
        } else {
          boost::intrusive_ptr<UniTensor_base> out(new BlockUniTensor());
          this->_impl = out;
        }
      } else {
        boost::intrusive_ptr<UniTensor_base> out(new DenseUniTensor());
        this->_impl = out;
      }
      this->_impl->Init(bonds, in_labels, rowrank, dtype, device, is_diag, false, name);
      return *this;
    }
 
    UniTensor &set_name_(const std::string &in) {
      this->_impl->set_name_(in);
      return *this;
    }
 
    [[deprecated("Please use set_name_(const std::string &in) instead.")]] UniTensor &set_name(
      const std::string &in) {
      return this->set_name_(in);
    }
 
    UniTensor &set_label_(const cytnx_int64 &idx, const std::string &new_label) {
      this->_impl->set_label_(idx, new_label);
      return *this;
    }
 
    UniTensor &set_label_(const cytnx_int64 &idx, const char *new_label) {
      this->_impl->set_label_(idx, std::string(new_label));
      return *this;
    }
 
    UniTensor &set_label_(const std::string &old_label, const std::string &new_label) {
      this->_impl->set_label_(old_label, new_label);
      return *this;
    }
 
    UniTensor &set_label_(const char *old_label, const std::string &new_label) {
      this->_impl->set_label_(std::string(old_label), new_label);
      return *this;
    }
 
    UniTensor &set_label_(const std::string &old_label, const char *new_label) {
      this->_impl->set_label_(old_label, std::string(new_label));
      return *this;
    }
 
    UniTensor &set_label_(const char *old_label, const char *new_label) {
      this->_impl->set_label_(std::string(old_label), std::string(new_label));
      return *this;
    }
 
    [[deprecated(
      "Please use set_label_(const cytnx_int64 &idx, const std::string &new_label) "
      "instead.")]] UniTensor &
      set_label(const cytnx_int64 &idx, const std::string &new_label) {
      return this->set_label_(idx, new_label);
    }
 
    [[deprecated(
      "Please use set_label_(const cytnx_int64 &idx, const char *new_label) "
      "instead.")]] UniTensor &
      set_label(const cytnx_int64 &idx, const char *new_label) {
      return this->set_label_(idx, new_label);
    }
 
    [[deprecated(
      "Please use set_label_(const std::string &old_label, const std::string "
      "&new_label) instead.")]] UniTensor &
      set_label(const std::string &old_label, const std::string &new_label) {
      return this->set_label_(old_label, new_label);
    }
 
    [[deprecated(
      "Please use set_label_(const char *old_label, const std::string &new_label) "
      "instead.")]] UniTensor &
      set_label(const char *old_label, const std::string &new_label) {
      return this->set_label_(old_label, new_label);
    }
 
    [[deprecated(
      "Please use set_label_(const std::string &old_label, const char *new_label) "
      "instead.")]] UniTensor &
      set_label(const std::string &old_label, const char *new_label) {
      return this->set_label_(old_label, new_label);
    }
 
    [[deprecated(
      "Please use set_label_(const char *old_label, const char *new_label) "
      "instead.")]] UniTensor &
      set_label(const char *old_label, const char *new_label) {
      return this->set_label_(old_label, new_label);
    }
 
    /*
    UniTensor &change_label(const cytnx_int64 &old_label, const cytnx_int64 &new_label){
        this->_impl->change_label(old_label,new_label);
        return *this;
    }
    */
 
    [[deprecated(
      "Please use "
      "UniTensor &relabel_(const std::vector<std::string> &new_labels) "
      "instead.")]] UniTensor &
      set_labels(const std::vector<std::string> &new_labels) {
      this->_impl->set_labels(new_labels);
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "UniTensor &relabel_(const std::initializer_list<char *> &new_labels) "
      "instead.")]] UniTensor &
      set_labels(const std::initializer_list<char *> &new_labels) {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      this->_impl->set_labels(vs);
      return *this;
    }
 
    UniTensor &set_rowrank_(const cytnx_uint64 &new_rowrank) {
      this->_impl->set_rowrank_(new_rowrank);
      return *this;
    }
 
    UniTensor set_rowrank(const cytnx_uint64 &new_rowrank) const {
      UniTensor out;
      out._impl = this->_impl->set_rowrank(new_rowrank);
      return out;
    }
 
    template <class T>
    T &item() {
      if (this->is_blockform()) {
        cytnx_error_msg(
          !this->is_scalar(),
          "[ERROR] Cannot use item on a non-scalar UniTensor with Symmetry.\n suggestion: use "
          "get_block()/get_blocks() first.%s",
          "\n");
        return this->_impl->get_block_(0).item<T>();
      }
 
      DenseUniTensor *tmp = static_cast<DenseUniTensor *>(this->_impl.get());
      return tmp->_block.item<T>();
    }
 
    Scalar::Sproxy item() const {
      if (this->is_blockform()) {
        cytnx_error_msg(
          !this->is_scalar(),
          "[ERROR] Cannot use item on a non-scalar UniTensor with Symmetry.\n suggestion: use "
          "get_block()/get_blocks() first.%s",
          "\n");
        return this->_impl->get_block_(0).item();
      }
 
      DenseUniTensor *tmp = static_cast<DenseUniTensor *>(this->_impl.get());
      return tmp->_block.item();
    }
    cytnx_uint64 Nblocks() const { return this->_impl->Nblocks(); }
 
    cytnx_uint64 rank() const { return this->_impl->rank(); }
 
    cytnx_uint64 rowrank() const { return this->_impl->rowrank(); }
 
    unsigned int dtype() const { return this->_impl->dtype(); }
 
    int uten_type() const { return this->_impl->uten_type(); }
 
    bool is_void() const { return this->uten_type() == UTenType.Void; }
 
    bool is_scalar() const { return !this->is_void() && this->rank() == 0; }
 
    cytnx_uint64 size() const {
      if (this->is_void()) return 0;
      cytnx_uint64 result = 1;
      for (const cytnx_uint64 extent : this->shape()) result *= extent;
      return result;
    }
 
    bool is_empty() const { return !this->is_void() && this->size() == 0; }
 
    int device() const { return this->_impl->device(); }
 
    std::string name() const { return this->_impl->name(); }
 
    std::string dtype_str() const { return this->_impl->dtype_str(); }
 
    std::string device_str() const { return this->_impl->device_str(); }
 
    std::string uten_type_str() const { return this->_impl->uten_type_str(); }
 
    bool is_contiguous() const { return this->_impl->is_contiguous(); }
 
    bool is_diag() const { return this->_impl->is_diag(); }
 
    bool is_tag() const { return this->_impl->is_tag(); }
 
    std::vector<Symmetry> syms() const { return this->_impl->syms(); }
 
    const bool &is_braket_form() const { return this->_impl->is_braket_form(); }
 
    const std::vector<std::string> &labels() const { return this->_impl->labels(); }
    cytnx_int64 get_index(std::string label) const { return this->_impl->get_index(label); }
 
    const std::vector<Bond> &bonds() const { return this->_impl->bonds(); }
 
    const Bond &bond_(const cytnx_uint64 &idx) const { return this->_impl->bond_(idx); }
    Bond &bond_(const cytnx_uint64 &idx) { return this->_impl->bond_(idx); }
 
    const Bond &bond_(const std::string &label) const { return this->_impl->bond_(label); }
    Bond &bond_(const std::string &label) { return this->_impl->bond_(label); }
 
    Bond bond(const cytnx_uint64 &idx) const { return this->_impl->bond_(idx).clone(); }
    Bond bond(const std::string &label) const { return this->_impl->bond_(label).clone(); }
 
    std::vector<cytnx_uint64> shape() const { return this->_impl->shape(); }
 
    std::vector<bool> signflip() const { return this->_impl->signflip(); }
 
    bool is_blockform() const { return this->_impl->is_blockform(); }
 
    UniTensor &to_(const int &device) {
      this->_impl->to_(device);
      return *this;
    }
 
    UniTensor to(const int &device) const {
      UniTensor out;
      out._impl = this->_impl->to(device);
      return out;
    }
 
    UniTensor clone() const {
      UniTensor out;
      out._impl = this->_impl->clone();
      return out;
    }
 
    UniTensor &relabel_(const std::vector<std::string> &new_labels) {
      this->_impl->relabel_(new_labels);
      return *this;
    }
    [[deprecated(
      "Please use "
      "UniTensor &relabel_(const std::vector<std::string> &new_labels) "
      "instead.")]] UniTensor &
      relabels_(const std::vector<std::string> &new_labels) {
      this->_impl->relabels_(new_labels);
      return *this;
    }
 
    UniTensor relabel(const std::vector<std::string> &new_labels) const {
      UniTensor out;
      out._impl = this->_impl->relabel(new_labels);
      return out;
    }
    [[deprecated(
      "Please use "
      "UniTensor relabel(const std::vector<std::string> &new_labels) const "
      "instead.")]] UniTensor
      relabels(const std::vector<std::string> &new_labels) const {
      UniTensor out;
      out._impl = this->_impl->relabels(new_labels);
      return out;
    }
 
    UniTensor relabel(const std::initializer_list<char *> &new_labels) const {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      UniTensor out;
      out._impl = this->_impl->relabel(vs);
      return out;
    }
    [[deprecated(
      "Please use "
      "UniTensor relabel(const std::initializer_list<char *> &new_labels) const "
      "instead.")]] UniTensor
      relabels(const std::initializer_list<char *> &new_labels) const {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      UniTensor out;
      out._impl = this->_impl->relabels(vs);
      return out;
    }
    UniTensor &relabel_(const std::initializer_list<char *> &new_labels) {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      this->_impl->relabel_(vs);
      return *this;
    }
    [[deprecated(
      "Please use "
      "UniTensor &relabel_(const std::initializer_list<char *> &new_labels) "
      "instead.")]] UniTensor &
      relabels_(const std::initializer_list<char *> &new_labels) {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      this->_impl->relabels_(vs);
      return *this;
    }
 
    UniTensor relabel(const std::vector<std::string> &old_labels,
                      const std::vector<std::string> &new_labels) const {
      UniTensor out;
      out._impl = this->_impl->relabel(old_labels, new_labels);
      return out;
    }
    [[deprecated(
      "Please use "
      "UniTensor relabel(const std::vector<std::string> &old_labels, const "
      "std::vector<std::string> &new_labels) const "
      "instead.")]] UniTensor
      relabels(const std::vector<std::string> &old_labels,
               const std::vector<std::string> &new_labels) const {
      UniTensor out;
      out._impl = this->_impl->relabels(old_labels, new_labels);
      return out;
    }
 
    UniTensor &relabel_(const std::vector<std::string> &old_labels,
                        const std::vector<std::string> &new_labels) {
      this->_impl->relabel_(old_labels, new_labels);
      return *this;
    }
    [[deprecated(
      "Please use "
      "UniTensor &relabel_(const std::vector<std::string> &old_labels, const "
      "std::vector<std::string> &new_labels) "
      "instead.")]] UniTensor &
      relabels_(const std::vector<std::string> &old_labels,
                const std::vector<std::string> &new_labels) {
      this->_impl->relabels_(old_labels, new_labels);
      return *this;
    }
 
    UniTensor relabel(const std::initializer_list<char *> &old_labels,
                      const std::initializer_list<char *> &new_labels) const {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      std::vector<char *> old_lbls(old_labels);
      std::vector<std::string> vs_old(old_lbls.size());
      transform(old_lbls.begin(), old_lbls.end(), vs_old.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      return this->relabel(vs_old, vs);
    }
 
    [[deprecated(
      "Please use "
      "UniTensor relabel(const std::initializer_list<char *> &old_labels, const "
      "std::initializer_list<char *> &new_labels) const "
      "instead.")]] UniTensor
      relabels(const std::initializer_list<char *> &old_labels,
               const std::initializer_list<char *> &new_labels) const {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      std::vector<char *> old_lbls(old_labels);
      std::vector<std::string> vs_old(old_lbls.size());
      transform(old_lbls.begin(), old_lbls.end(), vs_old.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      return this->relabels(vs_old, vs);
    }
 
    UniTensor &relabel_(const std::initializer_list<char *> &old_labels,
                        const std::initializer_list<char *> &new_labels) {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      std::vector<char *> old_lbls(old_labels);
      std::vector<std::string> vs_old(old_lbls.size());
      transform(old_lbls.begin(), old_lbls.end(), vs_old.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      this->relabel_(vs_old, vs);
      return *this;
    }
    [[deprecated(
      "Please use "
      "UniTensor &relabel_(const std::initializer_list<char *> &old_labels, const "
      "std::initializer_list<char *> &new_labels) "
      "instead.")]] UniTensor &
      relabels_(const std::initializer_list<char *> &old_labels,
                const std::initializer_list<char *> &new_labels) {
      std::vector<char *> new_lbls(new_labels);
      std::vector<std::string> vs(new_lbls.size());
      transform(new_lbls.begin(), new_lbls.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      std::vector<char *> old_lbls(old_labels);
      std::vector<std::string> vs_old(old_lbls.size());
      transform(old_lbls.begin(), old_lbls.end(), vs_old.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      this->relabels_(vs_old, vs);
      return *this;
    }
 
    UniTensor relabel(const cytnx_int64 &inx, const std::string &new_label) const {
      UniTensor out;
      out._impl = this->_impl->relabel(inx, new_label);
      return out;
    }
    UniTensor &relabel_(const cytnx_int64 &inx, const std::string &new_label) {
      this->_impl->relabel_(inx, new_label);
      return *this;
    }
 
    UniTensor &relabel_(const std::string &old_label, const std::string &new_label) {
      this->_impl->relabel_(old_label, new_label);
      return *this;
    }
 
    UniTensor relabel(const std::string &old_label, const std::string &new_label) const {
      UniTensor out;
      out._impl = this->_impl->relabel(old_label, new_label);
      return out;
    }
 
    UniTensor astype(const unsigned int &dtype) const {
      UniTensor out;
      if (this->dtype() == dtype) {
        out._impl = this->_impl;
      } else {
        out._impl = this->_impl->astype(dtype);
      }
      return out;
    }
 
    UniTensor permute(const std::vector<cytnx_int64> &mapper,
                      const cytnx_int64 &rowrank = -1) const {
      UniTensor out;
      out._impl = this->_impl->permute(mapper, rowrank);
      return out;
    }
 
    UniTensor permute(const std::vector<std::string> &mapper,
                      const cytnx_int64 &rowrank = -1) const {
      UniTensor out;
      out._impl = this->_impl->permute(mapper, rowrank);
      return out;
    }
 
    UniTensor permute(const std::initializer_list<char *> &mapper,
                      const cytnx_int64 &rowrank = -1) const {
      std::vector<char *> mprs = mapper;
      std::vector<std::string> vs(mprs.size());
      transform(mprs.begin(), mprs.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      return this->permute(vs, rowrank);
    }
 
    UniTensor &permute_(const std::vector<cytnx_int64> &mapper, const cytnx_int64 &rowrank = -1) {
      this->_impl->permute_(mapper, rowrank);
      return *this;
    }
 
    UniTensor &permute_(const std::vector<std::string> &mapper, const cytnx_int64 &rowrank = -1) {
      this->_impl->permute_(mapper, rowrank);
      return *this;
    }
 
    UniTensor permute_nosignflip(const std::vector<cytnx_int64> &mapper,
                                 const cytnx_int64 &rowrank = -1) const {
      UniTensor out;
      out._impl = this->_impl->permute_nosignflip(mapper, rowrank);
      return out;
    }
 
    UniTensor permute_nosignflip(const std::vector<std::string> &mapper,
                                 const cytnx_int64 &rowrank = -1) const {
      UniTensor out;
      out._impl = this->_impl->permute_nosignflip(mapper, rowrank);
      return out;
    }
 
    UniTensor permute_nosignflip(const std::initializer_list<char *> &mapper,
                                 const cytnx_int64 &rowrank = -1) const {
      std::vector<char *> mprs = mapper;
      std::vector<std::string> vs(mprs.size());
      transform(mprs.begin(), mprs.end(), vs.begin(),
                [](char *x) -> std::string { return std::string(x); });
 
      return this->permute_nosignflip(vs, rowrank);
    }
 
    UniTensor &permute_nosignflip_(const std::vector<cytnx_int64> &mapper,
                                   const cytnx_int64 &rowrank = -1) {
      this->_impl->permute_nosignflip_(mapper, rowrank);
      return *this;
    }
 
    UniTensor &permute_nosignflip_(const std::vector<std::string> &mapper,
                                   const cytnx_int64 &rowrank = -1) {
      this->_impl->permute_nosignflip_(mapper, rowrank);
      return *this;
    }
 
    UniTensor twist(const std::string &label) const {
      UniTensor out = this->clone();
      out._impl->twist_(label);
      return out;
    }
    UniTensor twist(const cytnx_int64 &idx) const {
      UniTensor out = this->clone();
      out._impl->twist_(idx);
      return out;
    }
    UniTensor &twist_(const std::string &label) {
      this->_impl->twist_(label);
      return *this;
    }
    UniTensor &twist_(const cytnx_int64 &idx) {
      this->_impl->twist_(idx);
      return *this;
    }
 
    UniTensor fermion_twists() const {
      UniTensor out = this->clone();
      out._impl->fermion_twists_();
      return out;
    }
    UniTensor &fermion_twists_() {
      this->_impl->fermion_twists_();
      return *this;
    }
 
    UniTensor contiguous() const {
      UniTensor out;
      out._impl = this->_impl->contiguous();
      return out;
    }
 
    UniTensor &contiguous_() {
      this->_impl = this->_impl->contiguous_();
      return *this;
    }
 
    UniTensor apply() const {
      UniTensor out;
      out._impl = this->_impl->apply();
      return out;
    }
 
    UniTensor &apply_() {
      this->_impl = this->_impl->apply_();
      return *this;
    }
 
    void print_diagram(const bool &bond_info = false) const {
      this->_impl->print_diagram(bond_info);
    }
 
    void print_blocks(const bool &full_info = true) const { this->_impl->print_blocks(full_info); }
 
    void print_block(const cytnx_int64 &idx, const bool &full_info = true) const {
      this->_impl->print_block(idx, full_info);
    }
 
    UniTensor &group_basis_() {
      this->_impl->group_basis_();
      return *this;
    }
 
    UniTensor group_basis() const {
      UniTensor out = this->clone();
      out.group_basis_();
      return out;
    }
 
    template <class T>
    T &at(const std::vector<cytnx_uint64> &locator) {
      if (this->uten_type() == UTenType.Block || this->uten_type() == UTenType.BlockFermionic) {
        // [NEW] this will not check if it exists, if it is not then error will throw!
        T aux;
        return this->_impl->at_for_sparse(locator, aux);
      } else if (this->uten_type() == UTenType.Dense) {
        return this->get_block_().at<T>(locator);
      } else {
        cytnx_error_msg(this->uten_type() == UTenType.Void,
                        "[ERROR] UniTensor is not initialized and of type Void.%s", "\n");
        cytnx_error_msg(
          this->uten_type() == UTenType.Sparse,
          "[ERROR] SparseUniTensor is deprecated. Use BlockUniTensor or LinOp instead.%s", "\n");
        cytnx_error_msg(true, "[ERROR] UniTensor type '%s' not supported\n",
                        this->uten_type_str().c_str());
      }
    }
 
    template <class T>
    const T &at(const std::vector<cytnx_uint64> &locator) const {
      if (this->uten_type() == UTenType.Block || this->uten_type() == UTenType.BlockFermionic) {
        // [NEW] this will not check if it exists, if it is not then error will throw!
        T aux;
        return this->_impl->at_for_sparse(locator, aux);
      } else if (this->uten_type() == UTenType.Dense) {
        return this->get_block_().at<T>(locator);
      } else {
        cytnx_error_msg(this->uten_type() == UTenType.Void,
                        "[ERROR] UniTensor is not initialized and of type Void.%s", "\n");
        cytnx_error_msg(
          this->uten_type() == UTenType.Sparse,
          "[ERROR] SparseUniTensor is deprecated. Use BlockUniTensor or LinOp instead.%s", "\n");
        cytnx_error_msg(true, "[ERROR] UniTensor type '%s' not supported\n",
                        this->uten_type_str().c_str());
      }
    }
 
    template <class T>
    const T &at(const std::vector<std::string> &labels,
                const std::vector<cytnx_uint64> &locator) const {
      // giving label <-> locator one to one corresponding, return the element:
      cytnx_error_msg(locator.size() != labels.size(),
                      "[ERROR][at] length of list should be the same for label and locator.%s",
                      "\n");
      cytnx_error_msg(
        labels.size() != this->rank(),
        "[ERROR][at] length of lists must be the same as UniTensor.rank (# of legs)%s", "\n");
      std::vector<cytnx_uint64> new_locator(this->rank());
      cytnx_uint64 new_loc;
      for (int i = 0; i < labels.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), labels[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR] label:%s does not exist in current UniTensor.\n",
                        labels[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_locator[new_loc] = locator[i];
      }
      return this->at<T>(new_locator);
    }
    template <class T>
    T &at(const std::vector<std::string> &labels, const std::vector<cytnx_uint64> &locator) {
      // giving label <-> locator one to one corresponding, return the element:
      cytnx_error_msg(locator.size() != labels.size(),
                      "[ERROR][at] length of list should be the same for label and locator.%s",
                      "\n");
      cytnx_error_msg(
        labels.size() != this->rank(),
        "[ERROR][at] length of lists must be the same as UniTensor.rank (# of legs)%s", "\n");
      std::vector<cytnx_uint64> new_locator(this->rank());
      cytnx_uint64 new_loc;
      for (int i = 0; i < labels.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), labels[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR] label:%s does not exist in current UniTensor.\n",
                        labels[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_locator[new_loc] = locator[i];
      }
      return this->at<T>(new_locator);
    }
 
    const Scalar::Sproxy at(const std::vector<cytnx_uint64> &locator) const {
      if (this->uten_type() == UTenType.Block || this->uten_type() == UTenType.BlockFermionic) {
        return this->_impl->at_for_sparse(locator);
      } else if (this->uten_type() == UTenType.Dense) {
        return this->get_block_().at(locator);
      } else {
        cytnx_error_msg(this->uten_type() == UTenType.Void,
                        "[ERROR] UniTensor is not initialized and of type Void.%s", "\n");
        cytnx_error_msg(
          this->uten_type() == UTenType.Sparse,
          "[ERROR] SparseUniTensor is deprecated. Use BlockUniTensor or LinOp instead.%s", "\n");
        cytnx_error_msg(true, "[ERROR] UniTensor type '%s' not supported\n",
                        this->uten_type_str().c_str());
      }
    }
 
    Scalar::Sproxy at(const std::vector<cytnx_uint64> &locator) {
      if (this->uten_type() == UTenType.Block || this->uten_type() == UTenType.BlockFermionic) {
        return this->_impl->at_for_sparse(locator);
      } else if (this->uten_type() == UTenType.Dense) {
        return this->get_block_().at(locator);
      } else {
        cytnx_error_msg(this->uten_type() == UTenType.Void,
                        "[ERROR] UniTensor is not initialized and of type Void.%s", "\n");
        cytnx_error_msg(
          this->uten_type() == UTenType.Sparse,
          "[ERROR] SparseUniTensor is deprecated. Use BlockUniTensor or LinOp instead.%s", "\n");
        cytnx_error_msg(true, "[ERROR] UniTensor type '%s' not supported\n",
                        this->uten_type_str().c_str());
      }
    }
 
    Scalar::Sproxy at(const std::vector<std::string> &labels,
                      const std::vector<cytnx_uint64> &locator) {
      // giving label <-> locator one to one corresponding, return the element:
      cytnx_error_msg(locator.size() != labels.size(),
                      "[ERROR][at] length of list should be the same for label and locator.%s",
                      "\n");
      cytnx_error_msg(
        labels.size() != this->rank(),
        "[ERROR][at] length of lists must be the same as UniTensor.rank (# of legs)%s", "\n");
      std::vector<cytnx_uint64> new_locator(this->rank());
      cytnx_uint64 new_loc;
      for (int i = 0; i < labels.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), labels[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR] label:%s does not exist in current UniTensor.\n",
                        labels[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_locator[new_loc] = locator[i];
      }
      return this->at(new_locator);
    }
 
    const Scalar::Sproxy at(const std::vector<std::string> &labels,
                            const std::vector<cytnx_uint64> &locator) const {
      // giving label <-> locator one to one corresponding, return the element:
      cytnx_error_msg(locator.size() != labels.size(),
                      "[ERROR][at] length of list should be the same for label and locator.%s",
                      "\n");
      cytnx_error_msg(
        labels.size() != this->rank(),
        "[ERROR][at] length of lists must be the same as UniTensor.rank (# of legs)%s", "\n");
      std::vector<cytnx_uint64> new_locator(this->rank());
      cytnx_uint64 new_loc;
      for (int i = 0; i < labels.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), labels[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR] label:%s does not exist in current UniTensor.\n",
                        labels[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_locator[new_loc] = locator[i];
      }
      return this->at(new_locator);
    }
 
    // return a copy of the block with the given block index
    Tensor get_block(const cytnx_uint64 &idx = 0) const { return this->_impl->get_block(idx); };
    //================================
    // return a clone of block
    Tensor get_block(const std::vector<cytnx_int64> &qidx, const bool &force = false) const {
      return this->_impl->get_block(qidx, force);
    }
 
    Tensor get_block(const std::vector<std::string> &labels, const std::vector<cytnx_int64> &qidx,
                     const bool &force = false) const {
      cytnx_error_msg(
        labels.size() != qidx.size(),
        "[ERROR][get_block] length of lists must be the same for both lables and qnidices%s", "\n");
      cytnx_error_msg(labels.size() != this->rank(),
                      "[ERROR][get_block] length of lists must be the rank (# of legs)%s", "\n");
 
      std::vector<cytnx_int64> loc_id(this->rank());
      std::vector<cytnx_int64> new_qidx(this->rank());
 
      cytnx_uint64 new_loc;
      std::vector<cytnx_uint64> new_order(this->rank());
      for (int i = 0; i < labels.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), labels[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR][get_block] label:%s does not exists in current Tensor.\n",
                        labels[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_qidx[new_loc] = qidx[i];
        new_order[i] = new_loc;
      }
      auto out = this->_impl->get_block(new_qidx, force);
      if (out.dtype() != Type.Void) out.permute_(new_order);
      return out;
    }
 
    Tensor get_block(const std::initializer_list<cytnx_int64> &qnum,
                     const bool &force = false) const {
      std::vector<cytnx_int64> tmp = qnum;
      return get_block(tmp, force);
    }
 
    Tensor get_block(const std::vector<cytnx_uint64> &qnum, const bool &force = false) const {
      std::vector<cytnx_int64> iqnum(qnum.begin(), qnum.end());
      return this->_impl->get_block(iqnum, force);
    }
 
    Tensor get_block(const std::vector<std::string> &labels, const std::vector<cytnx_uint64> &qidx,
                     const bool &force = false) const {
      std::vector<cytnx_int64> iqnum(qidx.begin(), qidx.end());
      return this->get_block(labels, iqnum, force);
    }
 
    const Tensor &get_block_(const cytnx_uint64 &idx = 0) const {
      return this->_impl->get_block_(idx);
    }
 
    Tensor &get_block_(const cytnx_uint64 &idx = 0) { return this->_impl->get_block_(idx); }
 
    Tensor &get_block_(const std::vector<cytnx_int64> &qidx, const bool &force = false) {
      return this->_impl->get_block_(qidx, force);
    }
 
    // developer note: Tensor is not the same object (Thus Tensor instead of Tensor& ),
    //                 since we permute! but they have shared data memory.
    Tensor get_block_(const std::vector<std::string> &labels, const std::vector<cytnx_int64> &qidx,
                      const bool &force = false) {
      cytnx_error_msg(
        labels.size() != qidx.size(),
        "[ERROR][get_block] length of lists must be the same for both lables and qnidices%s", "\n");
      cytnx_error_msg(labels.size() != this->rank(),
                      "[ERROR][get_block] length of lists must be the rank (# of legs)%s", "\n");
 
      std::vector<cytnx_int64> loc_id(this->rank());
      std::vector<cytnx_int64> new_qidx(this->rank());
 
      cytnx_uint64 new_loc;
      std::vector<cytnx_uint64> new_order(this->rank());
      for (int i = 0; i < labels.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), labels[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR][get_block] label:%s does not exists in current Tensor.\n",
                        labels[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_qidx[new_loc] = qidx[i];
        new_order[i] = new_loc;
      }
      auto out = this->_impl->get_block_(new_qidx, force);
      if (out.dtype() != Type.Void) {
        out = out.permute(new_order);
      }
      return out;
    }
 
    Tensor &get_block_(const std::initializer_list<cytnx_int64> &qidx, const bool &force = false) {
      std::vector<cytnx_int64> tmp = qidx;
      return get_block_(tmp, force);
    }
 
    Tensor &get_block_(const std::vector<cytnx_uint64> &qidx, const bool &force = false) {
      std::vector<cytnx_int64> iqidx(qidx.begin(), qidx.end());
      return get_block_(iqidx, force);
    }
 
    Tensor get_block_(const std::vector<std::string> &labels, const std::vector<cytnx_uint64> &qidx,
                      const bool &force = false) {
      std::vector<cytnx_int64> iqidx(qidx.begin(), qidx.end());
      return get_block_(labels, iqidx, force);
    }
    //================================
 
    // this only work for non-symm tensor. return a shared view of block
    const Tensor &get_block_(const std::vector<cytnx_int64> &qidx,
                             const bool &force = false) const {
      return this->_impl->get_block_(qidx, force);
    }
 
    const Tensor &get_block_(const std::initializer_list<cytnx_int64> &qidx,
                             const bool &force = false) const {
      std::vector<cytnx_int64> tmp = qidx;
      return this->_impl->get_block_(tmp, force);
    }
 
    const Tensor &get_block_(const std::vector<cytnx_uint64> &qidx,
                             const bool &force = false) const {
      std::vector<cytnx_int64> iqidx(qidx.begin(), qidx.end());
      return get_block_(iqidx, force);
    }
 
    //================================
    //[dev]
    std::vector<Tensor> get_blocks() const { return this->_impl->get_blocks(); }
 
    //[dev]
    const std::vector<Tensor> &get_blocks_(const bool &silent = false) const {
      return this->_impl->get_blocks_(silent);
    }
 
    //[dev]
    std::vector<Tensor> &get_blocks_(const bool &silent = false) {
      return this->_impl->get_blocks_(silent);
    }
 
    UniTensor &put_block(const Tensor &in, const cytnx_uint64 &idx = 0) {
      this->_impl->put_block(in, idx);
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "UniTensor &put_block(const Tensor &in, const cytnx_uint64 &idx) "
      "instead.")]] UniTensor &
      put_block(const Tensor &in, const cytnx_uint64 &idx, const bool &force) {
      this->_impl->put_block(in, idx);
      return *this;
    }
 
    UniTensor &put_block(const Tensor &in_tens, const std::vector<cytnx_int64> &qidx) {
      this->_impl->put_block(in_tens, qidx);
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "UniTensor &put_block(const Tensor &in_tens, const std::vector<cytnx_int64> &qidx) "
      "instead.")]] UniTensor &
      put_block(const Tensor &in_tens, const std::vector<cytnx_int64> &qidx, const bool &force) {
      this->_impl->put_block(in_tens, qidx);
      return *this;
    }
 
    UniTensor &put_block(Tensor &in, const std::vector<std::string> &lbls,
                         const std::vector<cytnx_int64> &qidx) {
      cytnx_error_msg(
        lbls.size() != qidx.size(),
        "[ERROR][put_block] length of lists must be the same for both lables and qnidices%s", "\n");
      cytnx_error_msg(lbls.size() != this->rank(),
                      "[ERROR][put_block] length of lists must be the rank (# of legs)%s", "\n");
 
      std::vector<cytnx_int64> loc_id(this->rank());
      std::vector<cytnx_int64> new_qidx(this->rank());
 
      cytnx_uint64 new_loc;
      // std::vector<cytnx_uint64> new_order(this->rank());
      std::vector<cytnx_uint64> inv_order(this->rank());
      for (int i = 0; i < lbls.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), lbls[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR][put_block] label:%s does not exists in current Tensor.\n",
                        lbls[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_qidx[new_loc] = qidx[i];
        // new_order[i] = new_loc;
        inv_order[new_loc] = i;
      }
      this->_impl->put_block(in.permute(inv_order), new_qidx);
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "UniTensor &put_block(Tensor &in, const std::vector<std::string> &lbls, const "
      "std::vector<cytnx_int64> &qidx) "
      "instead.")]] UniTensor &
      put_block(Tensor &in, const std::vector<std::string> &lbls,
                const std::vector<cytnx_int64> &qidx, const bool &force) {
      return put_block(in, lbls, qidx);
    }
 
    UniTensor &put_block_(Tensor &in, const cytnx_uint64 &idx = 0) {
      this->_impl->put_block_(in, idx);
      return *this;
    }
 
    UniTensor &put_block_(Tensor &in, const std::vector<cytnx_int64> &qidx) {
      this->_impl->put_block_(in, qidx);
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "UniTensor &put_block_(Tensor &in, const std::vector<cytnx_int64> &qidx) "
      "instead.")]] UniTensor &
      put_block_(Tensor &in, const std::vector<cytnx_int64> &qidx, const bool &force) {
      this->_impl->put_block_(in, qidx);
      return *this;
    }
 
    UniTensor &put_block_(Tensor &in, const std::vector<std::string> &lbls,
                          const std::vector<cytnx_int64> &qidx) {
      cytnx_error_msg(
        lbls.size() != qidx.size(),
        "[ERROR][put_block_] length of lists must be the same for both lables and qnidices%s",
        "\n");
      cytnx_error_msg(lbls.size() != this->rank(),
                      "[ERROR][put_block_] length of lists must be the rank (# of legs)%s", "\n");
 
      std::vector<cytnx_int64> loc_id(this->rank());
      std::vector<cytnx_int64> new_qidx(this->rank());
 
      cytnx_uint64 new_loc;
      std::vector<cytnx_uint64> new_order(this->rank());
      std::vector<cytnx_uint64> inv_order(this->rank());
      for (int i = 0; i < lbls.size(); i++) {
        auto res = std::find(this->_impl->_labels.begin(), this->_impl->_labels.end(), lbls[i]);
        cytnx_error_msg(res == this->_impl->_labels.end(),
                        "[ERROR][put_block_] label:%s does not exists in current Tensor.\n",
                        lbls[i].c_str());
        new_loc = std::distance(this->_impl->_labels.begin(), res);
        new_qidx[new_loc] = qidx[i];
        new_order[i] = new_loc;
        inv_order[new_loc] = i;
      }
      in.permute_(inv_order);
      this->_impl->put_block_(in, new_qidx);
      in.permute_(new_order);
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "UniTensor &put_block_(Tensor &in, const std::vector<std::string> &lbls, const "
      "std::vector<cytnx_int64> &qidx) "
      "instead.")]] UniTensor &
      put_block_(Tensor &in, const std::vector<std::string> &lbls,
                 const std::vector<cytnx_int64> &qidx, const bool &force) {
      return put_block_(in, lbls, qidx);
    }
 
    UniTensor get(const std::vector<Accessor> &accessors) const {
      UniTensor out;
      out._impl = this->_impl->get(accessors);
      return out;
    }
 
    UniTensor operator[](const std::vector<cytnx::Accessor> &accessors) const {
      UniTensor out;
      out._impl = this->_impl->get(accessors);
      return out;
    }
    UniTensor operator[](const std::initializer_list<cytnx::Accessor> &accessors) const {
      std::vector<cytnx::Accessor> acc_in = accessors;
      return this->get(acc_in);
    }
    UniTensor operator[](const std::vector<cytnx_int64> &accessors) const {
      std::vector<cytnx::Accessor> acc_in;
      for (cytnx_int64 i = 0; i < accessors.size(); i++) {
        acc_in.push_back(cytnx::Accessor(accessors[i]));
      }
      return this->get(acc_in);
    }
    UniTensor operator[](const std::initializer_list<cytnx_int64> &accessors) const {
      std::vector<cytnx_int64> acc_in = accessors;
      return (*this)[acc_in];
    }
 
    UniTensor &set(const std::vector<Accessor> &accessors, const Tensor &rhs) {
      this->_impl->set(accessors, rhs);
      return *this;
    }
    UniTensor &set(const std::vector<Accessor> &accessors, const UniTensor &rhs) {
      cytnx_error_msg(
        rhs.uten_type() != UTenType.Dense,
        "[ERROR] Cannot set elements from UniTensor with symmetry. Use at() instead.%s", "\n");
      cytnx_error_msg(this->is_diag(), "[ERROR] Cannot set UniTensor with is_diag=True.%s", "\n");
      cytnx_error_msg(rhs.is_diag(),
                      "[ERROR] Cannot set UniTensor. incoming UniTensor is_diag=True.%s", "\n");
 
      this->_impl->set(accessors, rhs.get_block());
      return *this;
    }
 
    UniTensor reshape(const std::vector<cytnx_int64> &new_shape, const cytnx_uint64 &rowrank = 0) {
      UniTensor out;
      out._impl = this->_impl->reshape(new_shape, rowrank);
      return out;
    }
 
    UniTensor &reshape_(const std::vector<cytnx_int64> &new_shape,
                        const cytnx_uint64 &rowrank = 0) {
      this->_impl->reshape_(new_shape, rowrank);
      return *this;
    }
 
    UniTensor to_dense() {
      UniTensor out;
      out._impl = this->_impl->to_dense();
      return out;
    }
 
    UniTensor &to_dense_() {
      this->_impl->to_dense_();
      return *this;
    }
 
    [[deprecated(
      "Please use "
      "combineBond_(const std::vector<std::string> &indicators, const bool &force) "
      "instead.")]] void
      combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force,
                   const bool &by_label) {
      this->_impl->combineBonds(indicators, force, by_label);
    }
 
    [[deprecated(
      "Please use "
      "combineBond_(const std::vector<std::string> &indicators, const bool &force) "
      "instead.")]] void
      combineBonds(const std::vector<std::string> &indicators, const bool &force = false) {
      this->_impl->combineBonds(indicators, force);
    }
 
    [[deprecated(
      "Please use "
      "combineBond_(const std::vector<std::string> &indicators, const bool &force) "
      "instead.")]] void
      combineBonds(const std::vector<cytnx_int64> &indicators, const bool &force = false) {
      this->_impl->combineBonds(indicators, force);
    }
 
    UniTensor &combineBond_(const std::vector<std::string> &indicators, const bool &force = false) {
      this->_impl->combineBond(indicators, force);
      return *this;
    }
 
    UniTensor combineBond(const std::vector<std::string> &indicators,
                          const bool &force = false) const {
      UniTensor out = this->clone();
      out.combineBond_(indicators, force);
      return out;
    }
 
    UniTensor contract(const UniTensor &inR, const bool &mv_elem_self = false,
                       const bool &mv_elem_rhs = false) const {
      UniTensor out;
      out._impl = this->_impl->contract(inR._impl, mv_elem_self, mv_elem_rhs);
      return out;
    }
 
 
    std::vector<Bond> getTotalQnums(const bool physical = false) const {
      return this->_impl->getTotalQnums(physical);
    }
 
    std::vector<std::vector<cytnx_int64>> get_blocks_qnums() const {
      return this->_impl->get_blocks_qnums();
    }
 
    bool same_data(const UniTensor &rhs) const {
      // check same type:
      if (this->_impl->uten_type() != rhs._impl->uten_type()) return false;
 
      return this->_impl->same_data(rhs._impl);
    }
 
    UniTensor &Add_(const UniTensor &rhs) {
      this->_impl->Add_(rhs._impl);
      return *this;
    }
 
    UniTensor &Mul_(const UniTensor &rhs) {
      this->_impl->Mul_(rhs._impl);
      return *this;
    }
 
    UniTensor &Sub_(const UniTensor &rhs) {
      this->_impl->Sub_(rhs._impl);
      return *this;
    }
 
    UniTensor &Div_(const UniTensor &rhs) {
      this->_impl->Div_(rhs._impl);
      return *this;
    }
 
    UniTensor &Add_(const Scalar &rhs) {
      this->_impl->Add_(rhs);
      return *this;
    }
 
    UniTensor &Mul_(const Scalar &rhs) {
      this->_impl->Mul_(rhs);
      return *this;
    }
 
    UniTensor &Sub_(const Scalar &rhs) {
      this->_impl->Sub_(rhs);
      return *this;
    }
 
    UniTensor &Div_(const Scalar &rhs) {
      this->_impl->Div_(rhs);
      return *this;
    }
 
    UniTensor Add(const UniTensor &rhs) const;
 
    UniTensor Add(const Scalar &rhs) const;
 
    UniTensor Mul(const UniTensor &rhs) const;
 
    UniTensor Mul(const Scalar &rhs) const;
 
    UniTensor Div(const UniTensor &rhs) const;
 
    UniTensor Div(const Scalar &rhs) const;
 
    UniTensor Sub(const UniTensor &rhs) const;
 
    UniTensor Sub(const Scalar &rhs) const;
 
    [[deprecated("use norm() (returns Scalar) instead")]] Tensor Norm() const {
      return this->_impl->Norm();
    };
 
    Scalar norm() const { return this->_impl->Norm().item(); };
 
    UniTensor &operator+=(const UniTensor &rhs) {
      this->Add_(rhs);
      return *this;
    }
 
    UniTensor &operator-=(const UniTensor &rhs) {
      this->Sub_(rhs);
      return *this;
    }
 
    UniTensor &operator/=(const UniTensor &rhs) {
      this->Div_(rhs);
      return *this;
    }
 
    UniTensor &operator*=(const UniTensor &rhs) {
      this->Mul_(rhs);
      return *this;
    }
 
    UniTensor &operator+=(const Scalar &rhs) {
      this->Add_(rhs);
      return *this;
    }
 
    UniTensor &operator-=(const Scalar &rhs) {
      this->Sub_(rhs);
      return *this;
    }
 
    UniTensor &operator/=(const Scalar &rhs) {
      this->Div_(rhs);
      return *this;
    }
 
    UniTensor &operator*=(const Scalar &rhs) {
      this->Mul_(rhs);
      return *this;
    }
 
    UniTensor Inv(double clip = -1.) const;
 
    UniTensor &Inv_(double clip = -1.);
 
    UniTensor Conj() const {
      UniTensor out;
      out._impl = this->_impl->Conj();
      return out;
    }
 
    UniTensor &Conj_() {
      this->_impl->Conj_();
      return *this;
    }
 
    UniTensor Transpose() const {
      UniTensor out;
      out._impl = this->_impl->Transpose();
      return out;
    }
 
    UniTensor &Transpose_() {
      this->_impl->Transpose_();
      return *this;
    }
 
    UniTensor normalize() const {
      UniTensor out;
      out._impl = this->_impl->normalize();
      return out;
    }
 
    UniTensor &normalize_() {
      this->_impl->normalize_();
      return *this;
    }
 
    UniTensor Trace(const std::string &a, const std::string &b) const {
      UniTensor out;
      out._impl = this->_impl->Trace(a, b);
      return out;
    }
 
    UniTensor Trace(const cytnx_int64 &a = 0, const cytnx_int64 &b = 1) const {
      UniTensor out;
      out._impl = this->_impl->Trace(a, b);
      return out;
    }
 
    UniTensor &Trace_(const std::string &a, const std::string &b) {
      this->_impl->Trace_(a, b);
      return *this;
    }
 
    UniTensor &Trace_(const cytnx_int64 &a = 0, const cytnx_int64 &b = 1) {
      this->_impl->Trace_(a, b);
      return *this;
    }
 
    UniTensor Dagger() const {
      UniTensor out;
      out._impl = this->_impl->Dagger();
      return out;
    }
 
    UniTensor &Dagger_() {
      this->_impl->Dagger_();
      return *this;
    }
 
    UniTensor &tag_() {
      this->_impl->tag_();
      return *this;
    }
 
    [[deprecated("Please use tag_() instead.")]] UniTensor &tag() { return this->tag_(); }
 
    UniTensor Pow(const double &p) const;
 
    UniTensor &Pow_(const double &p);
 
    bool elem_exists(const std::vector<cytnx_uint64> &locator) const {
      return this->_impl->elem_exists(locator);
    }
 
    template <class T>
    [[deprecated("Use at() instead.")]] T get_elem(const std::vector<cytnx_uint64> &locator) const {
      return this->at<T>(locator);
    }
 
    template <class T2>
    [[deprecated("Use at() instead.")]] UniTensor &set_elem(
      const std::vector<cytnx_uint64> &locator, const T2 &rc) {
      // cytnx_error_msg(true,"[ERROR] invalid type%s","\n");
      this->at(locator) = rc;
      return *this;
    }
 
    void Save(const std::string &fname) const;
 
    void Save(const char *fname) const;
 
    static UniTensor Load(const std::string &fname);
 
    static UniTensor Load(const char *fname);
 
    UniTensor &truncate_(const std::string &label, const cytnx_uint64 &dim) {
      this->_impl->truncate_(label, dim);
      return *this;
    }
 
    UniTensor &truncate_(const cytnx_int64 &bond_idx, const cytnx_uint64 &dim) {
      this->_impl->truncate_(bond_idx, dim);
      return *this;
    }
 
    UniTensor truncate(const std::string &label, const cytnx_uint64 &dim) const {
      UniTensor out = this->clone();
      out.truncate_(label, dim);
      return out;
    }
 
    UniTensor truncate(const cytnx_int64 &bond_idx, const cytnx_uint64 &dim) const {
      UniTensor out = this->clone();
      out.truncate_(bond_idx, dim);
      return out;
    }
 
    const std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) const {
      return this->_impl->get_qindices(bidx);
    }
    std::vector<cytnx_uint64> &get_qindices(const cytnx_uint64 &bidx) {
      return this->_impl->get_qindices(bidx);
    }
 
    const vec2d<cytnx_uint64> &get_itoi() const { return this->_impl->get_itoi(); }
    vec2d<cytnx_uint64> &get_itoi() { return this->_impl->get_itoi(); }
 
    void _Load(std::fstream &f);
    void _Save(std::fstream &f) const;
 
    UniTensor &convert_from_(const UniTensor &rhs, bool force = false, cytnx_double tol = 0.) {
      this->_impl->from_(rhs._impl, force, tol);
      return *this;
    }
 
    [[deprecated(
      "Please use convert_from_(const UniTensor &rhs, bool force, cytnx_double tol) "
      "instead.")]] UniTensor &
      convert_from(const UniTensor &rhs, bool force = false, cytnx_double tol = 0.) {
      return this->convert_from_(rhs, force, tol);
    }
 
    // Generators:
    static UniTensor zeros(const std::vector<cytnx_uint64> &shape,
                           const std::vector<std::string> &in_labels = {},
                           unsigned int dtype = Type.Double, int device = Device.cpu,
                           const std::string &name = "") {
      return UniTensor(cytnx::zeros(shape, dtype, device), false, -1, in_labels, name);
    }
    static UniTensor zeros(std::initializer_list<cytnx_uint64> shape,
                           const std::vector<std::string> &in_labels = {},
                           unsigned int dtype = Type.Double, int device = Device.cpu,
                           const std::string &name = "") {
      return zeros(std::vector<cytnx_uint64>(shape), in_labels, dtype, device, name);
    }
 
    static UniTensor identity(cytnx_uint64 dim, const std::vector<std::string> &in_labels = {},
                              cytnx_bool is_diag = false, unsigned int dtype = Type.Double,
                              int device = Device.cpu, const std::string &name = "") {
      if (is_diag) {
        return UniTensor(cytnx::ones({dim}, dtype, device), is_diag, -1, in_labels, name);
      } else {
        return UniTensor(cytnx::identity(dim, dtype, device), is_diag, -1, in_labels, name);
      }
    }
 
    static UniTensor eye(cytnx_uint64 dim, const std::vector<std::string> &in_labels = {},
                         cytnx_bool is_diag = false, unsigned int dtype = Type.Double,
                         int device = Device.cpu, const std::string &name = "") {
      return identity(dim, in_labels, is_diag, dtype, device, name);
    }
 
    static UniTensor ones(const std::vector<cytnx_uint64> &shape,
                          const std::vector<std::string> &in_labels = {},
                          unsigned int dtype = Type.Double, int device = Device.cpu,
                          const std::string &name = "") {
      return UniTensor(cytnx::ones(shape, dtype, device), false, -1, in_labels, name);
    }
    static UniTensor ones(std::initializer_list<cytnx_uint64> shape,
                          const std::vector<std::string> &in_labels = {},
                          unsigned int dtype = Type.Double, int device = Device.cpu,
                          const std::string &name = "") {
      return ones(std::vector<cytnx_uint64>(shape), in_labels, dtype, device, name);
    }
 
    static UniTensor arange(cytnx_int64 Nelem, const std::vector<std::string> &in_labels = {},
                            const std::string &name = "") {
      return UniTensor(cytnx::arange(Nelem), false, -1, in_labels, name);
    }
 
    static UniTensor arange(cytnx_double start, cytnx_double end, cytnx_double step = 1,
                            const std::vector<std::string> &in_labels = {},
                            unsigned int dtype = Type.Double, int device = Device.cpu,
                            const std::string &name = "") {
      return UniTensor(cytnx::arange(start, end, step, dtype, device), false, -1, in_labels, name);
    }
 
    static UniTensor linspace(cytnx_double start, cytnx_double end, cytnx_uint64 Nelem,
                              bool endpoint = true, const std::vector<std::string> &in_labels = {},
                              unsigned int dtype = Type.Double, int device = Device.cpu,
                              const std::string &name = "") {
      return UniTensor(cytnx::linspace(start, end, Nelem, endpoint, dtype, device), false, -1,
                       in_labels, name);
    }
 
    // Random Generators:
    static UniTensor normal(const cytnx_uint64 &Nelem, const double &mean, const double &std,
                            const std::vector<std::string> &in_labels = {},
                            const unsigned int &seed = cytnx::random::__static_random_device(),
                            const unsigned int &dtype = Type.Double, const int &device = Device.cpu,
                            const std::string &name = "");
 
    static UniTensor normal(const std::vector<cytnx_uint64> &shape, const double &mean,
                            const double &std, const std::vector<std::string> &in_labels = {},
                            const unsigned int &seed = cytnx::random::__static_random_device(),
                            const unsigned int &dtype = Type.Double, const int &device = Device.cpu,
                            const std::string &name = "");
 
    static UniTensor uniform(const cytnx_uint64 &Nelem, const double &low, const double &high,
                             const std::vector<std::string> &in_labels = {},
                             const unsigned int &seed = cytnx::random::__static_random_device(),
                             const unsigned int &dtype = Type.Double,
                             const int &device = Device.cpu, const std::string &name = "");
 
    static UniTensor uniform(const std::vector<cytnx_uint64> &shape, const double &low,
                             const double &high, const std::vector<std::string> &in_labels = {},
                             const unsigned int &seed = cytnx::random::__static_random_device(),
                             const unsigned int &dtype = Type.Double,
                             const int &device = Device.cpu, const std::string &name = "");
 
    // Inplace Random Generators:
    UniTensor &normal_(const double &mean, const double &std,
                       const unsigned int &seed = cytnx::random::__static_random_device());
 
    UniTensor &uniform_(const double &low = 0, const double &high = 1,
                        const unsigned int &seed = cytnx::random::__static_random_device());
 
  };  // class UniTensor
 
  std::ostream &operator<<(std::ostream &os, const UniTensor &in);
 
  UniTensor Contract(const UniTensor &inL, const UniTensor &inR, const bool &cacheL = false,
                     const bool &cacheR = false);
 
  UniTensor Contract(const std::vector<UniTensor> &TNs, const std::string &order,
                     const bool &optimal);
 
  [[deprecated(
    "Please use "
    "UniTensor Contract(const std::vector<UniTensor> &TNs, const std::string &order, const bool "
    "&optimal) "
    "instead.")]] UniTensor
    Contracts(const std::vector<UniTensor> &TNs, const std::string &order, const bool &optimal);
 
  void _resolve_CT(std::vector<UniTensor> &TNlist);
  template <class... T>
  void _resolve_CT(std::vector<UniTensor> &TNlist, const UniTensor &in, const T &...args) {
    TNlist.push_back(in);
    _resolve_CT(TNlist, args...);
  }
 
  template <class... T>
  UniTensor Contract(const UniTensor &in, const T &...args, const std::string &order,
                     const bool &optimal) {
    std::vector<UniTensor> TNlist;
    _resolve_CT(TNlist, in, args...);
    return Contract(TNlist, order, optimal);
  }
 
  template <class... T>
  [[deprecated(
    "Please use "
    "UniTensor Contract(const UniTensor &in, const T &...args, const std::string &order, const "
    "bool &optimal) "
    "instead.")]] UniTensor
    Contracts(const UniTensor &in, const T &...args, const std::string &order,
              const bool &optimal) {
    std::vector<UniTensor> TNlist;
    _resolve_CT(TNlist, in, args...);
    return Contracts(TNlist, order, optimal);
  }
 
}  // namespace cytnx
 
#endif  // BACKEND_TORCH
 
#endif  // CYTNX_UNITENSOR_H_