rsvd_impl.hpp

Go to the documentation of this file.

/*
  mptensor - Parallel Library for Tensor Network Methods
 
  Copyright 2016 Satoshi Morita
 
  mptensor is free software: you can redistribute it and/or modify it
  under the terms of the GNU Lesser General Public License as
  published by the Free Software Foundation, either version 3 of the
  License, or (at your option) any later version.
 
  mptensor is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
 
  You should have received a copy of the GNU Lesser General Public
  License along with mptensor.  If not, see
  <https://www.gnu.org/licenses/>.
*/
 
#ifndef _TENSOR_RSVD_IMPL_HPP_
#define _TENSOR_RSVD_IMPL_HPP_
 
#include <cassert>
#include <vector>
 
#include "complex.hpp"
#include "index.hpp"
#include "matrix.hpp"
#include "rsvd.hpp"
#include "tensor.hpp"
 
namespace mptensor {
 
namespace random_tensor {
 
template <typename C>
C uniform_dist();
void set_seed(unsigned int seed);
 
template <typename tensor_t>
void fill(tensor_t &t) {
  const size_t n = t.local_size();
  for (size_t i = 0; i < n; ++i) {
    t[i] = uniform_dist<typename tensor_t::value_type>();
  }
}
 
}  // namespace random_tensor
 
 
template <template <typename> class Matrix, typename C>
int rsvd(const Tensor<Matrix, C> &a, const Axes &axes_row, const Axes &axes_col,
         Tensor<Matrix, C> &u, std::vector<double> &s, Tensor<Matrix, C> &vt,
         const size_t target_rank, const size_t oversamp) {
  assert(axes_row.size() > 0);
  assert(axes_col.size() > 0);
  assert(debug::check_svd_axes(axes_row, axes_col, a.rank()));
 
  const size_t rank = a.rank();
  const size_t rank_row = axes_row.size();
  const size_t rank_col = axes_col.size();
  int info;
 
  Axes axes = axes_row + axes_col;
  Tensor<Matrix, C> a_t = transpose(a, axes, rank_row);
  const Shape &shape = a_t.shape();
 
  Tensor<Matrix, C> q;
  {
    Shape shape_omega;
    shape_omega.resize(rank_col + 1);
    for (size_t i = 0; i < rank_col; ++i) shape_omega[i] = shape[i + rank_row];
    shape_omega[rank_col] = target_rank + oversamp;
 
    // Tensor<Matrix,C> omega = random_tensor<C>(a.get_comm(), shape_omega);
    Tensor<Matrix, C> omega(a.get_comm(), shape_omega, rank_col);
    random_tensor::fill(omega);
 
    Tensor<Matrix, C> r;
    qr(tensordot(a_t, omega, range(rank_row, rank), range(0, rank_col)),
       range(0, rank_row), Axes(rank_row), q, r);
  }
 
  info = svd(tensordot(conj(q), a_t, range(0, rank_row), range(0, rank_row)),
             Axes(0), range(1, rank_col + 1), u, s, vt);
  s.resize(target_rank);
  u = slice(u, 1, 0, target_rank);
  vt = slice(vt, 0, 0, target_rank);
  u = tensordot(q, u, Axes(rank_row), Axes(0));
  return info;
};
 
template <template <typename> class Matrix, typename C>
int rsvd(const Tensor<Matrix, C> &a, const Axes &axes_row, const Axes &axes_col,
         Tensor<Matrix, C> &u, std::vector<double> &s, Tensor<Matrix, C> &vt,
         const size_t target_rank) {
  return rsvd(a, axes_row, axes_col, u, s, vt, target_rank, target_rank);
};
 
 
template <template <typename> class Matrix, typename C, typename Func1,
          typename Func2>
int rsvd(Func1 &multiply_row, Func2 &multiply_col, const Shape &shape_row,
         const Shape &shape_col, Tensor<Matrix, C> &u, std::vector<double> &s,
         Tensor<Matrix, C> &vt, const size_t target_rank,
         const size_t oversamp) {
  const size_t rank_row = shape_row.size();
  const size_t rank_col = shape_col.size();
 
  int info;
  Tensor<Matrix, C> q;
  {
    Shape shape_omega = shape_col;
    shape_omega.resize(rank_col + 1);
    shape_omega[rank_col] = target_rank + oversamp;
 
    Tensor<Matrix, C> omega(u.get_comm(), shape_omega, rank_col);
    random_tensor::fill(omega);
 
    Tensor<Matrix, C> r;
    qr(multiply_col(omega), range(0, rank_row), Axes(rank_row), q, r);
  }
 
  info = svd(multiply_row(conj(q)), Axes(0), range(1, rank_col + 1), u, s, vt);
  s.resize(target_rank);
  u = slice(u, 1, 0, target_rank);
  vt = slice(vt, 0, 0, target_rank);
  u = tensordot(q, u, Axes(rank_row), Axes(0));
 
  return info;
}
 
template <template <typename> class Matrix, typename C, typename Func1,
          typename Func2>
int rsvd(Func1 &multiply_row, Func2 &multiply_col, const Shape &shape_row,
         const Shape &shape_col, Tensor<Matrix, C> &u, std::vector<double> &s,
         Tensor<Matrix, C> &vt, const size_t target_rank) {
  return rsvd(multiply_row, multiply_col, shape_row, shape_col, u, s, vt,
              target_rank, target_rank);
};
 
 
inline void set_seed(unsigned int seed) { random_tensor::set_seed(seed); };
 
}  // namespace mptensor
 
#endif  // _TENSOR_RSVD_IMPL_HPP_