two_site_split.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <uni20/symmetry/block_tensor_contract.hpp>
#include <uni20/symmetry/block_tensor_linear.hpp>
#include <uni20/symmetry/block_tensor_repartition.hpp>
#include <uni20/symmetry/block_tensor_space_traits.hpp>
#include <uni20/symmetry/block_tensor_svd.hpp>
#include <uni20/tensor_network/finite_chain.hpp>
 
#include <concepts>
#include <cstddef>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace uni20::tensor_network
{
 
enum class MpsSweepDirection
{
  left_to_right, 
  right_to_left  
};
 
template <class FirstSite, class SecondSite, class SingularValues, class Truncation> struct TwoSiteMpsSplit
{
    FirstSite first_site;
    SecondSite second_site;
    SingularValues singular_values;
    Truncation truncation;
    MpsSweepDirection direction;
};
 
template <class SingularValues, class Truncation> struct InstalledMpsBond
{
    SingularValues singular_values;
    Truncation truncation;
    MpsSweepDirection direction;
};
 
namespace detail
{
 
template <class Space>
concept SplitBondSpace =
    BlockTensorSpace<Space> && BlockTensorSpaceTraits<std::remove_cvref_t<Space>>::has_block_coordinate &&
    BlockTensorSpaceTraits<std::remove_cvref_t<Space>>::has_dense_axis;
 
template <class Space>
concept SplitPhysicalSpace =
    BlockTensorSpace<Space> && BlockTensorSpaceTraits<std::remove_cvref_t<Space>>::has_block_coordinate &&
    !BlockTensorSpaceTraits<std::remove_cvref_t<Space>>::has_dense_axis;
 
template <class Storage, class Scalar>
concept MpsSiteStorage = SparseBlockStorage<Storage> && !DiagonalBlockStorage<Storage> &&
                         LocalBlockStorageFor<Storage, std::remove_cv_t<Scalar>, 3, 2> &&
                         !AsyncLocalBlockStorageFor<Storage, std::remove_cv_t<Scalar>, 3, 2>;
 
template <class Tensor, std::size_t Index>
using split_domain_space_t = typename block_tensor_domain_t<Tensor>::template space_type<Index>;
 
template <class Tensor, std::size_t Index>
using split_codomain_space_t = typename block_tensor_codomain_t<Tensor>::template space_type<Index>;
 
template <class Tensor>
concept TwoSiteSvdCenter =
    BlockTensorView<Tensor> && (block_tensor_domain_t<Tensor>::size() == 3) &&
    (block_tensor_codomain_t<Tensor>::size() == 1) && (block_tensor_type_t<Tensor>::key_coordinate_count() == 4) &&
    (block_tensor_type_t<Tensor>::dense_block_order() == 2) && SplitBondSpace<split_domain_space_t<Tensor, 0>> &&
    SplitPhysicalSpace<split_domain_space_t<Tensor, 1>> && SplitPhysicalSpace<split_domain_space_t<Tensor, 2>> &&
    SplitBondSpace<split_codomain_space_t<Tensor, 0>>;
 
template <SparseBlockStorage OutputStorage, BlockTensorView Source>
auto materialize_block_tensor_copy(Source const& source)
{
  using output_type = BlockTensor<block_tensor_value_t<Source>, block_tensor_domain_t<Source>,
                                  block_tensor_codomain_t<Source>, OutputStorage>;
  using key_type = typename output_type::key_type;
  std::vector<key_type> keys(source.stored_keys().begin(), source.stored_keys().end());
  auto result = [&] {
    if constexpr (requires {
                    output_type(uninitialized, source.symmetry(), source.domain(), source.codomain(), keys,
                                source.allocation_context());
                  })
      return output_type(uninitialized, source.symmetry(), source.domain(), source.codomain(), std::move(keys),
                         source.allocation_context());
    else
      return output_type(uninitialized, source.symmetry(), source.domain(), source.codomain(), std::move(keys));
  }();
  for (std::size_t ordinal = 0; ordinal < result.stored_block_count(); ++ordinal)
  {
    auto output_block = result.block_by_ordinal(ordinal);
    auto input_block = source.block_by_ordinal(ordinal);
    uni20::copy(output_block, input_block);
  }
  return result;
}
 
} // namespace detail
 
template <detail::TwoSiteSvdCenter Center, class Measurements, class Event>
  requires uni20::LapackScalar<block_tensor_value_t<Center>> && performance::BatchMeasurementPolicy<Measurements, Event>
[[nodiscard]] auto decompose_two_site_center(Center const& center, linalg::SvdOptions options,
                                             Measurements& measurements, Event batch_event)
{
  auto matrix = repartition<MorphismSide::Domain, BoundaryEnd::Right>(center);
  return block_svd(matrix, options, measurements, batch_event);
}
 
template <detail::TwoSiteSvdCenter Center>
  requires uni20::LapackScalar<block_tensor_value_t<Center>>
[[nodiscard]] auto decompose_two_site_center(Center const& center, linalg::SvdOptions options = {})
{
  performance::NoMeasurements measurements;
  return decompose_two_site_center(center, options, measurements, nullptr);
}
 
template <SparseBlockStorage OutputStorage = SeparateSparseBlockStorage<>, class Decomposition>
  requires detail::MpsSiteStorage<OutputStorage, typename Decomposition::scalar_type>
[[nodiscard]] auto materialize_two_site_mps_split(Decomposition const& decomposition,
                                                  BlockSvdSelection<typename Decomposition::real_type> const& selection,
                                                  MpsSweepDirection direction,
                                                  BlockSvdMaterializationOptions options = {})
{
  static_assert(Decomposition::domain_type::size() == 2, "two-site MPS decomposition requires a two-factor domain");
  static_assert(Decomposition::codomain_type::size() == 2, "two-site MPS decomposition requires a two-factor codomain");
  if (selection.state_ids().empty()) throw std::invalid_argument("two-site MPS split requires a retained state");
  if (direction != MpsSweepDirection::left_to_right && direction != MpsSweepDirection::right_to_left)
    throw std::invalid_argument("two-site MPS split has an invalid sweep direction");
 
  bool const absorb_left = direction == MpsSweepDirection::left_to_right;
  bool const absorb_right = direction == MpsSweepDirection::right_to_left;
  auto factors = uni20::detail::materialize_svd_with_absorption(decomposition, selection, std::move(options),
                                                                absorb_left, absorb_right);
  if (direction == MpsSweepDirection::left_to_right)
  {
    auto first_site = detail::materialize_block_tensor_copy<OutputStorage>(factors.right_singular_vectors_adjoint);
    auto second_view = repartition<MorphismSide::Codomain, BoundaryEnd::Right>(factors.left_singular_vectors);
    auto second_site = detail::materialize_block_tensor_copy<OutputStorage>(second_view);
    return TwoSiteMpsSplit{.first_site = std::move(first_site),
                           .second_site = std::move(second_site),
                           .singular_values = std::move(factors.singular_values),
                           .truncation = std::move(factors.truncation),
                           .direction = direction};
  }
 
  auto first_site = detail::materialize_block_tensor_copy<OutputStorage>(factors.right_singular_vectors_adjoint);
  auto second_view = repartition<MorphismSide::Codomain, BoundaryEnd::Right>(factors.left_singular_vectors);
  auto second_site = detail::materialize_block_tensor_copy<OutputStorage>(second_view);
  return TwoSiteMpsSplit{.first_site = std::move(first_site),
                         .second_site = std::move(second_site),
                         .singular_values = std::move(factors.singular_values),
                         .truncation = std::move(factors.truncation),
                         .direction = direction};
}
 
template <typename Scalar, Space Bond, Space Physical, BlockTensorStorage SiteStorage, class Decomposition>
  requires std::same_as<Bond, BlockSpace> && detail::MpsSiteStorage<SiteStorage, Scalar>
[[nodiscard]] auto replace_two_site_from_svd(FiniteMps<Scalar, Bond, Physical, SiteStorage>& mps,
                                             std::size_t first_index, Decomposition const& decomposition,
                                             BlockSvdSelection<typename Decomposition::real_type> const& selection,
                                             MpsSweepDirection direction, BlockSvdMaterializationOptions options = {})
{
  auto split = materialize_two_site_mps_split<SiteStorage>(decomposition, selection, direction, std::move(options));
  mps.replace_pair(first_index, std::move(split.first_site), std::move(split.second_site));
  return InstalledMpsBond{.singular_values = std::move(split.singular_values),
                          .truncation = std::move(split.truncation),
                          .direction = direction};
}
 
} // namespace uni20::tensor_network