https://github.com/ITensor/ITensor
Tip revision: 7c0a20602668796c706c5911d4203582d6bed9e6 authored by Miles Stoudenmire on 22 October 2013, 11:46:06 UTC
Fixed complex case of tieIndices, fixes #46.
Fixed complex case of tieIndices, fixes #46.
Tip revision: 7c0a206
iqindex.cc
//
// Distributed under the ITensor Library License, Version 1.0.
// (See accompanying LICENSE file.)
//
#include "iqindex.h"
using namespace std;
using boost::format;
using boost::array;
using boost::shared_ptr;
using boost::make_shared;
//
// IQIndexDat
//
class IQIndexDat
{
public:
typedef vector<IndexQN>
StorageT;
typedef StorageT::iterator
iterator;
typedef StorageT::const_iterator
const_iterator;
IQIndexDat() { }
IQIndexDat(const Index& i1, const QN& q1,
const Index& i2 = Index::Null(), const QN& q2 = QN(),
const Index& i3 = Index::Null(), const QN& q3 = QN(),
const Index& i4 = Index::Null(), const QN& q4 = QN(),
const Index& i5 = Index::Null(), const QN& q5 = QN(),
const Index& i6 = Index::Null(), const QN& q6 = QN(),
const Index& i7 = Index::Null(), const QN& q7 = QN(),
const Index& i8 = Index::Null(), const QN& q8 = QN());
explicit
IQIndexDat(vector<IndexQN>& ind_qn);
IQIndexDat(istream& s);
const StorageT&
indices() const { return iq_; }
int
size() { return iq_.size(); }
const Index&
index(int i) { return iq_[i-1]; }
const QN&
qn(int i) { return iq_[i-1].qn; }
iterator
begin() { return iq_.begin(); }
iterator
end() { return iq_.end(); }
const_iterator
begin() const { return iq_.begin(); }
const_iterator
end() const { return iq_.end(); }
void
write(ostream& s) const;
void
read(istream& s);
static const IQIndexDatPtr& Null();
void
makeCopyOf(const IQIndexDat& other);
private:
//////////////////
StorageT iq_;
/////////////////
//Disallow copying using =
void
operator=(const IQIndexDat&);
};
IQIndexDat::
IQIndexDat(const Index& i1, const QN& q1,
const Index& i2, const QN& q2,
const Index& i3, const QN& q3,
const Index& i4, const QN& q4,
const Index& i5, const QN& q5,
const Index& i6, const QN& q6,
const Index& i7, const QN& q7,
const Index& i8, const QN& q8)
{
iq_.push_back(IndexQN(i1,q1));
if(i2 != Index::Null())
iq_.push_back(IndexQN(i2,q2));
if(i3 != Index::Null())
iq_.push_back(IndexQN(i3,q3));
if(i4 != Index::Null())
iq_.push_back(IndexQN(i4,q4));
if(i5 != Index::Null())
iq_.push_back(IndexQN(i5,q5));
if(i6 != Index::Null())
iq_.push_back(IndexQN(i6,q6));
if(i7 != Index::Null())
iq_.push_back(IndexQN(i7,q7));
if(i8 != Index::Null())
iq_.push_back(IndexQN(i8,q8));
}
IQIndexDat::
IQIndexDat(StorageT& ind_qn)
{
iq_.swap(ind_qn);
}
void IQIndexDat::
makeCopyOf(const IQIndexDat& other)
{
iq_ = other.iq_;
}
IQIndexDat::
IQIndexDat(istream& s)
{ read(s); }
void IQIndexDat::
write(ostream& s) const
{
size_t size = iq_.size();
s.write((char*)&size,sizeof(size));
Foreach(const IndexQN& x, iq_)
{
x.write(s);
}
}
void IQIndexDat::
read(istream& s)
{
size_t size; s.read((char*)&size,sizeof(size));
iq_.resize(size);
Foreach(IndexQN& x, iq_)
{
x.read(s);
}
}
const IQIndexDatPtr& IQIndexDat::
Null()
{
static IQIndexDatPtr Null_ = boost::make_shared<IQIndexDat>(Index::Null(),QN());
return Null_;
}
//
// IQIndex Methods
//
#ifdef DEBUG
#define IQINDEX_CHECK_NULL if(pd == 0) Error("IQIndex is null");
#else
#define IQINDEX_CHECK_NULL
#endif
const IQIndexDat::StorageT& IQIndex::
indices() const
{
IQINDEX_CHECK_NULL
return pd->indices();
}
int IQIndex::
nindex() const
{
IQINDEX_CHECK_NULL
return (int) pd->size();
}
const Index& IQIndex::
index(int i) const
{
IQINDEX_CHECK_NULL
#ifdef DEBUG
if(i > nindex())
{
Print(nindex());
Print(i);
Error("IQIndex::index arg out of range");
}
#endif
return pd->index(i);
}
const QN& IQIndex::
qn(int i) const
{
IQINDEX_CHECK_NULL
#ifdef DEBUG
if(i > nindex())
{
Print(nindex());
Print(i);
Error("IQIndex::qn arg out of range");
}
#endif
return pd->qn(i);
}
IQIndex::
IQIndex()
:
dir_(Neither)
{ }
IQIndex::
IQIndex(const string& name,
const Index& i1, const QN& q1,
Arrow dir)
:
Index(name,i1.m(),i1.type(),i1.primeLevel()),
dir_(dir),
pd(boost::make_shared<IQIndexDat>(i1,q1))
{
}
IQIndex::
IQIndex(const string& name,
const Index& i1, const QN& q1,
const Index& i2, const QN& q2,
Arrow dir)
:
Index(name,i1.m()+i2.m(),i1.type(),i1.primeLevel()),
dir_(dir),
pd(boost::make_shared<IQIndexDat>(i1,q1,i2,q2))
{
#ifdef DEBUG
if(i2.type() != i1.type())
{
Print(i2);
Error("Indices must have the same type");
}
#endif
}
IQIndex::
IQIndex(const string& name,
const Index& i1, const QN& q1,
const Index& i2, const QN& q2,
const Index& i3, const QN& q3,
Arrow dir)
:
Index(name,i1.m()+i2.m()+i3.m(),i1.type(),i1.primeLevel()),
dir_(dir),
pd(boost::make_shared<IQIndexDat>(i1,q1,i2,q2,i3,q3))
{
#ifdef DEBUG
if(i2.type() != i1.type()
|| i3.type() != i1.type())
{
Print(i2);
Print(i3);
Error("Indices must have the same type");
}
#endif
}
IQIndex::
IQIndex(const string& name,
const Index& i1, const QN& q1,
const Index& i2, const QN& q2,
const Index& i3, const QN& q3,
const Index& i4, const QN& q4,
Arrow dir)
:
Index(name,i1.m()+i2.m()+i3.m()+i4.m(),i1.type(),i1.primeLevel()),
dir_(dir),
pd(boost::make_shared<IQIndexDat>(i1,q1,i2,q2,i3,q3,i4,q4))
{
#ifdef DEBUG
if(i2.type() != i1.type()
|| i3.type() != i1.type()
|| i4.type() != i1.type())
{
Print(i2);
Print(i3);
Print(i4);
Error("Indices must have the same type");
}
#endif
}
IQIndex::
IQIndex(const string& name,
const Index& i1, const QN& q1,
const Index& i2, const QN& q2,
const Index& i3, const QN& q3,
const Index& i4, const QN& q4,
const Index& i5, const QN& q5,
Arrow dir)
:
Index(name,i1.m()+i2.m()+i3.m()+i4.m()+i5.m(),i1.type(),i1.primeLevel()),
dir_(dir),
pd(new IQIndexDat(i1,q1,i2,q2,i3,q3,i4,q4,i5,q5))
{
#ifdef DEBUG
if(i2.type() != i1.type()
|| i3.type() != i1.type()
|| i4.type() != i1.type()
|| i5.type() != i1.type())
{
Print(i2);
Print(i3);
Print(i4);
Print(i5);
Error("Indices must have the same type");
}
#endif
}
IQIndex::
IQIndex(const string& name,
const Index& i1, const QN& q1,
const Index& i2, const QN& q2,
const Index& i3, const QN& q3,
const Index& i4, const QN& q4,
const Index& i5, const QN& q5,
const Index& i6, const QN& q6,
Arrow dir)
:
Index(name,i1.m()+i2.m()+i3.m()+i4.m()+i5.m()+i6.m(),i1.type(),i1.primeLevel()),
dir_(dir),
pd(new IQIndexDat(i1,q1,i2,q2,i3,q3,i4,q4,i5,q5,i6,q6))
{
#ifdef DEBUG
if(i2.type() != i1.type()
|| i3.type() != i1.type()
|| i4.type() != i1.type()
|| i5.type() != i1.type()
|| i6.type() != i1.type())
{
Print(i2);
Print(i3);
Print(i4);
Print(i5);
Print(i6);
Error("Indices must have the same type");
}
#endif
}
int
totalM(const IQIndexDat::StorageT& storage)
{
int tm = 0;
Foreach(const IndexQN& iq, storage)
{
tm += iq.m();
#ifdef DEBUG
if(iq.type() != storage.front().type())
Error("Indices must have the same type");
#endif
}
return tm;
}
IQIndex::
IQIndex(const string& name,
IQIndexDat::StorageT& ind_qn,
Arrow dir, int plev)
:
Index(name,totalM(ind_qn),ind_qn.front().type(),plev),
dir_(dir),
pd(new IQIndexDat(ind_qn))
{
}
IQIndex::
IQIndex(const Index& index, const IQIndexDatPtr& pdat)
:
Index(index),
dir_(In),
pd(pdat)
{ }
void IQIndex::
conj() { dir_ = -dir_; }
void IQIndex::
write(ostream& s) const
{
IQINDEX_CHECK_NULL
Index::write(s);
s.write((char*)&dir_,sizeof(dir_));
pd->write(s);
}
void IQIndex::
read(istream& s)
{
Index::read(s);
s.read((char*)&dir_,sizeof(dir_));
pd = boost::make_shared<IQIndexDat>();
pd->read(s);
}
string
showm(const IQIndex& I)
{
#ifdef DEBUG
if(I.isNull()) Error("Null IQIndex");
#endif
string res = " ";
ostringstream oh;
oh << I.m() << " | ";
Foreach(const IndexQN& iq, I.indices())
{
oh << iq.qn << ":" << iq.m() << " ";
}
return oh.str();
}
void IQIndex::
primeLevel(int val)
{
solo();
Index::primeLevel(val);
Foreach(IndexQN& iq, *pd)
iq.primeLevel(val);
}
void IQIndex::
prime(int inc)
{
solo();
Index::prime(inc);
Foreach(IndexQN& iq, *pd)
iq.prime(inc);
}
void IQIndex::
prime(IndexType type, int inc)
{
solo();
Index::prime(type,inc);
Foreach(IndexQN& iq, *pd)
iq.prime(type,inc);
}
void IQIndex::
mapprime(int plevold, int plevnew, IndexType type)
{
solo();
Index::mapprime(plevold,plevnew,type);
Foreach(IndexQN& iq, *pd)
iq.mapprime(plevold,plevnew,type);
}
void IQIndex::
noprime(IndexType type)
{
solo();
Index::noprime(type);
Foreach(IndexQN& iq, *pd)
iq.noprime(type);
}
void IQIndex::
solo()
{
IQINDEX_CHECK_NULL
if(!pd.unique())
{
const IQIndexDat& olddat = *pd;
pd = boost::make_shared<IQIndexDat>();
pd->makeCopyOf(olddat);
}
}
const IQIndex& IQIndex::
Null()
{
static const IQIndex Null_(Index::Null(),IQIndexDat::Null());
return Null_;
}
void
calc_ind_ii(const IQIndexVal& iv, int& j, int& ii)
{
j = 1;
ii = iv.i;
while(ii > iv.index(j).m())
{
ii -= iv.index(j).m();
++j;
}
}
IQIndexVal::
IQIndexVal()
:
IQIndex(IQIndex::Null()),
i(0)
{ }
IQIndexVal::
IQIndexVal(const IQIndex& iqindex, int i_)
:
IQIndex(iqindex),
i(i_)
{
#ifdef DEBUG
if(i > m() || i < 1)
{
Print(iqindex);
Print(i);
Error("IQIndexVal: i out of range");
}
#endif
}
IndexQN IQIndexVal::
indexqn() const
{
int j,ii;
calc_ind_ii(*this,j,ii);
return IndexQN(IQIndex::index(j),IQIndex::qn(j));
}
const QN& IQIndexVal::
qn() const
{
int j,ii;
calc_ind_ii(*this,j,ii);
return IQIndex::qn(j);
}
bool IQIndexVal::
operator==(const IQIndexVal& other) const
{
return (IQIndex::operator==(other) && i == other.i);
}
IQIndexVal::
operator IndexVal() const
{
return IndexVal(Index(*this),i);
}
IndexVal IQIndexVal::
blockIndexVal() const
{
if(*this == IQIndexVal::Null())
return IndexVal::Null();
int j,ii;
calc_ind_ii(*this,j,ii);
return IndexVal(this->index(j),ii);
}
/*
IQIndexVal::
operator ITensor() const
{
return ITensor(IndexVal(iqind,i));
}
*/
IQIndexVal IQIndex::
operator()(int n) const
{
return IQIndexVal(*this,n);
}
bool
hasindex(const IQIndex& J, const Index& i)
{
Foreach(const Index& j, J.indices())
{
if(j == i) return true;
}
return false;
}
int
offset(const IQIndex& I, const Index& i)
{
int os = 0;
Foreach(const IndexQN& iq, I.indices())
{
if(iq == i) return os;
os += iq.m();
}
Print(I);
Print(i);
Error("Index not contained in IQIndex");
return 0;
}
QN
qn(const IQIndex& I, const Index& i)
{
Foreach(const IndexQN& jq, I.indices())
{
if(jq == i)
return jq.qn;
}
cout << I << "\n";
cout << "i = " << i << endl;
Error("IQIndex does not contain given index.");
return QN();
}
Index
findByQN(const IQIndex& I, const QN& qn)
{
Foreach(const IndexQN& jq, I.indices())
{
if(jq.qn == qn)
return jq;
}
cout << I << "\n";
cout << "qn = " << qn << endl;
Error("IQIndex does not contain given QN block.");
return Index();
}
ostream&
operator<<(ostream &o, const IQIndex& I)
{
if(I.isNull())
{
o << "IQIndex: (null)";
return o;
}
o << "IQIndex: " << Index(I) << " <" << I.dir() << ">" << endl;
for(int j = 1; j <= I.nindex(); ++j)
o << " " << I.index(j) SP I.qn(j) << "\n";
return o;
}
std::ostream&
operator<<(std::ostream &s, const IndexQN& x)
{
const Index& i = x;
return s << "IndexQN: " << i
<< " (" << x.qn << ")\n";
}
std::ostream&
operator<<(std::ostream& s, const IQIndexVal& iv)
{
const IQIndex& I = iv;
return s << "IQIndexVal: i = " << iv.i << ", " << I << "\n";
}