https://github.com/ITensor/ITensor
Tip revision: 6b54fa39438c296f80049dc937ca966746eded68 authored by Miles Stoudenmire on 11 April 2016, 18:02:29 UTC
Updated localop.h diag method to correctly use delta function to tie indices. Fixes #96 - thanks to Xiongjie Yu for reporting it.
Updated localop.h diag method to correctly use delta function to tie indices. Fixes #96 - thanks to Xiongjie Yu for reporting it.
Tip revision: 6b54fa3
itensor.ih
//
// Distributed under the ITensor Library License, Version 1.2
// (See accompanying LICENSE file.)
//
#ifndef __ITENSOR_ITENSOR_IH_
#define __ITENSOR_ITENSOR_IH_
//
// Template Method Implementations
//
namespace itensor {
//template<>
//template <typename... IVals>
//ITensor::
//ITensorT(const IndexVal& iv1,
// const IVals&... rest)
// : scale_(1.)
// {
// const size_t size = 1+sizeof...(rest);
// auto ivs = std::array<IndexVal,size>{{iv1,rest...}};
// std::array<Index,size> inds;
// for(size_t j = 0; j < size; ++j) inds[j] = ivs[j].index;
// is_ = IndexSet(inds);
// store_ = newITData<DenseReal>(area(is_),0.);
// set(iv1,rest...,1.);
// }
// TODO: implement special code for doing this operation
// more efficiently
inline ITensor&
operator*=(ITensor & T, IndexVal const& iv) { return T *= setElt(iv); }
ITensor inline
operator*(IndexVal const& iv, ITensor const& B)
{
auto A = setElt(iv);
A *= B;
return A;
}
ITensor inline
operator*(IndexVal const& iv1, IndexVal const& iv2)
{
auto t = setElt(iv1);
return (t *= iv2);
}
ITensor inline
operator*(IndexVal const& iv1, Real val)
{
auto res = setElt(iv1);
res *= val;
return res;
}
template<typename... Inds>
ITensor
delta(Index const& i1,
Inds const&... inds)
{
auto is = IndexSet(i1,inds...);
auto len = minM(is);
return ITensor(std::move(is),DiagReal(len,1.));
}
template<typename Container, typename... Inds, class>
ITensor
diagTensor(Container const& C,
Index const& i1,
Inds &&... inds)
{
auto is = IndexSet(i1,std::forward<Inds>(inds)...);
#ifdef DEBUG
using size_type = decltype(C.size());
//Compute min of all index dimensions
auto minm = i1.m();
for(const auto& ind : is)
if(ind.m() < minm) minm = ind.m();
if(C.size() != size_type(minm))
{
println("minm = ",minm);
println("C.size() = ",C.size());
Error("Wrong size of data in diagonal ITensor constructor");
}
#endif
using value_type = typename Container::value_type;
return ITensor(std::move(is),Diag<value_type>(C.begin(),C.end()));
}
template<typename IndexT>
ITensorT<IndexT>
randomTensor(const IndexSetT<IndexT>& inds)
{
return random(ITensorT<IndexT>{inds});
}
template<typename V, typename... Indxs>
TenRef<Range1,V>
ordered_impl(ITensor & T, Indxs&&... inds)
{
T.scaleTo(1.);
auto getData = [](Dense<V> & d) -> DataRange<V>
{
return makeDataRange(d.data(),d.size());
};
auto d = applyFunc(getData,T.store());
constexpr auto size = sizeof...(inds);
auto indarr = std::array<Index,size>{{static_cast<Index>(inds)...}};
auto P = Permutation(T.r());
calcPerm(T.inds(),indarr,P);
return makeRef(d,permuteRangeTo<Range1>(T.inds(),P));
}
template<typename... Indxs>
TensorRef1
ordered(ITensor & T, Indxs&&... inds)
{
if(isComplex(T)) Error("ITensor is complex; use orderedC instead of ordered");
return ordered_impl<Real>(T,std::forward<Indxs>(inds)...);
}
template<typename... Indxs>
CTensorRef1
orderedC(ITensor & T, Indxs&&... inds)
{
if(!isComplex(T)) Error("ITensor is real; use ordered instead of orderedC");
return ordered_impl<Cplx>(T,std::forward<Indxs>(inds)...);
}
} //namespace itensor
#endif