ITensor/itensor/index.cc at unordered_map · ITensor/ITensor

779 lines (657 loc) · 15.3 KB
// Copyright 2018 The Simons Foundation, Inc. - All Rights Reserved.
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//    http://www.apache.org/licenses/LICENSE-2.0
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "itensor/index.h"
#include "itensor/util/readwrite.h"
#include "itensor/util/print_macro.h"
namespace itensor {
using std::string;
using std::stringstream;
putprimes(string s, int plev)
    stringstream str;
    str << s;
    if(plev < 0) Error("Negative prime level");
    if(plev > 3)
        str << "'" << plev;
        for(int i = 1; i <= plev; ++i) str << "\'";
    return str.str();
// class Index
Index::id_type Index::
generateID()
    static thread_local Index::IDGenerator G;
    return G();
    id_(0),
    dim_(1),
    tags_(TagSet("0"))
    if(primeLevel() < 0) setPrime(0);
Index(long m, 
      TagSet const& t)
  : id_(generateID()),
    dim_(m),
    tags_(t)
    if(primeLevel() < 0) setPrime(0);
Index(id_type id,
      long dim, 
      Arrow dir, 
      TagSet const& ts)
  : id_(id),
    dim_(dim),
    dir_(dir),
    tags_(ts)
    if(primeLevel() < 0) setPrime(0);
Index& Index::
setPrime(int plev) 
    tags_.setPrime(plev);
#ifdef DEBUG
    if(this->primeLevel() < 0)
        Error("Negative primeLevel");
    return *this;
Index& Index::
noPrime()  
    tags_.setPrime(0);
    return *this;
Index& Index::
prime(int inc) 
    tags_.prime(inc);
#ifdef DEBUG
    if(this->primeLevel() < 0)
        Error("Negative primeLevel");
    return *this;
operator bool() const { return (id_!=0); }
IndexVal Index::
operator()(long val) const
    return IndexVal(*this,val);
IndexVal Index::
operator=(long val) const 
    return operator()(val); 
operator==(Index const& i1, Index const& i2)
    return (id(i1) == id(i2)) && (tags(i1) == tags(i2));
operator!=(Index const& i1, Index const& i2)
    return not operator==(i1,i2);
Index::id_type
id(Index const& I) { return I.id(); }
dim(Index const& I) { return I.dim(); }
dim(IndexVal const& I) { return I.dim(); }
primeLevel(Index const& I) { return I.primeLevel(); }
tags(const Index& I) { return I.tags(); }
addTags(Index I, TagSet const& t) { I.addTags(t); return I; }
removeTags(Index I, TagSet const& t) { I.removeTags(t); return I; }
setTags(Index I, TagSet const& t) { I.setTags(t); return I; }
noTags(Index I) { I.noTags(); return I; }
replaceTags(Index I, TagSet const& tsold, TagSet const& tsnew) { I.replaceTags(tsold,tsnew); return I; }
prime(Index I, int plinc) { I.prime(plinc); return I; }
setPrime(Index I, int plev) { I.setPrime(plev); return I; }
noPrime(Index I) { I.noPrime(); return I; }
addTags(IndexVal I, TagSet const& t) { I.addTags(t); return I; }
removeTags(IndexVal I, TagSet const& t) { I.removeTags(t); return I; }
setTags(IndexVal I, TagSet const& t) { I.setTags(t); return I; }
noTags(IndexVal I) { I.noTags(); return I; }
replaceTags(IndexVal I, TagSet const& tsold, TagSet const& tsnew) { I.replaceTags(tsold,tsnew); return I; }
prime(IndexVal I, int plinc) { I.prime(plinc); return I; }
setPrime(IndexVal I, int plev) { I.setPrime(plev); return I; }
noPrime(IndexVal I) { I.noPrime(); return I; }
// Check if Index I contains the tags tsmatch.
hasTags(Index I, const TagSet& tsmatch) { return hasTags(tags(I),tsmatch); }
hasQNs(Index const& I) { return nblock(I)!=0; }
removeQNs(Index I) { if(hasQNs(I)) I.removeQNs(); return I; }
Index const&
index(IndexVal const& res) { return res.index; }
val(IndexVal const& res) { return res.val; }
hasQNs(IndexVal const& iv) { return hasQNs(index(iv)); }
dag(Index res) { res.dag(); return res; }
dag(IndexVal res) { res.dag(); return res; }
qn(IndexVal const& iv) { return iv.qn(); }
dir(Index const& res) { return res.dir(); }
dir(IndexVal const& res) { return dir(index(res)); }
//TODO: what is this for? It doesn't make as much sense with
//tags, but I guess we can compare tags with inequalities
operator>(Index const& i1, Index const& i2)
    if(dim(i1) == dim(i2)) 
        if(id(i1) == id(i2)) return primeLevel(i1) > primeLevel(i2);
        return id(i1) > id(i2);
    return dim(i1) > dim(i2);
operator<(Index const& i1, Index const& i2)
    if(dim(i1) == dim(i2)) 
        if(id(i1) == id(i2)) return primeLevel(i1) < primeLevel(i2);
        return id(i1) < id(i2);
    return dim(i1) < dim(i2);
std::ostream& 
operator<<(std::ostream & s, Index const& I)
    s << "(";
    s << dim(I);
    if(Global::showIDs()) 
        s << "|id=" << (id(I) % 1000);
        //s << "," << id(I);
    auto ts = tags(I);
    if(size(ts) > 0)
      s << "|";
      auto ts = tags(I);
      for(auto i : range(size(ts)))
        s << ts[i];
        if( i < (size(ts)-1) ) s << ",";
    s << ")"; 
    if(primeLevel(I) > 0) 
        if(primeLevel(I) > 3)
            s << "'" << primeLevel(I);
        else
            for(int n = 1; n <= primeLevel(I); ++n)
                s << "'";
    else if(primeLevel(I) < 0) 
        s << "'" << primeLevel(I) << " (WARNING: prime level of Index is negative, not well defined)";
    if(hasQNs(I))
        s << " <" << I.dir() << ">\n";
        for(auto j : range1(nblock(I)))
            s << "  " << j << ": " << blocksize(I,j) << " " <<  qn(I,j);
            if(j != nblock(I)) s << "\n";
    return s;
IndexVal() 
    val(0) 
IndexVal(Index const& index_, 
         long val_) 
    index(index_),
    val(val_)
#ifdef DEBUG
    if(!index) Error("IndexVal initialized with default initialized Index");
    //Can also use IndexVal's to indicate prime increments:
    //if(val_ < 1 || val_ > dim(index))
    //    {
    //    println("val = ",val_);
    //    println("index = ",index);
    //    Error("val out of range");
    //    }
operator==(IndexVal const& iv1, IndexVal const& iv2)
    return (index(iv1) == index(iv2) && iv1.val == iv2.val);
operator!=(IndexVal const& iv1, IndexVal const& iv2)
    return not operator==(iv1,iv2);
operator==(Index const& I, IndexVal const& iv)
    return index(iv) == I;
operator==(IndexVal const& iv, Index const& I)
    return index(iv) == I;
showDim(Index const& I) { return format("dim=%d",dim(I)); }
std::ostream& 
operator<<(std::ostream& s, IndexVal const& iv)
    return s << "IndexVal: val = " << iv.val 
             << ", ind = " << index(iv);
add(Args            & args, 
    Args::Name const& name, 
    TagSet     const& ts) 
    args.add(name,ts); 
getTagSet(Args       const& args, 
          Args::Name const& name)
    if(!args.defined(name)) Error(format("Name %s not found in Args",name));
    return TagSet(args.getString(name));
getTagSet(const Args& args, 
          const Args::Name& name, 
          TagSet const& default_val)
    if(!args.defined(name)) return default_val; 
    return TagSet(args.getString(name));
struct IndSector
    long sector = 0l;
    long sind   = 0l;
    IndSector(long sec, long si) : sector(sec), sind(si) { }
sectorInfo(IndexVal const& iv)
    auto is = IndSector(1,iv.val);
    while(is.sind > blocksize(index(iv),is.sector))
        is.sind -= blocksize(index(iv),is.sector);
        is.sector += 1;
    return is;
QN const& IndexVal::
qn() const 
    auto is = sectorInfo(*this);
    return index.qn(is.sector);
IndexVal&  IndexVal::
dag() { index.dag(); return *this; }
class IQIndexDat
    public:
    using storage = std::vector<QNInt>;
    using iterator = storage::iterator;
    using const_iterator = storage::const_iterator;
    private:
    storage iq_;
    public:
    IQIndexDat() { }
    explicit
    IQIndexDat(storage const& ind_qn) 
      : iq_(ind_qn)
    explicit
    IQIndexDat(storage&& ind_qn) 
      : iq_(std::move(ind_qn))
    //Disallow copying
    IQIndexDat(IQIndexDat const&) = delete;
    operator=(IQIndexDat const&) = delete;
    setStore(storage && iq) { iq_ = std::move(iq); }
    storage const&
    inds() const { return iq_; }
    size() const { return iq_.size(); }
    blocksize(long i) { return iq_[i-1].second; }
    blocksize0(long i) { return iq_[i].second; }
    QN const&
    qn(long i) { return iq_[i-1].first; }
    iterator
    begin() { return iq_.begin(); }
    iterator
    end() { return iq_.end(); }
    const_iterator
    begin() const { return iq_.begin(); }
    const_iterator
    end()   const { return iq_.end(); }
    storage const&
    store() const { return iq_; }
#ifdef DEBUG
#define IQINDEX_CHECK_NULL if(pd == 0) throw ITError("IQIndex storage unallocated");
#define IQINDEX_CHECK_NULL
directSum(Index const& i,
          Index const& j,
          Args const& args)
  auto tags = getTagSet(args,"Tags","sum");
  if(not hasQNs(i) && not hasQNs(j))
    auto dim_ij = dim(i) + dim(j);
    if(dim_ij <= 0) dim_ij = 1;
    return Index(dim_ij,tags);
#ifdef DEBUG
    if( dir(i) != dir(j) ) Error("In directSum(Index, Index), input indices must have same arrow direction");
    auto nblock_i = nblock(i);
    auto nblock_j = nblock(j);
    auto siq = stdx::reserve_vector<QNInt>(nblock_i+nblock_j);
    for(auto iqn : range1(nblock_i))
        siq.emplace_back(qn(i,iqn),blocksize(i,iqn));
    for(auto jqn : range1(nblock_j))
        siq.emplace_back(qn(j,jqn),blocksize(j,jqn));
#ifdef DEBUG
    if(siq.empty()) Error("siq is empty in plussers");
    return Index(std::move(siq),dir(i),tags);
  return Index();
sim(Index const& I)
    Index J;
    J.sim(I);
    return J;
write(std::ostream & s, QNInt const& q)
    write(s,q.first);
    write(s,q.second);
read(std::istream & s, QNInt & q)
    read(s,q.first);
    read(s,q.second);
write(std::ostream & s, IQIndexDat const& d) 
    write(s,d.store()); 
read(std::istream & s, IQIndexDat & d) 
    IQIndexDat::storage store;
    read(s,store); 
    d.setStore(std::move(store));
long Index::
nblock() const 
    if(not pd) return 0;
    return static_cast<long>(pd->size());
nblock(Index const& i) { return i.nblock(); }
QN const& Index::
qn(long i) const 
    IQINDEX_CHECK_NULL
#ifdef DEBUG
    if(i > nblock())
        Print(nblock());
        Print(i);
        Error("IQIndex::qn arg out of range");
    return pd->qn(i);
qn(Index const& i, long b) { return i.qn(b); }
long Index::
blocksize(long i) const 
    IQINDEX_CHECK_NULL
#ifdef DEBUG
    if(i > nblock())
        Print(nblock());
        Print(i);
        Error("Index::blocksize arg out of range");
    return pd->blocksize(i);
blocksize(Index const& i, long b) { return i.blocksize(b); }
long Index::
blocksize0(long i) const 
    IQINDEX_CHECK_NULL
#ifdef DEBUG
    if(i >= nblock())
        Print(nblock());
        Print(i);
        Error("Index::blocksize0 arg out of range");
    return pd->blocksize0(i);
void Index::
makeStorage(qnstorage && qi)
    pd = std::make_shared<IQIndexDat>(std::move(qi));
totalM(Index::qnstorage const& storage)
    long tm = 0;
    for(auto& iq : storage)
        tm += iq.second;
    return tm;
Index(qnstorage && ind_qn, 
      TagSet const& ts)
  : Index(totalM(ind_qn),ts)
    dir_ = Out;
    makeStorage(std::move(ind_qn));
    if(primeLevel() < 0) setPrime(0);
Index(qnstorage && ind_qn, 
      Arrow dir, 
      TagSet const& ts)
  : Index(totalM(ind_qn),ts)
    dir_ = dir;
    makeStorage(std::move(ind_qn));
    if(primeLevel() < 0) setPrime(0);
QNblock(Index const& I,
        QN const& Q)
    for(auto n : range1(nblock(I)))
        if(qn(I,n) == Q) return n;
    if(not hasQNs(I)) Error("Index does not contain any QN blocks");
    println("I = ",I);
    println("Q = ",Q);
    Error("Index does not contain given QN block.");
    return 0l;
QNblockSize(Index const& I, 
            QN const& Q)
    return blocksize(I,QNblock(I,Q));
void Index::
write(std::ostream& s) const 
    if(!bool(*this)) Error("Index::write: Index is default initialized");
    itensor::write(s,tags_);
    itensor::write(s,id_);
    itensor::write(s,dim_);
    itensor::write(s,dir_);
    if(pd) itensor::write(s,*pd);
    else itensor::write(s,IQIndexDat());
Index& Index::
read(std::istream& s)
    itensor::read(s,tags_);
    if(Global::read32BitIDs())
        using ID32 = std::mt19937::result_type;
        ID32 oldid = 0;
        itensor::read(s,oldid);
        id_ = oldid;
        itensor::read(s,id_);
    itensor::read(s,dim_);
    itensor::read(s,dir_);
    IQIndexDat dat;
    itensor::read(s,dat);
    if(dat.size()>0)
        pd = std::make_shared<IQIndexDat>(std::move(dat.store()));
#ifdef DEBUG
    if(tags_.primeLevel() < 0) Error("Negative primeLevel");
    return *this;
isFermionic(Index const& I)
    for(auto i : range1(I.nblock()))
        for(auto& qnum : I.qn(i).store())
            if(isFermionic(qnum)) return true;
    return false;
#ifdef ITENSOR_USE_HDF5
h5_write(h5::group parent, std::string const& name, Index const& I)
    auto g = parent.create_group(name);
    h5_write_attribute(g,"type","Index",true);
    h5_write_attribute(g,"version",long(1));
    h5_write(g,"id",static_cast<unsigned long>(I.id()));
    h5_write(g,"dim",long(I.dim()));
    h5_write(g,"dir",long(I.dir()));
    h5_write(g,"tags",I.tags());
h5_read(h5::group parent, std::string const& name, Index & I)
    auto g = parent.open_group(name);
    auto type = h5_read_attribute<string>(g,"type");
    if(type != "Index") Error("Group does not contain Index data in HDF5 file");
    auto id = h5_read<unsigned long>(g,"id");
    auto dim = h5_read<long>(g,"dim");
    auto dir = h5_read<long>(g,"dir");
    auto tags = h5_read<TagSet>(g,"tags");
    I = Index(id,dim,toArrow(dir),tags);
} //namespace itensor
Provide feedback

Saved searches

Use saved searches to filter your results more quickly

FilesExpand file tree

index.cc

Latest commit

History

index.cc

File metadata and controls
