Realization of SVD compression for MPO
======================================

From the construction of the SVD based compression scheme for MPO, we try to understand the convention used for the tensor library. This could serve as an example on how you develop your own MPO/MPS operation based on the tensor library.

In :code:`mpo_impl.h`, the class mpo_impl contains the following functions

.. code-block :: c++ 
  :linenos: 
  
   std::pair<size_t, magnitude_type> svdCompressLeft(int maxbonddimension = -1, magnitude_type tol = 1.0e-12, int verbose = 0);

   std::pair<size_t, magnitude_type> svdCompressRight(int maxbonddimension = -1, magnitude_type tol = 1.0e-12, int verbose = 0);

   std::pair<size_t, magnitude_type> svdCompress(int maxbonddimension=1000000, magnitude_type tol = 1.0e-12, int verbose = 0);


Function :code:`svdCompress` would call :code:`svdCompressLeft` and :code:`svdCompressRight`.

.. code-block :: c++ 
  :linenos: 

   template <typename _Tensor4>
   std::pair<size_t, typename mpo_impl<_Tensor4>::magnitude_type>
   mpo_impl<_Tensor4>::svdCompress(int maxbonddimension, magnitude_type tol, int verbose)
   {
       svdCompressLeft(-1, tol, verbose);
       return svdCompressRight(maxbonddimension, tol, verbose);
   }




.. code-block :: c++
  :linenos:

   template <typename _Tensor4>
   std::pair<size_t, typename mpo_impl<_Tensor4>::magnitude_type>
   mpo_impl<_Tensor4>::svdCompressLeft(int maxbonddimension, magnitude_type tol, int verbose)
   {
       size_t L = super_type::size(), i;
       assert(L > 1);
       typedef std::pair<size_t, magnitude_type> bonderror_type;
       std::vector<bonderror_type> bonderror;
       std::array<size_t, 1> Is2 = {2}, Is1 = {1}, Is0 = {0};
       std::array<size_t, 4> per1 = {0, 1, 3, 2}, per2 = {1, 0, 2, 3};
       for (i = 0; i < L - 1; ++i) {
           if (verbose >= 2) {
               std::cout << "Compress mpo by balanced svd from left to right on site: " << i << '\n';
           }
           auto svdResults = gc::svdCompress(operator[](i), Is2, maxbonddimension, tol, verbose);
           bonderror.push_back(std::get<3>(svdResults));
           if (std::get<1>(svdResults).size() == 0) {
               if (verbose) {
                   std::cout << "mpo becomes zero after cutoff.\n";
               }
               for (size_t j = 0; j < L; ++j) {
                   operator[](j).clear();
               }
               break;
           }
           typename std::remove_reference<decltype(std::get<2>(svdResults))>::type sm;
           diag(std::get<1>(svdResults), sm);
           std::get<2>(svdResults) = contract(sm, Is1, std::get<2>(svdResults), Is0);
           // operator[](i) = contract(std::get<0>(svdResults), Is3, sm, Is0);
           operator[](i) = gc::permute(std::get<0>(svdResults), per1);
           operator[](i + 1) = contract(std::get<2>(svdResults), Is1, operator[](i + 1), Is1);
           operator[](i + 1) = gc::permute(operator[](i + 1), per2);
       }
       bonderror_type bet = bonderror[0];
       for (i = 1; i < bonderror.size(); ++i) {
           bet = std::make_pair(std::max(bet.first, bonderror[i].first), bet.second + bonderror[i].second);
       }
       return bet;
   }