diff --git a/include/Accessor.hpp b/include/Accessor.hpp
index dbb0fc749..b497f75e2 100644
--- a/include/Accessor.hpp
+++ b/include/Accessor.hpp
@@ -6,8 +6,8 @@
 #include <vector>
 #include <cstring>
 #include <string>
-#include <iostream>
 #include <initializer_list>
+#include <ostream>
 
 namespace cytnx {
 
diff --git a/include/Bond.hpp b/include/Bond.hpp
index 4e2417901..3de7ae725 100644
--- a/include/Bond.hpp
+++ b/include/Bond.hpp
@@ -1,14 +1,16 @@
 #ifndef CYTNX_BOND_H_
 #define CYTNX_BOND_H_
 
-#include "Type.hpp"
-#include "cytnx_error.hpp"
-#include "Symmetry.hpp"
-#include <initializer_list>
-#include <vector>
+#include <algorithm>
 #include <fstream>
+#include <initializer_list>
 #include <map>
-#include <algorithm>
+#include <ostream>
+#include <vector>
+
+#include "Symmetry.hpp"
+#include "Type.hpp"
+#include "cytnx_error.hpp"
 #include "intrusive_ptr_base.hpp"
 #include "utils/vec_clone.hpp"
 
diff --git a/include/Device.hpp b/include/Device.hpp
index 99747e107..418a5e2d0 100644
--- a/include/Device.hpp
+++ b/include/Device.hpp
@@ -1,8 +1,9 @@
 #ifndef CYTNX_DEVICE_H_
 #define CYTNX_DEVICE_H_
 
-#include <vector>
+#include <iostream>
 #include <string>
+#include <vector>
 
 namespace cytnx {
 
@@ -18,7 +19,10 @@ namespace cytnx {
     int Ncpus;
     std::vector<std::vector<bool>> CanAccessPeer;
     Device_class();
-    void Print_Property();
+    void Print_Property() { this->Print_Property(std::cout); }
+    /// @cond
+    void Print_Property(std::ostream &os);
+    /// @endcond
     std::string getname(const int &device_id);
     ~Device_class();
     // void cudaDeviceSynchronize();
diff --git a/include/Gncon.hpp b/include/Gncon.hpp
index d83c1c14f..3b7cf0032 100644
--- a/include/Gncon.hpp
+++ b/include/Gncon.hpp
@@ -1,16 +1,18 @@
 #ifndef CYTNX_GNCON_H_
 #define CYTNX_GNCON_H_
 
-#include "Type.hpp"
-#include "cytnx_error.hpp"
+#include <fstream>
 #include <initializer_list>
-#include <vector>
+#include <iostream>
 #include <map>
-#include <fstream>
-#include "intrusive_ptr_base.hpp"
-#include "utils/utils.hpp"
+#include <vector>
+
+#include "Type.hpp"
 #include "UniTensor.hpp"
 #include "contraction_tree.hpp"
+#include "cytnx_error.hpp"
+#include "intrusive_ptr_base.hpp"
+#include "utils/utils.hpp"
 
 #ifdef BACKEND_TORCH
 #else
@@ -69,7 +71,7 @@ namespace cytnx {
       return true;
     }
 
-    // void print_Gncon() const;
+    // void print_Gncon(std::ostream &os = std::cout) const;
 
     // void PreConstruct(bool force = true);
 
@@ -93,7 +95,7 @@ namespace cytnx {
     virtual void clear();
     virtual std::string getOptimalOrder();
     virtual UniTensor Launch(const bool &optimal = false, const std::string &contract_order = "");
-    virtual void PrintNet(std::ostream &os);
+    virtual void PrintNet(std::ostream &os) const;
     virtual boost::intrusive_ptr<Gncon_base> clone();
     virtual void Savefile(const std::string &fname);
     virtual ~Gncon_base(){};
@@ -137,7 +139,7 @@ namespace cytnx {
       boost::intrusive_ptr<Gncon_base> out(tmp);
       return out;
     }
-    void PrintNet(std::ostream &os);
+    void PrintNet(std::ostream &os) const;
     void Savefile(const std::string &fname);
     ~RegularGncon(){};
   };
@@ -182,7 +184,7 @@ namespace cytnx {
       boost::intrusive_ptr<Gncon_base> out(tmp);
       return out;
     }
-    void PrintNet(std::ostream &os){};
+    void PrintNet(std::ostream &os) const {};
     void Savefile(const std::string &fname){};
     ~FermionGncon(){};
   };
@@ -352,7 +354,10 @@ namespace cytnx {
       out._impl = this->_impl->clone();
       return out;
     }
-    void PrintNet() { this->_impl->PrintNet(std::cout); }
+    void PrintNet() const { this->_impl->PrintNet(std::cout); }
+    /// @cond
+    void PrintNet(std::ostream &os) const { this->_impl->PrintNet(os); }
+    /// @endcond
 
     void Savefile(const std::string &fname) { this->_impl->Savefile(fname); }
   };
diff --git a/include/LinOp.hpp b/include/LinOp.hpp
index 1ed60f8b6..a09d9c297 100644
--- a/include/LinOp.hpp
+++ b/include/LinOp.hpp
@@ -1,17 +1,19 @@
 #ifndef CYTNX_LINOP_H_
 #define CYTNX_LINOP_H_
 
-#include "Type.hpp"
-#include "cytnx_error.hpp"
-#include <vector>
+#include <algorithm>
 #include <fstream>
 #include <functional>
+#include <iostream>
 #include <map>
 #include <utility>
-#include <algorithm>
-#include "intrusive_ptr_base.hpp"
+#include <vector>
+
 #include "Tensor.hpp"
+#include "Type.hpp"
 #include "UniTensor.hpp"
+#include "cytnx_error.hpp"
+#include "intrusive_ptr_base.hpp"
 
 #ifdef BACKEND_TORCH
 #else
@@ -156,7 +158,10 @@ namespace cytnx {
     int dtype() const { return this->_dtype; };
     cytnx_uint64 nx() const { return this->_nx; };
 
-    void _print();
+    void _print() { this->_print(std::cout); };
+    /// @cond
+    void _print(std::ostream &os);
+    /// @endcond
 
     /// @cond
     // this expose to interitance:
diff --git a/include/Network.hpp b/include/Network.hpp
index 75aa82db1..1f9491eeb 100644
--- a/include/Network.hpp
+++ b/include/Network.hpp
@@ -1,16 +1,18 @@
 #ifndef CYTNX_NETWORK_H_
 #define CYTNX_NETWORK_H_
 
-#include "Type.hpp"
-#include "cytnx_error.hpp"
+#include <fstream>
 #include <initializer_list>
-#include <vector>
+#include <iostream>
 #include <map>
-#include <fstream>
-#include "intrusive_ptr_base.hpp"
-#include "utils/utils.hpp"
+#include <vector>
+
+#include "Type.hpp"
 #include "UniTensor.hpp"
 #include "contraction_tree.hpp"
+#include "cytnx_error.hpp"
+#include "intrusive_ptr_base.hpp"
+#include "utils/utils.hpp"
 
 #ifdef BACKEND_TORCH
 #else
@@ -95,7 +97,7 @@ namespace cytnx {
       return true;
     }
 
-    // void print_network() const;
+    // void print_network(std::ostream &os = std::cout) const;
 
     // void PreConstruct(bool force = true);
 
@@ -135,7 +137,7 @@ namespace cytnx {
                            const std::vector<std::vector<std::string>> &labels,
                            const std::vector<std::string> &outlabel, const cytnx_int64 &outrk,
                            const std::string &order, const bool optim);
-    virtual void PrintNet(std::ostream &os);
+    virtual void PrintNet(std::ostream &os = std::cout) const;
     virtual boost::intrusive_ptr<Network_base> clone();
     virtual void Savefile(const std::string &fname);
     virtual ~Network_base(){};
@@ -194,7 +196,7 @@ namespace cytnx {
       boost::intrusive_ptr<Network_base> out(tmp);
       return out;
     }
-    void PrintNet(std::ostream &os);
+    void PrintNet(std::ostream &os = std::cout) const;
     void Savefile(const std::string &fname);
     ~RegularNetwork(){};
   };
@@ -245,7 +247,7 @@ namespace cytnx {
       boost::intrusive_ptr<Network_base> out(tmp);
       return out;
     }
-    void PrintNet(std::ostream &os){};
+    void PrintNet(std::ostream &os = std::cout) const {};
     void Savefile(const std::string &fname){};
     ~FermionNetwork(){};
   };
@@ -450,7 +452,10 @@ namespace cytnx {
       out._impl = this->_impl->clone();
       return out;
     }
-    void PrintNet() { this->_impl->PrintNet(std::cout); }
+    void PrintNet() const { this->_impl->PrintNet(std::cout); }
+    /// @cond
+    void PrintNet(std::ostream &os) const { this->_impl->PrintNet(os); }
+    /// @endcond
 
     void Savefile(const std::string &fname) { this->_impl->Savefile(fname); }
   };
diff --git a/include/Symmetry.hpp b/include/Symmetry.hpp
index ec9fc5b40..3c1a1972d 100644
--- a/include/Symmetry.hpp
+++ b/include/Symmetry.hpp
@@ -2,7 +2,7 @@
 #define CYTNX_SYMMETRY_H_
 
 #include <fstream>
-#include <ostream>
+#include <iostream>
 #include <string>
 #include <vector>
 
@@ -90,7 +90,7 @@ namespace cytnx {
     virtual fermionParity get_fermion_parity(const cytnx_int64 &in_qnum) const;
     virtual bool is_fermionic() const { return false; };
 
-    virtual void print_info() const;
+    virtual void print_info(std::ostream &os = std::cout) const;
     virtual std::string stype_str() const;
     // virtual std::vector<cytnx_int64>& combine_rule(const std::vector<cytnx_int64> &inL, const
     // std::vector<cytnx_int64> &inR);
@@ -118,7 +118,7 @@ namespace cytnx {
     void combine_rule_(cytnx_int64 &out, const cytnx_int64 &inL, const cytnx_int64 &inR,
                        const bool &is_reverse);
     void reverse_rule_(cytnx_int64 &out, const cytnx_int64 &in);
-    void print_info() const;
+    void print_info(std::ostream &os = std::cout) const;
     std::string stype_str() const override { return "U1"; };
   };
   ///@endcond
@@ -144,7 +144,7 @@ namespace cytnx {
     void combine_rule_(cytnx_int64 &out, const cytnx_int64 &inL, const cytnx_int64 &inR,
                        const bool &is_reverse);
     void reverse_rule_(cytnx_int64 &out, const cytnx_int64 &in);
-    void print_info() const;
+    void print_info(std::ostream &os = std::cout) const;
     std::string stype_str() const override { return "Z" + std::to_string(this->n); };
   };
   ///@endcond
@@ -169,7 +169,7 @@ namespace cytnx {
     void reverse_rule_(cytnx_int64 &out, const cytnx_int64 &in);
     fermionParity get_fermion_parity(const cytnx_int64 &in_qnum) const override;
     bool is_fermionic() const override { return true; };
-    void print_info() const;
+    void print_info(std::ostream &os = std::cout) const;
     std::string stype_str() const override { return "fP"; }
   };
   ///@endcond
@@ -194,7 +194,7 @@ namespace cytnx {
     void reverse_rule_(cytnx_int64 &out, const cytnx_int64 &in);
     fermionParity get_fermion_parity(const cytnx_int64 &in_qnum) const override;
     bool is_fermionic() const override { return true; };
-    void print_info() const;
+    void print_info(std::ostream &os = std::cout) const;
     std::string stype_str() const override { return "f#"; }
   };
   ///@endcond
@@ -524,7 +524,10 @@ namespace cytnx {
     /**
      * @brief Print the information of current Symmetry object.
      */
-    void print_info() const { this->_impl->print_info(); }
+    void print_info() const { this->_impl->print_info(std::cout); }
+    /// @cond
+    void print_info(std::ostream &os) const { this->_impl->print_info(os); }
+    /// @endcond
 
     /**
      * @brief the equality operator of the Symmetry object.
diff --git a/include/UniTensor.hpp b/include/UniTensor.hpp
index 899f3a582..a1fafd401 100644
--- a/include/UniTensor.hpp
+++ b/include/UniTensor.hpp
@@ -1,23 +1,24 @@
 #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 <algorithm>
+#include <fstream>
+#include <initializer_list>
 #include <iostream>
-#include <vector>
 #include <map>
+#include <random>
 #include <utility>
-#include <initializer_list>
-#include <fstream>
-#include <algorithm>
-#include "Symmetry.hpp"
+#include <vector>
+
 #include "Bond.hpp"
+#include "Device.hpp"
 #include "Generator.hpp"
-#include <random>
+#include "Symmetry.hpp"
+#include "Tensor.hpp"
+#include "Type.hpp"
+#include "cytnx_error.hpp"
+#include "intrusive_ptr_base.hpp"
+#include "utils/utils.hpp"
 
 #ifdef BACKEND_TORCH
 #else
@@ -273,9 +274,10 @@ namespace cytnx {
     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 void print_diagram(std::ostream &os, const bool &bond_info = false) const;
+    virtual void print_blocks(std::ostream &os, const bool &full_info = true) const;
+    virtual void print_block(std::ostream &os, const cytnx_int64 &idx,
+                             const bool &full_info = true) const;
 
     virtual boost::intrusive_ptr<UniTensor_base> astype(const unsigned int &dtype) const;
 
@@ -636,9 +638,9 @@ namespace cytnx {
       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;
+    void print_diagram(std::ostream &os, const bool &bond_info = false) const;
+    void print_blocks(std::ostream &os, const bool &full_info = true) const;
+    void print_block(std::ostream &os, 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,
@@ -1458,9 +1460,9 @@ namespace cytnx {
       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;
+    void print_diagram(std::ostream &os, const bool &bond_info = false) const;
+    void print_blocks(std::ostream &os, const bool &full_info = true) const;
+    void print_block(std::ostream &os, 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,
@@ -2224,9 +2226,9 @@ namespace cytnx {
     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;
+    void print_diagram(std::ostream &os, const bool &bond_info = false) const;
+    void print_blocks(std::ostream &os, const bool &full_info = true) const;
+    void print_block(std::ostream &os, 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,
@@ -3865,18 +3867,30 @@ namespace cytnx {
     }
 
     /**
-    @brief Plot the diagram of the UniTensor.
+    @brief Print the diagram of the UniTensor.
         @param[in] bond_info whether need to print the information of the bonds of the UniTensor.
     */
     void print_diagram(const bool &bond_info = false) const {
-      this->_impl->print_diagram(bond_info);
+      this->_impl->print_diagram(std::cout, bond_info);
+    }
+    /// @cond
+    void print_diagram(std::ostream &os, const bool &bond_info = false) const {
+      this->_impl->print_diagram(os, bond_info);
     }
+    /// @endcond
 
     /**
     @brief Print all blocks of the UniTensor.
         @param[in] full_info whether need to print the full information of the blocks
     */
-    void print_blocks(const bool &full_info = true) const { this->_impl->print_blocks(full_info); }
+    void print_blocks(const bool &full_info = true) const {
+      this->_impl->print_blocks(std::cout, full_info);
+    }
+    /// @cond
+    void print_blocks(std::ostream &os, const bool &full_info = true) const {
+      this->_impl->print_blocks(os, full_info);
+    }
+    /// @endcond
 
     /**
     @brief Print out the block of the UniTensor with a given block index number.
@@ -3884,8 +3898,13 @@ namespace cytnx {
         @param[in] full_info whether need to print the full information of the block
     */
     void print_block(const cytnx_int64 &idx, const bool &full_info = true) const {
-      this->_impl->print_block(idx, full_info);
+      this->_impl->print_block(std::cout, idx, full_info);
+    }
+    /// @cond
+    void print_block(std::ostream &os, const cytnx_int64 &idx, const bool &full_info = true) const {
+      this->_impl->print_block(os, idx, full_info);
     }
+    /// @endcond
 
     /**
     @brief Group the same quantum number basis together.
diff --git a/include/backend/Scalar.hpp b/include/backend/Scalar.hpp
index e6d1e7a94..19959b25d 100644
--- a/include/backend/Scalar.hpp
+++ b/include/backend/Scalar.hpp
@@ -2,16 +2,19 @@
 #define CYTNX_BACKEND_SCALAR_H_
 
 #ifndef BACKEND_TORCH
-  #include "Type.hpp"
-  #include "cytnx_error.hpp"
-  #include "intrusive_ptr_base.hpp"
-  #include <vector>
+
+  #include <cmath>
   #include <initializer_list>
-  #include <string>
   #include <iostream>
-  #include <cmath>
-  #include <type_traits>
   #include <limits>
+  #include <string>
+  #include <type_traits>
+  #include <vector>
+
+  #include "Type.hpp"
+  #include "cytnx_error.hpp"
+  #include "intrusive_ptr_base.hpp"
+
 namespace cytnx {
 
   ///@cond
@@ -2845,15 +2848,17 @@ namespace cytnx {
      */
     int dtype() const { return this->_impl->_dtype; }
 
-    // print()
     /**
-     * @brief Print the Scalar to the standard output.
+     * @brief Print the Scalar
      */
-    void print() const {
-      this->_impl->print(std::cout);
-      std::cout << std::string(" Scalar dtype: [") << Type.getname(this->_impl->_dtype)
-                << std::string("]") << std::endl;
+    void print() const { this->print(std::cout); }
+    /// @cond
+    void print(std::ostream &os) const {
+      this->_impl->print(os);
+      os << std::string(" Scalar dtype: [") << Type.getname(this->_impl->_dtype) << std::string("]")
+         << std::endl;
     }
+    /// @endcond
 
     // casting
     /// @brief The explicit casting operator of the Scalar class to cytnx::cytnx_double.
diff --git a/include/backend/Storage.hpp b/include/backend/Storage.hpp
index 678a1d2fd..97fa69c3c 100644
--- a/include/backend/Storage.hpp
+++ b/include/backend/Storage.hpp
@@ -63,8 +63,9 @@ namespace cytnx {
 
     void _cpy_bool(void *ptr, const std::vector<cytnx_bool> &vin);
 
-    void print();
-    void print_info();
+    void print(std::ostream &os = std::cout) const;
+    void print_info(std::ostream &os = std::cout) const;
+
     /*
         This function is design to check the type mismatch.
         Handy for developer to exclude the assign of double
@@ -186,7 +187,7 @@ namespace cytnx {
                               const std::vector<cytnx_uint64> &invmapper);
     virtual void PrintElem_byShape(std::ostream &os, const std::vector<cytnx_uint64> &shape,
                                    const std::vector<cytnx_uint64> &mapper = {});
-    virtual void print_elems();
+    virtual void print_elems(std::ostream &os = std::cout) const;
 
     virtual boost::intrusive_ptr<Storage_base> real();
     virtual boost::intrusive_ptr<Storage_base> imag();
@@ -258,7 +259,11 @@ namespace cytnx {
     boost::intrusive_ptr<Storage_base> to(const int &device);
     void PrintElem_byShape(std::ostream &os, const std::vector<cytnx_uint64> &shape,
                            const std::vector<cytnx_uint64> &mapper = {});
-    void print_elems();
+
+    void print_elems() const { this->print_elems(std::cout); }
+    /// @cond
+    void print_elems(std::ostream &os) const;
+    /// @endcond
 
     ~StorageImplementation();
 
@@ -768,10 +773,11 @@ namespace cytnx {
     @brief print the info of the Storage, including the device, dtype and size.
 
     */
-    void print_info() const { this->_impl->print_info(); }
+    void print_info() const { this->_impl->print_info(std::cout); }
     /// @cond
+    void print_info(std::ostream &os) const { this->_impl->print(os); }
     // this is a redundant function
-    void print() const { this->_impl->print(); }
+    void print(std::ostream &os = std::cout) const { this->_impl->print(os); }
     /// @endcond
 
     /**
diff --git a/include/stat.hpp b/include/stat.hpp
index 142bec071..0c2163586 100644
--- a/include/stat.hpp
+++ b/include/stat.hpp
@@ -1,13 +1,13 @@
 #ifndef CYTNX_STAT_H_
 #define CYTNX_STAT_H_
 
-#include "Type.hpp"
-#include "cytnx_error.hpp"
-
-#include "Tensor.hpp"
 #include <algorithm>
 #include <iostream>
 
+#include "Tensor.hpp"
+#include "Type.hpp"
+#include "cytnx_error.hpp"
+
 #ifdef BACKEND_TORCH
 #else
 
@@ -86,7 +86,7 @@ namespace cytnx {
       }
 
       void normalize();
-      void print() const;
+      void print(std::ostream &os = std::cout) const;
 
       const Storage &get_x() const {
         // get x
@@ -179,7 +179,7 @@ namespace cytnx {
       }
 
       void normalize();
-      void print() const;
+      void print(std::ostream &os = std::cout) const;
 
       const Storage &get_x() const {
         // get x
diff --git a/include/utils/print.hpp b/include/utils/print.hpp
index a765da634..c8021d992 100644
--- a/include/utils/print.hpp
+++ b/include/utils/print.hpp
@@ -1,14 +1,14 @@
 #ifndef CYTNX_UTILS_PRINT_H_
 #define CYTNX_UTILS_PRINT_H_
 
-#include "Type.hpp"
-#include "Tensor.hpp"
+#include <iostream>
 
 #include "Bond.hpp"
+#include "Network.hpp"
 #include "Symmetry.hpp"
+#include "Tensor.hpp"
+#include "Type.hpp"
 #include "UniTensor.hpp"
-#include "Network.hpp"
-#include <iostream>
 
 #ifdef BACKEND_TORCH
 #else
@@ -23,6 +23,11 @@ namespace cytnx {
     std::cout << ipt << std::endl;
   }
 
+  template <class T>
+  void print(std::ostream &os, const T &ipt) {
+    os << ipt << std::endl;
+  }
+
 }  // namespace cytnx
 
 #endif  // CYTNX_UTILS_PRINT_H_
diff --git a/pybind/network_py.cpp b/pybind/network_py.cpp
index 5f442d19f..fbecafa37 100644
--- a/pybind/network_py.cpp
+++ b/pybind/network_py.cpp
@@ -114,7 +114,9 @@ void network_binding(py::module &m) {
         return std::string("");
       },
       py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
-    .def("PrintNet", &Network::PrintNet)
+    .def(
+      "PrintNet", [](Network &self) { self.PrintNet(); },
+      py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
     .def_static(
       "Contract",
       [](const std::vector<UniTensor> &utensors, const std::string &Tout,
diff --git a/pybind/scalar_py.cpp b/pybind/scalar_py.cpp
index 6cd648864..005cdc0e6 100644
--- a/pybind/scalar_py.cpp
+++ b/pybind/scalar_py.cpp
@@ -116,7 +116,9 @@ void scalar_binding(py::module &m) {
     .def("iabs", &Scalar::iabs)
     .def("isqrt", &Scalar::isqrt)
     .def("dtype", &Scalar::dtype)
-    .def("print", &Scalar::print)
+    .def(
+      "print", [](const Scalar &self) { self.print(); },
+      py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
 
     // ---- Arithmetic in-place ----
     .def("__iadd__",
diff --git a/pybind/storage_py.cpp b/pybind/storage_py.cpp
index 2808b6672..631dc83c6 100644
--- a/pybind/storage_py.cpp
+++ b/pybind/storage_py.cpp
@@ -1,3 +1,4 @@
+#include <iostream>
 #include <map>
 #include <random>
 #include <string>
@@ -194,8 +195,9 @@ void storage_binding(py::module &m) {
     .def("__deepcopy__", &cytnx::Storage::clone)
     .def("size", &cytnx::Storage::size)
     .def("__len__", [](cytnx::Storage &self) { return self.size(); })
-    .def("print_info", &cytnx::Storage::print_info,
-         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
+    .def(
+      "print_info", [](cytnx::Storage &self) { self.print_info(); },
+      py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
     .def("set_zeros", &cytnx::Storage::set_zeros)
     .def("__eq__",
          [](cytnx::Storage &self, const cytnx::Storage &rhs) -> bool { return self == rhs; })
diff --git a/pybind/unitensor_py.cpp b/pybind/unitensor_py.cpp
index 508a461a4..6309fec02 100644
--- a/pybind/unitensor_py.cpp
+++ b/pybind/unitensor_py.cpp
@@ -588,12 +588,12 @@ void unitensor_binding(py::module &m) {
     .def("apply_", &UniTensor::apply_,
          "Apply fermionic signflips inplacely. Subsequently, signflip() returns False for all "
          "elements. Non-fermionic tensors are left unchanged. See also: apply()")
-    .def("print_diagram", &UniTensor::print_diagram, py::arg("bond_info") = false,
-         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
-    .def("print_blocks", &UniTensor::print_blocks, py::arg("full_info") = true,
-         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
-    .def("print_block", &UniTensor::print_block, py::arg("idx"), py::arg("full_info") = true,
-         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
+    .def("print_diagram", [](UniTensor &self, const bool &bond_info){ self.print_diagram(bond_info); },
+         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>(), py::arg("bond_info") = false)
+    .def("print_block", [](UniTensor &self, const cytnx_int64 &idx, const bool &full_info){ self.print_block(idx, full_info); },
+         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>(), py::arg("idx"), py::arg("full_info") = true)
+    .def("print_blocks", [](UniTensor &self, const bool &full_info){ self.print_blocks(full_info); },
+         py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>(), py::arg("full_info") = true)
 
     .def("group_basis_", &UniTensor::group_basis_)
     .def("group_basis", &UniTensor::group_basis)
diff --git a/src/Accessor.cpp b/src/Accessor.cpp
index f8bb9b6d7..4d9b30f7d 100644
--- a/src/Accessor.cpp
+++ b/src/Accessor.cpp
@@ -1,8 +1,10 @@
 #include "Accessor.hpp"
-#include "utils/str_utils.hpp"
-#include <iostream>
+
 #include <algorithm>
+#include <ostream>
 #include <utility>
+
+#include "utils/str_utils.hpp"
 #include "utils/vec_print.hpp"
 
 using namespace std;
diff --git a/src/BlockFermionicUniTensor.cpp b/src/BlockFermionicUniTensor.cpp
index 5aa99be1d..489eb4462 100644
--- a/src/BlockFermionicUniTensor.cpp
+++ b/src/BlockFermionicUniTensor.cpp
@@ -1,15 +1,20 @@
 #include "UniTensor.hpp"
+
+#include <map>
+#include <ostream>
+#include <stack>
+#include <vector>
+
+#include <boost/unordered_map.hpp>
+
 #include "Accessor.hpp"
-#include "utils/utils.hpp"
-#include "linalg.hpp"
 #include "Generator.hpp"
-#include <vector>
+#include "linalg.hpp"
+#include "utils/utils.hpp"
 #include "utils/vec_clone.hpp"
-#include "utils/vec_print.hpp"
 #include "utils/vec_concatenate.hpp"
-#include <map>
-#include <boost/unordered_map.hpp>
-#include <stack>
+#include "utils/vec_print.hpp"
+
 using namespace std;
 
 #ifdef BACKEND_TORCH
@@ -306,15 +311,14 @@ namespace cytnx {
     os << (std::string(" ") * (3 + Lmax + 5)) << std::string("-") * (4 + mL + mR + 5) << endl;
   }
 
-  void BlockFermionicUniTensor::print_block(const cytnx_int64 &idx, const bool &full_info) const {
+  void BlockFermionicUniTensor::print_block(std::ostream &os, const cytnx_int64 &idx,
+                                            const bool &full_info) const {
     //[21 Aug 2024] This is a copy from BlockUniTensor; sign structure is printed additionally
     cytnx_error_msg(
       (idx < 0) || (idx >= this->_blocks.size()),
       "[ERROR] index [%d] out of bound. should be >0 and < number of available blocks %d\n", idx,
       this->_blocks.size());
 
-    std::ostream &os = std::cout;
-
     os << "========================\n";
     if (this->_is_diag) os << " *is_diag: True\n";
     os << "BLOCK [#" << idx << "]\n";
@@ -356,10 +360,8 @@ namespace cytnx {
     }
   }
 
-  void BlockFermionicUniTensor::print_blocks(const bool &full_info) const {
+  void BlockFermionicUniTensor::print_blocks(std::ostream &os, const bool &full_info) const {
     //[21 Aug 2024] This is a copy from BlockUniTensor; sign structure is printed additionally
-    std::ostream &os = std::cout;
-
     os << "-------- start of print ---------\n";
     char *buffer = (char *)malloc(sizeof(char) * 10240);
     sprintf(buffer, "Tensor name: %s\n", this->_name.c_str());
@@ -391,7 +393,7 @@ namespace cytnx {
 
     // print each blocks with its qnum!
     for (int b = 0; b < this->_blocks.size(); b++) {
-      this->print_block(b, full_info);
+      this->print_block(os, b, full_info);
     }
 
     /*
@@ -418,34 +420,34 @@ namespace cytnx {
     free(buffer);
   }
 
-  void BlockFermionicUniTensor::print_diagram(const bool &bond_info) const {
+  void BlockFermionicUniTensor::print_diagram(std::ostream &os, const bool &bond_info) const {
     //[21 Aug 2024] This is a copy from BlockUniTensor; additionally, 'f' symbols in the corners of
     // the diagram are shown, and the sign structure is printed
     char *buffer = (char *)malloc(10240 * sizeof(char));
     unsigned int BUFFsize = 100;
 
     sprintf(buffer, "--fermionic UniTensor---%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Name : %s\n", this->_name.c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Rank : %d\n", this->_labels.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     // sprintf(buffer, "block_form  : true%s", "\n");
-    // std::cout << std::string(buffer);
+    // os << std::string(buffer);
     sprintf(buffer, "contiguous  : %s\n", this->is_contiguous() ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "valid blocks: %d\n", this->_blocks.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "is diag     : %s\n", this->is_diag() ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "on device   : %s\n", this->device_str().c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "fermion sign:");
     for (int i = 0; i < this->_signflip.size(); i++) {
       sprintf(buffer + strlen(buffer), this->_signflip[i] ? " -1" : " +1");
     }
     sprintf(buffer + strlen(buffer), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     cytnx_uint64 Nin = this->_rowrank;
     cytnx_uint64 Nout = this->_labels.size() - this->_rowrank;
@@ -481,12 +483,12 @@ namespace cytnx {
 
     std::string tmpss;
     sprintf(buffer, "%s row %s col %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s    f%sf    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     for (cytnx_uint64 i = 0; i < vl; i++) {
       sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
 
       if (i < Nin) {
         if (this->_bonds[i].type() == bondType::BD_KET)
@@ -530,21 +532,21 @@ namespace cytnx {
         sprintf(rlbl, "%s", tmpss.c_str());
       }
       sprintf(buffer, "   %s| %s     %s |%s\n", l, llbl, rlbl, r);
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
     }
     sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s    f%sf    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     if (bond_info) {
       for (cytnx_uint64 i = 0; i < this->_bonds.size(); i++) {
         // sprintf(buffer, "lbl:%d ", this->_labels[i]);
         sprintf(buffer, "lbl:%s ", this->_labels[i].c_str());
-        std::cout << std::string(buffer);
-        std::cout << this->_bonds[i] << std::endl;
+        os << std::string(buffer);
+        os << this->_bonds[i] << std::endl;
       }
     }
 
diff --git a/src/BlockUniTensor.cpp b/src/BlockUniTensor.cpp
index edaf0bbbe..aaea4e050 100644
--- a/src/BlockUniTensor.cpp
+++ b/src/BlockUniTensor.cpp
@@ -1,15 +1,20 @@
 #include "UniTensor.hpp"
+
+#include <map>
+#include <ostream>
+#include <stack>
+#include <vector>
+
+#include <boost/unordered_map.hpp>
+
 #include "Accessor.hpp"
-#include "utils/utils.hpp"
-#include "linalg.hpp"
 #include "Generator.hpp"
-#include <vector>
+#include "linalg.hpp"
+#include "utils/utils.hpp"
 #include "utils/vec_clone.hpp"
-#include "utils/vec_print.hpp"
 #include "utils/vec_concatenate.hpp"
-#include <map>
-#include <boost/unordered_map.hpp>
-#include <stack>
+#include "utils/vec_print.hpp"
+
 using namespace std;
 
 #ifdef BACKEND_TORCH
@@ -287,14 +292,13 @@ namespace cytnx {
 
     os << (std::string(" ") * (3 + Lmax + 5)) << std::string("-") * (4 + mL + mR + 5) << endl;
   }
-  void BlockUniTensor::print_block(const cytnx_int64 &idx, const bool &full_info) const {
+  void BlockUniTensor::print_block(std::ostream &os, const cytnx_int64 &idx,
+                                   const bool &full_info) const {
     cytnx_error_msg(
       (idx < 0) || (idx >= this->_blocks.size()),
       "[ERROR] index [%d] out of bound. should be >0 and < number of available blocks %d\n", idx,
       this->_blocks.size());
 
-    std::ostream &os = std::cout;
-
     os << "========================\n";
     if (this->_is_diag) os << " *is_diag: True\n";
     os << "BLOCK [#" << idx << "]\n";
@@ -335,9 +339,7 @@ namespace cytnx {
     }
   }
 
-  void BlockUniTensor::print_blocks(const bool &full_info) const {
-    std::ostream &os = std::cout;
-
+  void BlockUniTensor::print_blocks(std::ostream &os, const bool &full_info) const {
     os << "-------- start of print ---------\n";
     char *buffer = (char *)malloc(sizeof(char) * 10240);
     sprintf(buffer, "Tensor name: %s\n", this->_name.c_str());
@@ -361,7 +363,7 @@ namespace cytnx {
 
     // print each blocks with its qnum!
     for (int b = 0; b < this->_blocks.size(); b++) {
-      this->print_block(b, full_info);
+      this->print_block(os, b, full_info);
     }
 
     /*
@@ -388,26 +390,26 @@ namespace cytnx {
     free(buffer);
   }
 
-  void BlockUniTensor::print_diagram(const bool &bond_info) const {
+  void BlockUniTensor::print_diagram(std::ostream &os, const bool &bond_info) const {
     char *buffer = (char *)malloc(10240 * sizeof(char));
     unsigned int BUFFsize = 100;
 
     sprintf(buffer, "-----------------------%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Name : %s\n", this->_name.c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Rank : %d\n", this->_labels.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     // sprintf(buffer, "block_form  : true%s", "\n");
-    // std::cout << std::string(buffer);
+    // os << std::string(buffer);
     sprintf(buffer, "contiguous  : %s\n", this->is_contiguous() ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "valid blocks: %d\n", this->_blocks.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "is diag     : %s\n", this->is_diag() ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "on device   : %s\n", this->device_str().c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     cytnx_uint64 Nin = this->_rowrank;
     cytnx_uint64 Nout = this->_labels.size() - this->_rowrank;
@@ -443,12 +445,12 @@ namespace cytnx {
 
     std::string tmpss;
     sprintf(buffer, "%s row %s col %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s    -%s-    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     for (cytnx_uint64 i = 0; i < vl; i++) {
       sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
 
       if (i < Nin) {
         if (this->_bonds[i].type() == bondType::BD_KET)
@@ -492,21 +494,21 @@ namespace cytnx {
         sprintf(rlbl, "%s", tmpss.c_str());
       }
       sprintf(buffer, "   %s| %s     %s |%s\n", l, llbl, rlbl, r);
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
     }
     sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s    -%s-    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     if (bond_info) {
       for (cytnx_uint64 i = 0; i < this->_bonds.size(); i++) {
         // sprintf(buffer, "lbl:%d ", this->_labels[i]);
         sprintf(buffer, "lbl:%s ", this->_labels[i].c_str());
-        std::cout << std::string(buffer);
-        std::cout << this->_bonds[i] << std::endl;
+        os << std::string(buffer);
+        os << this->_bonds[i] << std::endl;
       }
     }
 
diff --git a/src/Bond.cpp b/src/Bond.cpp
index 989a2b681..42ea7a781 100644
--- a/src/Bond.cpp
+++ b/src/Bond.cpp
@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <filesystem>
+#include <ostream>
 
 #include "utils/utils.hpp"
 
diff --git a/src/DenseUniTensor.cpp b/src/DenseUniTensor.cpp
index 3dac3dba3..6a06c52b1 100644
--- a/src/DenseUniTensor.cpp
+++ b/src/DenseUniTensor.cpp
@@ -1,11 +1,14 @@
 #include "UniTensor.hpp"
-#include "utils/utils.hpp"
 
-#include "Generator.hpp"
-#include "linalg.hpp"
 #include <algorithm>
+#include <ostream>
 #include <utility>
 #include <vector>
+
+#include "Generator.hpp"
+#include "linalg.hpp"
+#include "utils/utils.hpp"
+
 typedef cytnx::Accessor ac;
 
 using namespace std;
@@ -416,8 +419,8 @@ namespace cytnx {
     }
   }
 
-  void DenseUniTensor::print_block(const cytnx_int64 &idx, const bool &full_info) const {
-    std::ostream &os = std::cout;
+  void DenseUniTensor::print_block(std::ostream &os, const cytnx_int64 &idx,
+                                   const bool &full_info) const {
     os << "-------- start of print ---------\n";
     char *buffer = (char *)malloc(sizeof(char) * 10240);
     sprintf(buffer, "Tensor name: %s\n", this->_name.c_str());
@@ -443,26 +446,26 @@ namespace cytnx {
     }
     free(buffer);
   }
-  void DenseUniTensor::print_blocks(const bool &full_info) const {
-    this->print_block(0, full_info);
+  void DenseUniTensor::print_blocks(std::ostream &os, const bool &full_info) const {
+    this->print_block(os, 0, full_info);
   }
 
-  void DenseUniTensor::print_diagram(const bool &bond_info) const {
+  void DenseUniTensor::print_diagram(std::ostream &os, const bool &bond_info) const {
     char *buffer = (char *)malloc(10240 * sizeof(char));
     unsigned int BUFFsize = 100;
 
     sprintf(buffer, "-----------------------%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Name : %s\n", this->_name.c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Rank : %ld\n", this->_labels.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "block_form  : False%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "is_diag     : %s\n", this->_is_diag ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "on device   : %s\n", this->device_str().c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     cytnx_uint64 Nin = this->_rowrank;
     cytnx_uint64 Nout = this->_labels.size() - this->_rowrank;
@@ -499,15 +502,15 @@ namespace cytnx {
 
     if (this->is_tag()) {
       sprintf(buffer, "braket_form : %s\n", this->_is_braket_form ? "True" : "False");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
 
       sprintf(buffer, "%s row %s col %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
       sprintf(buffer, "%s    -%s-    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
       for (cytnx_uint64 i = 0; i < vl; i++) {
         sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-        std::cout << std::string(buffer);
+        os << std::string(buffer);
 
         if (i < Nin) {
           if (this->_bonds[i].type() == bondType::BD_KET)
@@ -553,25 +556,25 @@ namespace cytnx {
           sprintf(rlbl, "%s", tmpss.c_str());
         }
         sprintf(buffer, "   %s| %s     %s |%s\n", l, llbl, rlbl, r);
-        std::cout << std::string(buffer);
+        os << std::string(buffer);
       }
       sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
       sprintf(buffer, "%s    -%s-    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
       // sprintf(buffer, "%s", "\n");
-      // std::cout << std::string(buffer);
+      // os << std::string(buffer);
 
     } else {
       sprintf(buffer, "%s     %s     %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
       for (cytnx_uint64 i = 0; i < vl; i++) {
         if (i == 0) {
           sprintf(buffer, "%s    /%s\\    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-          std::cout << std::string(buffer);
+          os << std::string(buffer);
         } else {
           sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-          std::cout << std::string(buffer);
+          os << std::string(buffer);
         }
 
         if (i < Nin) {
@@ -614,20 +617,20 @@ namespace cytnx {
           sprintf(rlbl, "%s", tmpss.c_str());
         }
         sprintf(buffer, "   %s| %s     %s |%s\n", l, llbl, rlbl, r);
-        std::cout << std::string(buffer);
+        os << std::string(buffer);
       }
       sprintf(buffer, "%s    \\%s/    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
       sprintf(buffer, "%s     %s     %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
     }
 
     if (bond_info) {
       for (cytnx_uint64 i = 0; i < this->_bonds.size(); i++) {
         // sprintf(buffer, "lbl:%ld ", this->_labels[i]);
         sprintf(buffer, "lbl:%s ", this->_labels[i].c_str());
-        std::cout << std::string(buffer);
-        std::cout << this->_bonds[i] << std::endl;
+        os << std::string(buffer);
+        os << this->_bonds[i] << std::endl;
       }
     }
 
diff --git a/src/Device.cpp b/src/Device.cpp
index 8215e131e..304251170 100644
--- a/src/Device.cpp
+++ b/src/Device.cpp
@@ -1,5 +1,6 @@
 #include "Device.hpp"
 
+#include <ostream>
 #include <thread>
 
 #include "cytnx_error.hpp"
@@ -59,42 +60,42 @@ namespace cytnx {
       return string("");
     }
   }
-  void Device_class::Print_Property() {
+  void Device_class::Print_Property(std::ostream& os) {
     char* buffer = (char*)malloc(sizeof(char) * 256);
 #ifdef UNI_GPU
-    cout << "=== CUDA support ===" << endl;
-    cout << ": Peer PCIE Access:" << endl;
-    cout << "   ";
+    os << "=== CUDA support ===" << std::endl;
+    os << ": Peer PCIE Access:" << std::endl;
+    os << "   ";
     for (int i = 0; i < this->Ngpus; i++) {
       sprintf(buffer, " %2d", i);
-      cout << string(buffer);
+      os << string(buffer);
     }
-    cout << endl;
+    os << std::endl;
 
-    cout << "   ";
+    os << "   ";
     for (int i = 0; i < this->Ngpus; i++) {
       sprintf(buffer, "%s", "---");
-      cout << string(buffer);
+      os << string(buffer);
     }
-    cout << endl;
+    os << std::endl;
 
     for (int i = 0; i < this->Ngpus; i++) {
       sprintf(buffer, "%2d|", i);
-      cout << string(buffer);
+      os << string(buffer);
       for (int j = 0; j < this->Ngpus; j++) {
         if (j == i) {
           sprintf(buffer, "%s", "  x");
-          cout << string(buffer);
+          os << string(buffer);
         } else {
           sprintf(buffer, "  %d", int(CanAccessPeer[i][j]));
-          cout << string(buffer);
+          os << string(buffer);
         }
       }
-      cout << endl;
+      os << std::endl;
     }
-    cout << "--------------------" << endl;
+    os << "--------------------" << std::endl;
 #else
-    cout << "=== No CUDA support ===" << endl;
+    os << "=== No CUDA support ===" << std::endl;
 #endif
 
     free(buffer);
diff --git a/src/FermionNetwork.cpp b/src/FermionNetwork.cpp
index fe919a8d8..57dbc67b4 100644
--- a/src/FermionNetwork.cpp
+++ b/src/FermionNetwork.cpp
@@ -2,7 +2,6 @@
 #include "Network.hpp"
 #include <stack>
 #include <algorithm>
-#include <iostream>
 
 using namespace std;
 
diff --git a/src/Generator.cpp b/src/Generator.cpp
index 7b0616643..4d9144542 100644
--- a/src/Generator.cpp
+++ b/src/Generator.cpp
@@ -1,10 +1,9 @@
 #include "Generator.hpp"
 
-#include "utils/utils.hpp"
+#include <cfloat>
 
 #include "linalg.hpp"
-#include <cfloat>
-#include <iostream>
+#include "utils/utils.hpp"
 
 #ifdef BACKEND_TORCH
 #else
diff --git a/src/Gncon_base.cpp b/src/Gncon_base.cpp
index 334783a0b..85de88b2d 100644
--- a/src/Gncon_base.cpp
+++ b/src/Gncon_base.cpp
@@ -1,6 +1,9 @@
-#include <typeinfo>
+
 #include "Gncon.hpp"
 
+#include <ostream>
+#include <typeinfo>
+
 #include "linalg.hpp"
 
 using namespace std;
@@ -48,7 +51,7 @@ namespace cytnx {
     cytnx_error_msg(true, "[ERROR][Gncon][clone] call from uninitialize Gncon. %s", "\n");
     return nullptr;
   }
-  void Gncon_base::PrintNet(std::ostream &os) {
+  void Gncon_base::PrintNet(std::ostream &os) const {
     os << "        [Empty Gncon]" << endl;
     os << "--- Please Load Gncon file ---\n";
     os << endl;
diff --git a/src/LinOp.cpp b/src/LinOp.cpp
index 9f58e87f4..3747db407 100644
--- a/src/LinOp.cpp
+++ b/src/LinOp.cpp
@@ -1,4 +1,7 @@
 #include "LinOp.hpp"
+
+#include <ostream>
+
 #include "Tensor.hpp"
 #include "utils/vec_print.hpp"
 
@@ -7,16 +10,16 @@
 
 namespace cytnx {
 
-  void LinOp::_print() {
+  void LinOp::_print(std::ostream &os) {
     if (this->_type == "mv_elem") {
-      std::cout << "elements :" << std::endl;
+      os << "elements :" << std::endl;
       for (auto it = this->_elems.begin(); it != this->_elems.end(); it++) {
-        std::cout << "row:" << it->first << std::endl;
+        os << "row:" << it->first << std::endl;
         // Tensor e_i = it->second.second;
-        std::cout << "col idx:" << std::endl;
-        std::cout << it->second.first << std::endl;
-        std::cout << "elem:" << std::endl;
-        std::cout << it->second.second << std::endl;
+        os << "col idx:" << std::endl;
+        os << it->second.first << std::endl;
+        os << "elem:" << std::endl;
+        os << it->second.second << std::endl;
       }
     }
   }
diff --git a/src/Network.cpp b/src/Network.cpp
index 4d852a8db..e0ac6f172 100644
--- a/src/Network.cpp
+++ b/src/Network.cpp
@@ -1,6 +1,8 @@
-#include <typeinfo>
 #include "Network.hpp"
 
+#include <ostream>
+#include <typeinfo>
+
 #include "linalg.hpp"
 
 using namespace std;
diff --git a/src/Network_base.cpp b/src/Network_base.cpp
index 6877fc8c7..1a54d6422 100644
--- a/src/Network_base.cpp
+++ b/src/Network_base.cpp
@@ -1,6 +1,8 @@
-#include <typeinfo>
 #include "Network.hpp"
 
+#include <ostream>
+#include <typeinfo>
+
 #include "linalg.hpp"
 
 using namespace std;
@@ -81,7 +83,7 @@ namespace cytnx {
     cytnx_error_msg(true, "[ERROR][Network][clone] call from uninitialized network. %s", "\n");
     return nullptr;
   }
-  void Network_base::PrintNet(std::ostream &os) {
+  void Network_base::PrintNet(std::ostream &os) const {
     os << "        [Empty Network]" << endl;
     os << "--- Please Load Network file ---\n";
     os << endl;
diff --git a/src/Physics.cpp b/src/Physics.cpp
index d0fd6687c..553e57baf 100644
--- a/src/Physics.cpp
+++ b/src/Physics.cpp
@@ -5,7 +5,6 @@
 
 #include "linalg.hpp"
 #include <cfloat>
-#include <iostream>
 #include <cmath>
 using namespace std;
 
diff --git a/src/RegularGncon.cpp b/src/RegularGncon.cpp
index 2dc4bd2aa..394666dc1 100644
--- a/src/RegularGncon.cpp
+++ b/src/RegularGncon.cpp
@@ -1,10 +1,11 @@
-#include <typeinfo>
 #include "Gncon.hpp"
 
-#include "search_tree.hpp"
-#include <stack>
 #include <algorithm>
 #include <iostream>
+#include <stack>
+#include <typeinfo>
+
+#include "search_tree.hpp"
 
 using namespace std;
 
@@ -134,17 +135,17 @@ namespace cytnx {
 
   /// This is debug function
   void print_gn(std::vector<vector<pair<string, string>>> &table, vector<string> &names,
-                map<string, cytnx_uint64> &name2pos) {
-    std::cout << "### table  ###" << std::endl;
+                map<string, cytnx_uint64> &name2pos, std::ostream &os = std::cout) {
+    os << "### table  ###" << std::endl;
     for (int i = 0; i < table.size(); i++) {
       for (int j = 0; j < table[i].size(); j++) {
-        std::cout << "original lbl, replaced lbl =  " << table[i][j].first << " ,"
-                  << table[i][j].second << std::endl;
+        os << "original lbl, replaced lbl =  " << table[i][j].first << " ," << table[i][j].second
+           << std::endl;
       }
     }
-    std::cout << "### tensor names  ###" << std::endl;
+    os << "### tensor names  ###" << std::endl;
     for (int i = 0; i < names.size(); i++) {
-      std::cout << names[i] << " pos : " << name2pos[names[i]] << std::endl;
+      os << names[i] << " pos : " << name2pos[names[i]] << std::endl;
     }
   }
 
@@ -385,9 +386,10 @@ namespace cytnx {
         TN_names.erase(it);
       }
       if (TN_names.size() != 0) {
-        cerr << "[ERROR] Following TNs appeared in ORDER line, but is not defined." << endl;
+        std::cerr << "[ERROR] The following TNs appeared in ORDER line, but are not defined."
+                  << std::endl;
         for (int i = 0; i < TN_names.size(); i++) {
-          cerr << "        " << TN_names[i] << endl;
+          std::cerr << "        " << TN_names[i] << std::endl;
         }
         cytnx_error_msg(true, "%s", "\n");
       }
@@ -544,7 +546,7 @@ namespace cytnx {
     this->PutUniTensor(idx, utensor);
   }
 
-  void RegularGncon::PrintNet(std::ostream &os) {
+  void RegularGncon::PrintNet(std::ostream &os) const {
     string status;
     os << "==== Gncon ====" << endl;
     if (this->tensors.size() == 0) {
diff --git a/src/RegularNetwork.cpp b/src/RegularNetwork.cpp
index 5c7600993..ac05c8f38 100644
--- a/src/RegularNetwork.cpp
+++ b/src/RegularNetwork.cpp
@@ -725,7 +725,7 @@ namespace cytnx {
     this->PutUniTensor(idx, utensor);
   }
 
-  void RegularNetwork::PrintNet(ostream &os) {
+  void RegularNetwork::PrintNet(ostream &os) const {
     string status;
     os << "==== Network ====" << endl;
     if (this->tensors.size() == 0) {
diff --git a/src/SparseUniTensor.cpp b/src/SparseUniTensor.cpp
index dcd78c812..e1ff237e9 100644
--- a/src/SparseUniTensor.cpp
+++ b/src/SparseUniTensor.cpp
@@ -1,13 +1,17 @@
 #include "UniTensor.hpp"
+
+#include <map>
+#include <ostream>
+#include <vector>
+
 #include "Accessor.hpp"
+#include "Generator.hpp"
+#include "UniTensor.hpp"
+#include "linalg.hpp"
 #include "utils/utils.hpp"
 #include "utils/utils_internal_interface.hpp"
-#include "linalg.hpp"
-#include "Generator.hpp"
-#include <vector>
-#include "utils/vec_print.hpp"
 #include "utils/vec_concatenate.hpp"
-#include <map>
+#include "utils/vec_print.hpp"
 
 using namespace std;
 namespace cytnx {
@@ -644,12 +648,12 @@ namespace cytnx {
 
   std::vector<Symmetry> SparseUniTensor::syms() const { return this->_bonds[0].syms(); }
 
-  void SparseUniTensor::print_block(const cytnx_int64 &idx, const bool &full_info) const {
+  void SparseUniTensor::print_block(std::ostream &os, const cytnx_int64 &idx,
+                                    const bool &full_info) const {
     cytnx_error_msg(true, "[ERROR] SpareUniTensor does not support individual printing.%s", "\n");
   }
 
-  void SparseUniTensor::print_blocks(const bool &full_info) const {
-    std::ostream &os = std::cout;
+  void SparseUniTensor::print_blocks(std::ostream &os, const bool &full_info) const {
     os << "-------- start of print ---------\n";
     char *buffer = (char *)malloc(sizeof(char) * 1024);
     sprintf(buffer, "Tensor name: %s\n", this->_name.c_str());
@@ -692,26 +696,26 @@ namespace cytnx {
     free(buffer);
   }
 
-  void SparseUniTensor::print_diagram(const bool &bond_info) const {
+  void SparseUniTensor::print_diagram(std::ostream &os, const bool &bond_info) const {
     char *buffer = (char *)malloc(1024 * sizeof(char));
     unsigned int BUFFsize = 100;
 
     sprintf(buffer, "-----------------------%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Name : %s\n", this->_name.c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "tensor Rank : %d\n", this->_labels.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "block_form  : true%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "contiguous  : %s\n", this->is_contiguous() ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "valid bocks : %d\n", this->_blocks.size());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "is diag   : %s\n", this->is_diag() ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "on device   : %s\n", this->device_str().c_str());
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     cytnx_uint64 Nin = this->_rowrank;
     cytnx_uint64 Nout = this->_labels.size() - this->_rowrank;
     cytnx_uint64 vl;
@@ -745,16 +749,16 @@ namespace cytnx {
     string M_dashes = string("-") * (1 + Space_Ldim_max + 5 + Space_Rdim_max + 1);
 
     sprintf(buffer, "braket_form : %s\n", this->_is_braket_form ? "True" : "False");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     std::string tmpss;
     sprintf(buffer, "%s row %s col %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s    -%s-    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     for (cytnx_uint64 i = 0; i < vl; i++) {
       sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
 
       if (i < Nin) {
         if (this->_bonds[i].type() == bondType::BD_KET)
@@ -798,21 +802,21 @@ namespace cytnx {
         sprintf(rlbl, "%s", tmpss.c_str());
       }
       sprintf(buffer, "   %s| %s     %s |%s\n", l, llbl, rlbl, r);
-      std::cout << std::string(buffer);
+      os << std::string(buffer);
     }
     sprintf(buffer, "%s    |%s|    %s", LallSpace.c_str(), MallSpace.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s    -%s-    %s", LallSpace.c_str(), M_dashes.c_str(), "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
     sprintf(buffer, "%s", "\n");
-    std::cout << std::string(buffer);
+    os << std::string(buffer);
 
     if (bond_info) {
       for (cytnx_uint64 i = 0; i < this->_bonds.size(); i++) {
         // sprintf(buffer, "lbl:%d ", this->_labels[i]);
         sprintf(buffer, "lbl:%s ", this->_labels[i].c_str());
-        std::cout << std::string(buffer);
-        std::cout << this->_bonds[i] << std::endl;
+        os << std::string(buffer);
+        os << this->_bonds[i] << std::endl;
       }
     }
 
diff --git a/src/Symmetry.cpp b/src/Symmetry.cpp
index 081ff1def..bfa3c0775 100644
--- a/src/Symmetry.cpp
+++ b/src/Symmetry.cpp
@@ -2,7 +2,7 @@
 
 #include <filesystem>
 #include <fstream>
-#include <iostream>
+#include <ostream>
 #include <string>
 #include <vector>
 
@@ -73,7 +73,7 @@ namespace cytnx {
     return EVEN;
   }
 
-  void cytnx::Symmetry_base::print_info() const {
+  void cytnx::Symmetry_base::print_info(std::ostream &os) const {
     cytnx_error_msg(true, "%s", "[ERROR][Internal] should not call Symmerty base!");
   }
 
@@ -102,13 +102,13 @@ namespace cytnx {
   }
   void cytnx::U1Symmetry::reverse_rule_(cytnx_int64 &out, const cytnx_int64 &in) { out = in * -1; }
 
-  void cytnx::U1Symmetry::print_info() const {
-    cout << "--------------------\n";
-    cout << "[Symmetry]" << endl;
-    cout << "type : Abelian, U1" << endl;
-    cout << "combine rule : Q1 + Q2" << endl;
-    cout << "reverse rule : Q*(-1) " << endl;
-    cout << "--------------------\n";
+  void cytnx::U1Symmetry::print_info(std::ostream &os) const {
+    os << "--------------------\n";
+    os << "[Symmetry]" << std::endl;
+    os << "type : Abelian, U1" << std::endl;
+    os << "combine rule : Q1 + Q2" << std::endl;
+    os << "reverse rule : Q*(-1) " << std::endl;
+    os << "--------------------\n";
   }
 
   ///========================
@@ -152,13 +152,13 @@ namespace cytnx {
     out = (this->n - in) % this->n;
   }
 
-  void cytnx::ZnSymmetry::print_info() const {
-    cout << "--------------------\n";
-    cout << "[Symmetry]" << endl;
-    cout << "type : Abelian, Z(" << this->n << ")" << endl;
-    cout << "combine rule : (Q1 + Q2)\%" << this->n << endl;
-    cout << "reverse rule : Q*(-1) " << endl;
-    cout << "--------------------\n";
+  void cytnx::ZnSymmetry::print_info(std::ostream &os) const {
+    os << "--------------------\n";
+    os << "[Symmetry]" << std::endl;
+    os << "type : Abelian, Z(" << this->n << ")" << std::endl;
+    os << "combine rule : (Q1 + Q2)\%" << this->n << std::endl;
+    os << "reverse rule : Q*(-1) " << std::endl;
+    os << "--------------------\n";
   }
 
   ///========================
@@ -206,13 +206,13 @@ namespace cytnx {
                       "\n");
   }
 
-  void cytnx::FermionParitySymmetry::print_info() const {
-    cout << "--------------------\n";
-    cout << "[Symmetry]" << endl;
-    cout << "type : fermionic, FermionParity" << endl;
-    cout << "combine rule : (Q1 + Q2)\2" << endl;
-    cout << "reverse rule : Q*(-1) " << endl;
-    cout << "--------------------\n";
+  void cytnx::FermionParitySymmetry::print_info(std::ostream &os) const {
+    os << "--------------------\n";
+    os << "[Symmetry]" << std::endl;
+    os << "type : fermionic, FermionParity" << std::endl;
+    os << "combine rule : (Q1 + Q2)\2" << std::endl;
+    os << "reverse rule : Q*(-1) " << std::endl;
+    os << "--------------------\n";
   }
 
   ///=========================
@@ -248,13 +248,13 @@ namespace cytnx {
       return ODD;
   }
 
-  void cytnx::FermionNumberSymmetry::print_info() const {
-    cout << "--------------------\n";
-    cout << "[Symmetry]" << endl;
-    cout << "type : fermionic, FermionNumber" << endl;
-    cout << "combine rule : Q1 + Q2" << endl;
-    cout << "reverse rule : Q*(-1) " << endl;
-    cout << "--------------------\n";
+  void cytnx::FermionNumberSymmetry::print_info(std::ostream &os) const {
+    os << "--------------------\n";
+    os << "[Symmetry]" << std::endl;
+    os << "type : fermionic, FermionNumber" << std::endl;
+    os << "combine rule : Q1 + Q2" << std::endl;
+    os << "reverse rule : Q*(-1) " << std::endl;
+    os << "--------------------\n";
   }
 
   //==================================================
@@ -315,8 +315,8 @@ namespace cytnx {
     this->Init(stype_in, n_in);
   }
 
-  ostream &operator<<(ostream &os, const Symmetry &in) {
-    in.print_info();
+  std::ostream &operator<<(std::ostream &os, const Symmetry &in) {
+    in.print_info(os);
     return os;
   }
 };  // namespace cytnx
diff --git a/src/Tensor.cpp b/src/Tensor.cpp
index 7d5d51544..d9fda4150 100644
--- a/src/Tensor.cpp
+++ b/src/Tensor.cpp
@@ -1,6 +1,7 @@
 #include "Tensor.hpp"
 
 #include <filesystem>
+#include <ostream>
 #include <typeinfo>
 
 #include "linalg.hpp"
diff --git a/src/UniTensor_base.cpp b/src/UniTensor_base.cpp
index 4e625b64b..140d5975f 100644
--- a/src/UniTensor_base.cpp
+++ b/src/UniTensor_base.cpp
@@ -1,8 +1,11 @@
 #include "UniTensor.hpp"
+
+#include <ostream>
+#include <string>
+
+#include "Network.hpp"
 #include "Tensor.hpp"
 #include "linalg.hpp"
-#include "Network.hpp"
-#include <string>
 
 #ifdef BACKEND_TORCH
 #else
@@ -236,15 +239,16 @@ namespace cytnx {
     return out;
   }
 
-  void UniTensor_base::print_diagram(const bool &bond_info) const {
+  void UniTensor_base::print_diagram(std::ostream &os, const bool &bond_info) const {
     cytnx_error_msg(
       true, "[ERROR] fatal internal, cannot call on an un-initialized UniTensor_base%s", "\n");
   }
-  void UniTensor_base::print_blocks(const bool &full_info) const {
+  void UniTensor_base::print_blocks(std::ostream &os, const bool &full_info) const {
     cytnx_error_msg(
       true, "[ERROR] fatal internal, cannot call on an un-initialized UniTensor_base%s", "\n");
   }
-  void UniTensor_base::print_block(const cytnx_int64 &idx, const bool &full_info) const {
+  void UniTensor_base::print_block(std::ostream &os, const cytnx_int64 &idx,
+                                   const bool &full_info) const {
     cytnx_error_msg(
       true, "[ERROR] fatal internal, cannot call on an un-initialized UniTensor_base%s", "\n");
   }
@@ -711,7 +715,7 @@ namespace cytnx {
 
   //-----------------------------------------
   std::ostream &operator<<(std::ostream &os, const UniTensor &in) {
-    in.print_blocks();
+    in.print_blocks(os);
     return os;
   }
 
diff --git a/src/backend/Scalar.cpp b/src/backend/Scalar.cpp
index 0f0303bcd..d9f180eaa 100644
--- a/src/backend/Scalar.cpp
+++ b/src/backend/Scalar.cpp
@@ -1,4 +1,7 @@
 #include "backend/Scalar.hpp"
+
+#include <ostream>
+
 #include "backend/Storage.hpp"
 
 namespace cytnx {
diff --git a/src/backend/Storage.cpp b/src/backend/Storage.cpp
index 8415d62a3..c7305dadc 100644
--- a/src/backend/Storage.cpp
+++ b/src/backend/Storage.cpp
@@ -1,7 +1,7 @@
 #include "backend/Storage.hpp"
 
 #include <filesystem>
-#include <iostream>
+#include <ostream>
 
 using namespace std;
 
@@ -10,7 +10,7 @@ namespace cytnx {
   Storage_init_interface __SII;
 
   std::ostream &operator<<(std::ostream &os, const Storage &in) {
-    in.print();
+    in.print(os);
     return os;
   }
 
diff --git a/src/backend/StorageImplementation.cpp b/src/backend/StorageImplementation.cpp
index 4d7bbf83b..4e525d1fe 100644
--- a/src/backend/StorageImplementation.cpp
+++ b/src/backend/StorageImplementation.cpp
@@ -4,7 +4,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <iomanip>
-#include <iostream>
+#include <ostream>
 #include <type_traits>
 #include <vector>
 
@@ -447,13 +447,13 @@ namespace cytnx {
   }
 
   template <typename DType>
-  void StorageImplementation<DType>::print_elems() {
+  void StorageImplementation<DType>::print_elems(std::ostream &os) const {
     DType *elem_ptr_ = reinterpret_cast<DType *>(this->start_);
-    cout << "[ ";
+    os << "[ ";
     for (unsigned long long cnt = 0; cnt < this->size_; cnt++) {
-      PrintValueAndSpace(std::cout, elem_ptr_[cnt]);
+      PrintValueAndSpace(os, elem_ptr_[cnt]);
     }
-    std::cout << "]" << std::endl;
+    os << "]" << std::endl;
   }
 
   template <typename DType>
diff --git a/src/backend/Storage_base.cpp b/src/backend/Storage_base.cpp
index 11bdfec5a..a6c8fb5b8 100644
--- a/src/backend/Storage_base.cpp
+++ b/src/backend/Storage_base.cpp
@@ -1,4 +1,7 @@
 #include "backend/Storage.hpp"
+
+#include <ostream>
+
 #include "utils_internal_interface.hpp"
 #include "utils/vec_print.hpp"
 
@@ -163,17 +166,17 @@ namespace cytnx {
     cytnx_error_msg(true, "%s", "[ERROR] call PrintElem_byShape directly on Void Storage.");
   }
 
-  void Storage_base::print_info() {
-    cout << "dtype : " << this->dtype_str() << endl;
-    cout << "device: " << Device.getname(this->device()) << endl;
-    cout << "size  : " << this->size() << endl;
+  void Storage_base::print_info(std::ostream &os) const {
+    os << "dtype : " << this->dtype_str() << std::endl;
+    os << "device: " << Device.getname(this->device()) << std::endl;
+    os << "size  : " << this->size() << std::endl;
   }
-  void Storage_base::print_elems() {
+  void Storage_base::print_elems(std::ostream &os) const {
     cytnx_error_msg(true, "%s", "[ERROR] call print_elems directly on Void Storage.");
   }
-  void Storage_base::print() {
-    this->print_info();
-    this->print_elems();
+  void Storage_base::print(std::ostream &os) const {
+    this->print_info(os);
+    this->print_elems(os);
   }
 
   // shadow new:
diff --git a/src/backend/algo_internal_cpu/Concate_internal.cpp b/src/backend/algo_internal_cpu/Concate_internal.cpp
index b5919c5de..b7c11a465 100644
--- a/src/backend/algo_internal_cpu/Concate_internal.cpp
+++ b/src/backend/algo_internal_cpu/Concate_internal.cpp
@@ -1,5 +1,4 @@
 #include "Concate_internal.hpp"
-#include <iostream>
 #include <algorithm>
 
 namespace cytnx {
diff --git a/src/backend/algo_internal_cpu/Concate_internal.hpp b/src/backend/algo_internal_cpu/Concate_internal.hpp
index 9f585e74c..1b7799d07 100644
--- a/src/backend/algo_internal_cpu/Concate_internal.hpp
+++ b/src/backend/algo_internal_cpu/Concate_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_ALGO_INTERNAL_CPU_CONCATE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/algo_internal_cpu/Split_internal.cpp b/src/backend/algo_internal_cpu/Split_internal.cpp
index 18e2bfc6c..72aaac767 100644
--- a/src/backend/algo_internal_cpu/Split_internal.cpp
+++ b/src/backend/algo_internal_cpu/Split_internal.cpp
@@ -1,5 +1,4 @@
 #include "Split_internal.hpp"
-#include <iostream>
 #include <algorithm>
 
 namespace cytnx {
diff --git a/src/backend/algo_internal_cpu/Split_internal.hpp b/src/backend/algo_internal_cpu/Split_internal.hpp
index 02247d2b8..c05cc24d9 100644
--- a/src/backend/algo_internal_cpu/Split_internal.hpp
+++ b/src/backend/algo_internal_cpu/Split_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_ALGO_INTERNAL_CPU_SPLIT_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/algo_internal_gpu/cuConcate_internal.hpp b/src/backend/algo_internal_gpu/cuConcate_internal.hpp
index c1b3a83a0..3ec681786 100644
--- a/src/backend/algo_internal_gpu/cuConcate_internal.hpp
+++ b/src/backend/algo_internal_gpu/cuConcate_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_ALGO_INTERNAL_GPU_CUCONCATE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/algo_internal_gpu/cuSplit_internal.hpp b/src/backend/algo_internal_gpu/cuSplit_internal.hpp
index dd2f8e62c..f8389618c 100644
--- a/src/backend/algo_internal_gpu/cuSplit_internal.hpp
+++ b/src/backend/algo_internal_gpu/cuSplit_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_ALGO_INTERNAL_GPU_CUSPLIT_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Abs_internal.hpp b/src/backend/linalg_internal_cpu/Abs_internal.hpp
index f76708819..5dbe7e049 100644
--- a/src/backend/linalg_internal_cpu/Abs_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Abs_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_ABS_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Axpy_internal.hpp b/src/backend/linalg_internal_cpu/Axpy_internal.hpp
index 2d404ff3e..fa7e5777a 100644
--- a/src/backend/linalg_internal_cpu/Axpy_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Axpy_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_AXPY_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Batch_matmul_internal.hpp b/src/backend/linalg_internal_cpu/Batch_matmul_internal.hpp
index 4fc225169..99d4b448d 100644
--- a/src/backend/linalg_internal_cpu/Batch_matmul_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Batch_matmul_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_BATCH_MATMUL_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Conj_inplace_internal.hpp b/src/backend/linalg_internal_cpu/Conj_inplace_internal.hpp
index 371c6e33e..bc0c86d89 100644
--- a/src/backend/linalg_internal_cpu/Conj_inplace_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Conj_inplace_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_CONJ_INPLACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Det_internal.hpp b/src/backend/linalg_internal_cpu/Det_internal.hpp
index d56435f5e..d713cf3ed 100644
--- a/src/backend/linalg_internal_cpu/Det_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Det_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_DET_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Diag_internal.hpp b/src/backend/linalg_internal_cpu/Diag_internal.hpp
index c227840c6..2be4c9fee 100644
--- a/src/backend/linalg_internal_cpu/Diag_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Diag_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_DIAG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Eig_internal.hpp b/src/backend/linalg_internal_cpu/Eig_internal.hpp
index 1feb84fc4..12ca810a1 100644
--- a/src/backend/linalg_internal_cpu/Eig_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Eig_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_EIG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Eigh_internal.hpp b/src/backend/linalg_internal_cpu/Eigh_internal.hpp
index ec421f44d..f3be91a07 100644
--- a/src/backend/linalg_internal_cpu/Eigh_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Eigh_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_EIGH_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Exp_internal.hpp b/src/backend/linalg_internal_cpu/Exp_internal.hpp
index 0ec02eb07..b65c6c041 100644
--- a/src/backend/linalg_internal_cpu/Exp_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Exp_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_EXP_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Gemm_Batch_internal.hpp b/src/backend/linalg_internal_cpu/Gemm_Batch_internal.hpp
index 1c99281bf..1fe51ea45 100644
--- a/src/backend/linalg_internal_cpu/Gemm_Batch_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Gemm_Batch_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_GEMM_BATCH_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Gemm_internal.hpp b/src/backend/linalg_internal_cpu/Gemm_internal.hpp
index 57690bc8c..397fba0f4 100644
--- a/src/backend/linalg_internal_cpu/Gemm_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Gemm_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_GEMM_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Ger_internal.hpp b/src/backend/linalg_internal_cpu/Ger_internal.hpp
index 258153e44..3ccfe5411 100644
--- a/src/backend/linalg_internal_cpu/Ger_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Ger_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_GER_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Gesvd_internal.hpp b/src/backend/linalg_internal_cpu/Gesvd_internal.hpp
index 2c304056c..3f931e6a9 100644
--- a/src/backend/linalg_internal_cpu/Gesvd_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Gesvd_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_GESVD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/InvM_inplace_internal.hpp b/src/backend/linalg_internal_cpu/InvM_inplace_internal.hpp
index 3c697ee35..634aad842 100644
--- a/src/backend/linalg_internal_cpu/InvM_inplace_internal.hpp
+++ b/src/backend/linalg_internal_cpu/InvM_inplace_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_INVM_INPLACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/Inv_inplace_internal.hpp b/src/backend/linalg_internal_cpu/Inv_inplace_internal.hpp
index 97a92b226..fc8f79e69 100644
--- a/src/backend/linalg_internal_cpu/Inv_inplace_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Inv_inplace_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_INV_INPLACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/Kron_internal.hpp b/src/backend/linalg_internal_cpu/Kron_internal.hpp
index 886724d41..071a36548 100644
--- a/src/backend/linalg_internal_cpu/Kron_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Kron_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_KRON_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/Lstsq_internal.hpp b/src/backend/linalg_internal_cpu/Lstsq_internal.hpp
index 5c93d0192..8791b1409 100644
--- a/src/backend/linalg_internal_cpu/Lstsq_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Lstsq_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_LSTSQ_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/Matmul_dg_internal.hpp b/src/backend/linalg_internal_cpu/Matmul_dg_internal.hpp
index 21c8c12ba..1627b9875 100644
--- a/src/backend/linalg_internal_cpu/Matmul_dg_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Matmul_dg_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_MATMUL_DG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/Matmul_internal.hpp b/src/backend/linalg_internal_cpu/Matmul_internal.hpp
index b29486024..07dca30dd 100644
--- a/src/backend/linalg_internal_cpu/Matmul_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Matmul_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_MATMUL_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/Matvec_internal.hpp b/src/backend/linalg_internal_cpu/Matvec_internal.hpp
index 53d95bb80..a81506468 100644
--- a/src/backend/linalg_internal_cpu/Matvec_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Matvec_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_MATVEC_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/MaxMin_internal.hpp b/src/backend/linalg_internal_cpu/MaxMin_internal.hpp
index 7f5a402bc..3a3dcbf2d 100644
--- a/src/backend/linalg_internal_cpu/MaxMin_internal.hpp
+++ b/src/backend/linalg_internal_cpu/MaxMin_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_MAXMIN_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Norm_internal.hpp b/src/backend/linalg_internal_cpu/Norm_internal.hpp
index 2fce1dc35..cac0004b3 100644
--- a/src/backend/linalg_internal_cpu/Norm_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Norm_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_NORM_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Outer_internal.hpp b/src/backend/linalg_internal_cpu/Outer_internal.hpp
index d3d403eef..55910e4c0 100644
--- a/src/backend/linalg_internal_cpu/Outer_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Outer_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_OUTER_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Pow_internal.hpp b/src/backend/linalg_internal_cpu/Pow_internal.hpp
index a4a08de20..3abf6577d 100644
--- a/src/backend/linalg_internal_cpu/Pow_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Pow_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_POW_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/QR_internal.hpp b/src/backend/linalg_internal_cpu/QR_internal.hpp
index cca40d9b3..0be7115cf 100644
--- a/src/backend/linalg_internal_cpu/QR_internal.hpp
+++ b/src/backend/linalg_internal_cpu/QR_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_QR_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Sdd_internal.hpp b/src/backend/linalg_internal_cpu/Sdd_internal.hpp
index 5e5f3151e..7312e3646 100644
--- a/src/backend/linalg_internal_cpu/Sdd_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Sdd_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_SDD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Trace_internal.hpp b/src/backend/linalg_internal_cpu/Trace_internal.hpp
index 68c7720c5..9434a701c 100644
--- a/src/backend/linalg_internal_cpu/Trace_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Trace_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_TRACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Tridiag_internal.cpp b/src/backend/linalg_internal_cpu/Tridiag_internal.cpp
index 588811805..d910748e7 100644
--- a/src/backend/linalg_internal_cpu/Tridiag_internal.cpp
+++ b/src/backend/linalg_internal_cpu/Tridiag_internal.cpp
@@ -1,7 +1,6 @@
 #include "Tridiag_internal.hpp"
 #include "cytnx_error.hpp"
 #include "backend/lapack_wrapper.hpp"
-#include <iostream>
 namespace cytnx {
 
   namespace linalg_internal {
diff --git a/src/backend/linalg_internal_cpu/Tridiag_internal.hpp b/src/backend/linalg_internal_cpu/Tridiag_internal.hpp
index e7494a21d..ea162431f 100644
--- a/src/backend/linalg_internal_cpu/Tridiag_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Tridiag_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_TRIDIAG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/Vectordot_internal.hpp b/src/backend/linalg_internal_cpu/Vectordot_internal.hpp
index 55d2b253b..6a81731ac 100644
--- a/src/backend/linalg_internal_cpu/Vectordot_internal.hpp
+++ b/src/backend/linalg_internal_cpu/Vectordot_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_VECTORDOT_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_cpu/iAdd_internal.cpp b/src/backend/linalg_internal_cpu/iAdd_internal.cpp
index 2a6bc708a..a69a0d1ce 100644
--- a/src/backend/linalg_internal_cpu/iAdd_internal.cpp
+++ b/src/backend/linalg_internal_cpu/iAdd_internal.cpp
@@ -1,7 +1,7 @@
 #include "iAdd_internal.hpp"
+
 #include "backend/utils_internal_interface.hpp"
 #include "utils/utils.hpp"
-#include <iostream>
 
 namespace cytnx {
 
diff --git a/src/backend/linalg_internal_cpu/iAdd_internal.hpp b/src/backend/linalg_internal_cpu/iAdd_internal.hpp
index ce9a30ddd..63d302524 100644
--- a/src/backend/linalg_internal_cpu/iAdd_internal.hpp
+++ b/src/backend/linalg_internal_cpu/iAdd_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_IADD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/iDiv_internal.cpp b/src/backend/linalg_internal_cpu/iDiv_internal.cpp
index b82c95a96..282abbe49 100644
--- a/src/backend/linalg_internal_cpu/iDiv_internal.cpp
+++ b/src/backend/linalg_internal_cpu/iDiv_internal.cpp
@@ -1,7 +1,6 @@
 #include "iDiv_internal.hpp"
 #include "backend/utils_internal_interface.hpp"
 #include "utils/utils.hpp"
-#include <iostream>
 
 namespace cytnx {
 
diff --git a/src/backend/linalg_internal_cpu/iDiv_internal.hpp b/src/backend/linalg_internal_cpu/iDiv_internal.hpp
index e28aaa93c..ec3d3afe7 100644
--- a/src/backend/linalg_internal_cpu/iDiv_internal.hpp
+++ b/src/backend/linalg_internal_cpu/iDiv_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_IDIV_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/iMul_internal.cpp b/src/backend/linalg_internal_cpu/iMul_internal.cpp
index 14e3f4f4c..603d613e2 100644
--- a/src/backend/linalg_internal_cpu/iMul_internal.cpp
+++ b/src/backend/linalg_internal_cpu/iMul_internal.cpp
@@ -1,7 +1,7 @@
 #include "iMul_internal.hpp"
+
 #include "backend/utils_internal_interface.hpp"
 #include "utils/utils.hpp"
-#include <iostream>
 
 namespace cytnx {
 
diff --git a/src/backend/linalg_internal_cpu/iMul_internal.hpp b/src/backend/linalg_internal_cpu/iMul_internal.hpp
index 8f2914728..9a62b5a7f 100644
--- a/src/backend/linalg_internal_cpu/iMul_internal.hpp
+++ b/src/backend/linalg_internal_cpu/iMul_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_IMUL_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_cpu/iSub_internal.cpp b/src/backend/linalg_internal_cpu/iSub_internal.cpp
index 14aa56c8f..262cd480a 100644
--- a/src/backend/linalg_internal_cpu/iSub_internal.cpp
+++ b/src/backend/linalg_internal_cpu/iSub_internal.cpp
@@ -1,7 +1,7 @@
 #include "iSub_internal.hpp"
+
 #include "backend/utils_internal_interface.hpp"
 #include "utils/utils.hpp"
-#include <iostream>
 
 namespace cytnx {
 
diff --git a/src/backend/linalg_internal_cpu/iSub_internal.hpp b/src/backend/linalg_internal_cpu/iSub_internal.hpp
index 8bed0af53..4ae154326 100644
--- a/src/backend/linalg_internal_cpu/iSub_internal.hpp
+++ b/src/backend/linalg_internal_cpu/iSub_internal.hpp
@@ -2,9 +2,9 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_CPU_ISUB_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
+
 #include "backend/Storage.hpp"
 #include "Type.hpp"
 
diff --git a/src/backend/linalg_internal_gpu/cuAbs_internal.hpp b/src/backend/linalg_internal_gpu/cuAbs_internal.hpp
index 566ce96e6..f14552ceb 100644
--- a/src/backend/linalg_internal_gpu/cuAbs_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuAbs_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUABS_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuAdd_internal.hpp b/src/backend/linalg_internal_gpu/cuAdd_internal.hpp
index 860bae9ea..d819ddba3 100644
--- a/src/backend/linalg_internal_gpu/cuAdd_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuAdd_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUADD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuConj_inplace_internal.hpp b/src/backend/linalg_internal_gpu/cuConj_inplace_internal.hpp
index eaeabee02..7f025394b 100644
--- a/src/backend/linalg_internal_gpu/cuConj_inplace_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuConj_inplace_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUCONJ_INPLACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuCpr_internal.hpp b/src/backend/linalg_internal_gpu/cuCpr_internal.hpp
index b9f79aae4..69b48a5cc 100644
--- a/src/backend/linalg_internal_gpu/cuCpr_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuCpr_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUCPR_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuDet_internal.hpp b/src/backend/linalg_internal_gpu/cuDet_internal.hpp
index 2ba94e67f..33cad328e 100644
--- a/src/backend/linalg_internal_gpu/cuDet_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuDet_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUDET_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuDiag_internal.hpp b/src/backend/linalg_internal_gpu/cuDiag_internal.hpp
index 4dd8786fe..de6327184 100644
--- a/src/backend/linalg_internal_gpu/cuDiag_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuDiag_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUDIAG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuDiv_internal.hpp b/src/backend/linalg_internal_gpu/cuDiv_internal.hpp
index a4bee2e61..e9be653ed 100644
--- a/src/backend/linalg_internal_gpu/cuDiv_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuDiv_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUDIV_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuEig_internal.hpp b/src/backend/linalg_internal_gpu/cuEig_internal.hpp
index 9532199c9..162560ce4 100644
--- a/src/backend/linalg_internal_gpu/cuEig_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuEig_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUEIG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuEigh_internal.hpp b/src/backend/linalg_internal_gpu/cuEigh_internal.hpp
index a242427fc..a2371265c 100644
--- a/src/backend/linalg_internal_gpu/cuEigh_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuEigh_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUEIGH_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuExp_internal.hpp b/src/backend/linalg_internal_gpu/cuExp_internal.hpp
index cf0c7677f..3219f9e1d 100644
--- a/src/backend/linalg_internal_gpu/cuExp_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuExp_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUEXP_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuGeSvd_internal.hpp b/src/backend/linalg_internal_gpu/cuGeSvd_internal.hpp
index 453d9436e..b7dc71b9b 100644
--- a/src/backend/linalg_internal_gpu/cuGeSvd_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuGeSvd_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUGESVD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuGemm_Batch_internal.hpp b/src/backend/linalg_internal_gpu/cuGemm_Batch_internal.hpp
index 09bce0c19..fb5e7a869 100644
--- a/src/backend/linalg_internal_gpu/cuGemm_Batch_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuGemm_Batch_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUGEMM_BATCH_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuGemm_internal.hpp b/src/backend/linalg_internal_gpu/cuGemm_internal.hpp
index 5dc66ebb4..c20a10166 100644
--- a/src/backend/linalg_internal_gpu/cuGemm_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuGemm_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUGEMM_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuGer_internal.hpp b/src/backend/linalg_internal_gpu/cuGer_internal.hpp
index e8fe3ef10..112d86a60 100644
--- a/src/backend/linalg_internal_gpu/cuGer_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuGer_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUGER_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuInvM_inplace_internal.hpp b/src/backend/linalg_internal_gpu/cuInvM_inplace_internal.hpp
index 014af906f..430f68e6e 100644
--- a/src/backend/linalg_internal_gpu/cuInvM_inplace_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuInvM_inplace_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUINVM_INPLACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuInv_inplace_internal.hpp b/src/backend/linalg_internal_gpu/cuInv_inplace_internal.hpp
index 01b1820de..dcdcd7d88 100644
--- a/src/backend/linalg_internal_gpu/cuInv_inplace_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuInv_inplace_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUINV_INPLACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuKron_internal.cuh b/src/backend/linalg_internal_gpu/cuKron_internal.cuh
index af7956587..8ee1575b8 100644
--- a/src/backend/linalg_internal_gpu/cuKron_internal.cuh
+++ b/src/backend/linalg_internal_gpu/cuKron_internal.cuh
@@ -9,7 +9,6 @@
 #include "cytnx_error.hpp"
 #include "Type.hpp"
 #include "backend/lapack_wrapper.hpp"
-#include <iostream>
 
 #ifndef __CUDACC__
   #error This file requires a CUDA compiler
diff --git a/src/backend/linalg_internal_gpu/cuMatmul_dg_internal.hpp b/src/backend/linalg_internal_gpu/cuMatmul_dg_internal.hpp
index a3387456d..6616fb4c4 100644
--- a/src/backend/linalg_internal_gpu/cuMatmul_dg_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuMatmul_dg_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUMATMUL_DG_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuMatmul_internal.hpp b/src/backend/linalg_internal_gpu/cuMatmul_internal.hpp
index 104ce35c7..8749c35fa 100644
--- a/src/backend/linalg_internal_gpu/cuMatmul_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuMatmul_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUMATMUL_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuMatvec_internal.hpp b/src/backend/linalg_internal_gpu/cuMatvec_internal.hpp
index 5a442df01..80f233e60 100644
--- a/src/backend/linalg_internal_gpu/cuMatvec_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuMatvec_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUMATVEC_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuMaxMin_internal.hpp b/src/backend/linalg_internal_gpu/cuMaxMin_internal.hpp
index 70db5e60b..4eb164013 100644
--- a/src/backend/linalg_internal_gpu/cuMaxMin_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuMaxMin_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUMAXMIN_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuMod_internal.hpp b/src/backend/linalg_internal_gpu/cuMod_internal.hpp
index 7f5cbe5eb..d7e8999ad 100644
--- a/src/backend/linalg_internal_gpu/cuMod_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuMod_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUMOD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuMul_internal.hpp b/src/backend/linalg_internal_gpu/cuMul_internal.hpp
index 8f393e4d7..0694aecc1 100644
--- a/src/backend/linalg_internal_gpu/cuMul_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuMul_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUMUL_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuNorm_internal.hpp b/src/backend/linalg_internal_gpu/cuNorm_internal.hpp
index 2ee6639a8..e9dcf4027 100644
--- a/src/backend/linalg_internal_gpu/cuNorm_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuNorm_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUNORM_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuOuter_internal.hpp b/src/backend/linalg_internal_gpu/cuOuter_internal.hpp
index d1e1d7bfe..e833727d3 100644
--- a/src/backend/linalg_internal_gpu/cuOuter_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuOuter_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUOUTER_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuPow_internal.hpp b/src/backend/linalg_internal_gpu/cuPow_internal.hpp
index 993e20776..db95fb425 100644
--- a/src/backend/linalg_internal_gpu/cuPow_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuPow_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUPOW_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.cu b/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.cu
index 794521096..a1e63d596 100644
--- a/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.cu
+++ b/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.cu
@@ -1,6 +1,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 
+#include <iostream>
 #include <unordered_map>
 #include <vector>
 #include <cassert>
diff --git a/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.hpp b/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.hpp
index d90c6dc5a..ee2470ec6 100644
--- a/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuQuantumGeSvd_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUQUANTUMGESVD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuQuantumQr_internal.hpp b/src/backend/linalg_internal_gpu/cuQuantumQr_internal.hpp
index c2d6eaaf8..db43c63df 100644
--- a/src/backend/linalg_internal_gpu/cuQuantumQr_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuQuantumQr_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUQUANTUMQR_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuSub_internal.hpp b/src/backend/linalg_internal_gpu/cuSub_internal.hpp
index 8dc83d66d..9943ea974 100644
--- a/src/backend/linalg_internal_gpu/cuSub_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuSub_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUSUB_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuSum_internal.hpp b/src/backend/linalg_internal_gpu/cuSum_internal.hpp
index dbe009b38..30bb0a722 100644
--- a/src/backend/linalg_internal_gpu/cuSum_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuSum_internal.hpp
@@ -4,7 +4,6 @@
 #include "boost/smart_ptr/intrusive_ptr.hpp"
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuSvd_internal.hpp b/src/backend/linalg_internal_gpu/cuSvd_internal.hpp
index 2fae72edf..748d3e97e 100644
--- a/src/backend/linalg_internal_gpu/cuSvd_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuSvd_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUSVD_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuTensordot_internal.hpp b/src/backend/linalg_internal_gpu/cuTensordot_internal.hpp
index 7fa1176d9..56a41bc3f 100644
--- a/src/backend/linalg_internal_gpu/cuTensordot_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuTensordot_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUTENSORDOT_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuTrace_internal.hpp b/src/backend/linalg_internal_gpu/cuTrace_internal.hpp
index 13fc06fd9..c458c06ca 100644
--- a/src/backend/linalg_internal_gpu/cuTrace_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuTrace_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUTRACE_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cuVectordot_internal.hpp b/src/backend/linalg_internal_gpu/cuVectordot_internal.hpp
index 70726eb28..7421775c9 100644
--- a/src/backend/linalg_internal_gpu/cuVectordot_internal.hpp
+++ b/src/backend/linalg_internal_gpu/cuVectordot_internal.hpp
@@ -2,7 +2,6 @@
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUVECTORDOT_INTERNAL_H_
 
 #include <assert.h>
-#include <iostream>
 #include <iomanip>
 #include <vector>
 #include "backend/Storage.hpp"
diff --git a/src/backend/linalg_internal_gpu/cudaMemcpyTruncation.hpp b/src/backend/linalg_internal_gpu/cudaMemcpyTruncation.hpp
index 9d860df2b..a88c017e3 100644
--- a/src/backend/linalg_internal_gpu/cudaMemcpyTruncation.hpp
+++ b/src/backend/linalg_internal_gpu/cudaMemcpyTruncation.hpp
@@ -1,7 +1,6 @@
 #ifndef CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUDAMEMCPYTRUNCATION_H_
 #define CYTNX_BACKEND_LINALG_INTERNAL_GPU_CUDAMEMCPYTRUNCATION_H_
 
-#include <iostream>
 #include <vector>
 
 #include "Tensor.hpp"
diff --git a/src/backend/linalg_internal_interface.hpp b/src/backend/linalg_internal_interface.hpp
index 1b50a6e25..f369c79e8 100644
--- a/src/backend/linalg_internal_interface.hpp
+++ b/src/backend/linalg_internal_interface.hpp
@@ -1,7 +1,6 @@
 #ifndef CYTNX_BACKEND_LINALG_INTERNAL_INTERFACE_H_
 #define CYTNX_BACKEND_LINALG_INTERNAL_INTERFACE_H_
 
-#include <iostream>
 #include <vector>
 
 #include "Type.hpp"
diff --git a/src/linalg/Conj.cpp b/src/linalg/Conj.cpp
index 504be5b73..53eced9c1 100644
--- a/src/linalg/Conj.cpp
+++ b/src/linalg/Conj.cpp
@@ -1,5 +1,4 @@
 #include "linalg.hpp"
-#include <iostream>
 #include "Tensor.hpp"
 #include "UniTensor.hpp"
 
diff --git a/src/linalg/Conj_.cpp b/src/linalg/Conj_.cpp
index c0edca8e9..e2352e46c 100644
--- a/src/linalg/Conj_.cpp
+++ b/src/linalg/Conj_.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 #include "Tensor.hpp"
-#include <iostream>
 
 #ifdef BACKEND_TORCH
 
diff --git a/src/linalg/Cpr.cpp b/src/linalg/Cpr.cpp
index aad46dead..7c63d7654 100644
--- a/src/linalg/Cpr.cpp
+++ b/src/linalg/Cpr.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 
-#include <iostream>
 #include "Tensor.hpp"
 
 #ifdef BACKEND_TORCH
diff --git a/src/linalg/Det.cpp b/src/linalg/Det.cpp
index ae796c508..9373e8d8c 100644
--- a/src/linalg/Det.cpp
+++ b/src/linalg/Det.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 
-#include <iostream>
 #include "Tensor.hpp"
 
 #ifdef BACKEND_TORCH
diff --git a/src/linalg/Dot.cpp b/src/linalg/Dot.cpp
index 6dce19e53..666e05480 100644
--- a/src/linalg/Dot.cpp
+++ b/src/linalg/Dot.cpp
@@ -1,5 +1,4 @@
 #include "linalg.hpp"
-#include <iostream>
 #include <vector>
 #include "Tensor.hpp"
 
diff --git a/src/linalg/Gemm.cpp b/src/linalg/Gemm.cpp
index e2b6d4de6..a62f64d17 100644
--- a/src/linalg/Gemm.cpp
+++ b/src/linalg/Gemm.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 
-#include <iostream>
 #include "Tensor.hpp"
 #include "Generator.hpp"
 
diff --git a/src/linalg/Gemm_Batch.cpp b/src/linalg/Gemm_Batch.cpp
index a2f060cc4..0813ceaa4 100644
--- a/src/linalg/Gemm_Batch.cpp
+++ b/src/linalg/Gemm_Batch.cpp
@@ -1,5 +1,4 @@
 #include "linalg.hpp"
-#include <iostream>
 #include "Tensor.hpp"
 #include "Generator.hpp"
 
diff --git a/src/linalg/Hosvd.cpp b/src/linalg/Hosvd.cpp
index aa91b8641..6f1588f96 100644
--- a/src/linalg/Hosvd.cpp
+++ b/src/linalg/Hosvd.cpp
@@ -2,7 +2,6 @@
 
 #include "Tensor.hpp"
 #include "UniTensor.hpp"
-#include <iostream>
 #include <vector>
 
 #ifdef BACKEND_TORCH
diff --git a/src/linalg/Kron.cpp b/src/linalg/Kron.cpp
index 71b38d39e..293fc462c 100644
--- a/src/linalg/Kron.cpp
+++ b/src/linalg/Kron.cpp
@@ -1,7 +1,6 @@
 #include "linalg.hpp"
 #include "utils/utils.hpp"
 #include "Device.hpp"
-#include <iostream>
 #include "Tensor.hpp"
 
 #ifdef BACKEND_TORCH
diff --git a/src/linalg/Lanczos_Exp.cpp b/src/linalg/Lanczos_Exp.cpp
index c29031507..299c43a6a 100644
--- a/src/linalg/Lanczos_Exp.cpp
+++ b/src/linalg/Lanczos_Exp.cpp
@@ -1,15 +1,16 @@
-#include "linalg.hpp"
-#include "Generator.hpp"
-#include "random.hpp"
-#include "Tensor.hpp"
-#include "LinOp.hpp"
-
 #include <cfloat>
-#include <vector>
 #include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+#include "Generator.hpp"
+#include "LinOp.hpp"
+#include "Tensor.hpp"
 #include "UniTensor.hpp"
+#include "linalg.hpp"
+#include "random.hpp"
 #include "utils/vec_print.hpp"
-#include <iomanip>
 
 #ifdef BACKEND_TORCH
 #else
diff --git a/src/linalg/Lanczos_Gnd.cpp b/src/linalg/Lanczos_Gnd.cpp
index 8bec21e4f..c796257a7 100644
--- a/src/linalg/Lanczos_Gnd.cpp
+++ b/src/linalg/Lanczos_Gnd.cpp
@@ -1,13 +1,13 @@
-#include "linalg.hpp"
-#include "Generator.hpp"
-#include "random.hpp"
-#include "Tensor.hpp"
-#include "LinOp.hpp"
-
 #include <cfloat>
+#include <iomanip>
+#include <iostream>
 #include <vector>
+
+#include "Generator.hpp"
+#include "LinOp.hpp"
 #include "Tensor.hpp"
-#include <iomanip>
+#include "linalg.hpp"
+#include "random.hpp"
 
 #ifdef BACKEND_TORCH
 #else
diff --git a/src/linalg/Lanczos_Gnd_Ut.cpp b/src/linalg/Lanczos_Gnd_Ut.cpp
index 16af56b08..24066f9f4 100644
--- a/src/linalg/Lanczos_Gnd_Ut.cpp
+++ b/src/linalg/Lanczos_Gnd_Ut.cpp
@@ -1,15 +1,16 @@
-#include "linalg.hpp"
-#include "Generator.hpp"
-#include "random.hpp"
-#include "Tensor.hpp"
-#include "LinOp.hpp"
-
 #include <cfloat>
-#include <vector>
 #include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+#include "Generator.hpp"
+#include "LinOp.hpp"
+#include "Tensor.hpp"
 #include "UniTensor.hpp"
+#include "linalg.hpp"
+#include "random.hpp"
 #include "utils/vec_print.hpp"
-#include <iomanip>
 
 #ifdef BACKEND_TORCH
 #else
diff --git a/src/linalg/Lstsq.cpp b/src/linalg/Lstsq.cpp
index aa08b3182..41f6d6335 100644
--- a/src/linalg/Lstsq.cpp
+++ b/src/linalg/Lstsq.cpp
@@ -1,5 +1,4 @@
 #include "linalg.hpp"
-#include <iostream>
 #include "Tensor.hpp"
 #include "Generator.hpp"
 
diff --git a/src/linalg/Matmul.cpp b/src/linalg/Matmul.cpp
index d36c14f4c..b95f04416 100644
--- a/src/linalg/Matmul.cpp
+++ b/src/linalg/Matmul.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 
-#include <iostream>
 #include "Tensor.hpp"
 
 #ifdef BACKEND_TORCH
diff --git a/src/linalg/Matmul_dg.cpp b/src/linalg/Matmul_dg.cpp
index 6b7095389..e3f2fe300 100644
--- a/src/linalg/Matmul_dg.cpp
+++ b/src/linalg/Matmul_dg.cpp
@@ -1,5 +1,4 @@
 #include "linalg.hpp"
-#include <iostream>
 
 #ifdef BACKEND_TORCH
 #else
diff --git a/src/linalg/Norm.cpp b/src/linalg/Norm.cpp
index cf9aa849a..a91503cc9 100644
--- a/src/linalg/Norm.cpp
+++ b/src/linalg/Norm.cpp
@@ -1,5 +1,4 @@
 #include "linalg.hpp"
-#include <iostream>
 #include "Tensor.hpp"
 #include "cytnx.hpp"
 #ifdef BACKEND_TORCH
diff --git a/src/linalg/Outer.cpp b/src/linalg/Outer.cpp
index 08667390a..f13aa2761 100644
--- a/src/linalg/Outer.cpp
+++ b/src/linalg/Outer.cpp
@@ -2,7 +2,6 @@
 #include "utils/utils.hpp"
 #include "Device.hpp"
 #include "Tensor.hpp"
-#include <iostream>
 
 #ifdef BACKEND_TORCH
 #else
diff --git a/src/linalg/Qdr.cpp b/src/linalg/Qdr.cpp
index 9075cb05f..f43c1960d 100644
--- a/src/linalg/Qdr.cpp
+++ b/src/linalg/Qdr.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 #include "Tensor.hpp"
-#include <iostream>
 #include <vector>
 #include "UniTensor.hpp"
 #include "algo.hpp"
diff --git a/src/linalg/Rand_isometry.cpp b/src/linalg/Rand_isometry.cpp
index f844a0ea9..9a85e36f9 100644
--- a/src/linalg/Rand_isometry.cpp
+++ b/src/linalg/Rand_isometry.cpp
@@ -1,4 +1,3 @@
-#include <iostream>
 #include <string>
 #include <vector>
 
diff --git a/src/linalg/Tridiag.cpp b/src/linalg/Tridiag.cpp
index b0e74e1f2..d643b5f6e 100644
--- a/src/linalg/Tridiag.cpp
+++ b/src/linalg/Tridiag.cpp
@@ -1,6 +1,5 @@
 #include "linalg.hpp"
 
-#include <iostream>
 #include "Tensor.hpp"
 
 #ifdef BACKEND_TORCH
diff --git a/src/search_tree.cpp b/src/search_tree.cpp
index c0e156753..f5909e309 100644
--- a/src/search_tree.cpp
+++ b/src/search_tree.cpp
@@ -1,4 +1,6 @@
 #include "search_tree.hpp"
+
+#include <iostream>
 #include <stack>
 
 using namespace std;
diff --git a/src/stat/histogram.cpp b/src/stat/histogram.cpp
index 3dc718dda..d1b83b382 100644
--- a/src/stat/histogram.cpp
+++ b/src/stat/histogram.cpp
@@ -35,11 +35,11 @@ namespace cytnx {
       this->vars = tmp.storage();
     }
 
-    void Histogram::print() const {
-      std::cout << "[Histogram 1D] Real" << std::endl;
-      std::cout << "Nbins: " << this->bins << std::endl;
-      std::cout << "bound: [ " << this->min << " , " << this->max << " ]\n";
-      std::cout << "current count: " << this->total_count << std::endl;
+    void Histogram::print(std::ostream &os) const {
+      os << "[Histogram 1D] Real" << std::endl;
+      os << "Nbins: " << this->bins << std::endl;
+      os << "bound: [ " << this->min << " , " << this->max << " ]\n";
+      os << "current count: " << this->total_count << std::endl;
     }
 
     //-----------[2d]
@@ -83,12 +83,12 @@ namespace cytnx {
       this->vars = tmp.storage();
     }
 
-    void Histogram2d::print() const {
-      std::cout << "[Histogram 2D] Real" << std::endl;
-      std::cout << "Nbins: [x= " << this->binx << " , y= " << this->biny << " ]\n";
-      std::cout << "bound,x: [ " << this->minx << " , " << this->maxx << " ]\n";
-      std::cout << "bound,y: [ " << this->miny << " , " << this->maxy << " ]\n";
-      std::cout << "current count: " << this->total_count << std::endl;
+    void Histogram2d::print(std::ostream &os) const {
+      os << "[Histogram 2D] Real" << std::endl;
+      os << "Nbins: [x= " << this->binx << " , y= " << this->biny << " ]\n";
+      os << "bound,x: [ " << this->minx << " , " << this->maxx << " ]\n";
+      os << "bound,y: [ " << this->miny << " , " << this->maxy << " ]\n";
+      os << "current count: " << this->total_count << std::endl;
     }
 
   }  // namespace stat
diff --git a/src/tn_algo/DMRG.cpp b/src/tn_algo/DMRG.cpp
index 8eba8d606..678e25c7b 100644
--- a/src/tn_algo/DMRG.cpp
+++ b/src/tn_algo/DMRG.cpp
@@ -1,12 +1,14 @@
-#include "tn_algo/DMRG.hpp"
+#include <iomanip>
+#include <iostream>
+#include <tuple>
+
 #include "Generator.hpp"
-#include "Network.hpp"
 #include "LinOp.hpp"
+#include "Network.hpp"
 #include "linalg.hpp"
-#include <tuple>
-#include <iomanip>
-#include <iostream>
+#include "tn_algo/DMRG.hpp"
 #include "utils/vec_print.hpp"
+
 using namespace std;
 
 #ifdef BACKEND_TORCH
diff --git a/src/tn_algo/MPS.cpp b/src/tn_algo/MPS.cpp
index e94e360c6..5c627e905 100644
--- a/src/tn_algo/MPS.cpp
+++ b/src/tn_algo/MPS.cpp
@@ -2,7 +2,7 @@
 
 #include <filesystem>
 #include <fstream>
-#include <iostream>
+#include <ostream>
 
 using namespace std;
 
diff --git a/src/utils/cutensornet.cu b/src/utils/cutensornet.cu
index 34f898f14..a3b0bd45d 100644
--- a/src/utils/cutensornet.cu
+++ b/src/utils/cutensornet.cu
@@ -1,9 +1,13 @@
-#include <stdlib.h>
 #include <stdio.h>
+#include <stdlib.h>
+
+#include <cassert>
+#include <iostream>
 #include <unordered_map>
 #include <vector>
-#include <cassert>
+
 #include <cytnx.hpp>
+
 #include "cutensornet.hpp"
 
 #ifdef UNI_GPU
diff --git a/src/utils/utils.cpp b/src/utils/utils.cpp
index 4d5258b28..d7696eee0 100644
--- a/src/utils/utils.cpp
+++ b/src/utils/utils.cpp
@@ -1,6 +1,5 @@
 #include "utils/utils.hpp"
 #include <vector>
-#include <iostream>
 #include <string>
 namespace cytnx {