Currently IQCombiner is allowed to combine IQIndices with different arrow directions. This may not actually make sense fundamentally and is probably pointing to a more correct algorithm. As a specific example, we should probably precede orthogonalization of MPOs (where we treat them as MPS's) by some operation that flips their bra arrows.

It should be easier to inherit from the Sweeps class to allow library users to implement their own customizations. Sweeps constructor can call a virtual method that sets sweep parameters etc.

Because of the way ITensors order their indices and access their data, reading in a chain of numbers directly into the ITensor's dat Vector makes commaInit the transpose of what a user would likely expect (the reverse of conventional matrix notation).

Page 6 of Singh and Vidal, arxiv:1307.1522 may be helpful.

ITensor::operator*= ends with a call to scaleOutNorm which ensures that any overall very large or small scale gets pulled into the LogNumber. For IQTensors it makes sense for all constituent ITensors to have the same LogNumber since they represent IQTensor blocks.

I found an extra PA in the current version of ./matrix/lapack_wrap.h
someone just saved the file with a typo.
Once eliminated everything has compiled fine on my linux box.
Hope this helps

Unlike in CMPS class, may be more appropriate to have MPO type objects handle projection of themselves. Return type can be typedef'd within MPO type object.

HermitianEigenvalues in utility.cc uses variable length arrays. They need to be replaced by std::vector. However, using static variables in this function and also functions in lapack_wrap.h breaks thread-safety. Is this not an issue?

In the older version of the code, even if an IQTensor had a non-trivial Virtual IQIndex (viqindex), its constituent ITensors did not. This was to prevent excessive, possibly inefficient, manipulation of Indices. However, now that m==1 Index objects are stored in a very lightweight fashion within an ITensor this behavior should be changed back.

Currently I[Q]Tensor::operator*= and operator/= do not check if this == &other and give an incorrect result if this is the case.

Currently accessing elements of ITensors (IQTensors) uses an overload of operator() taking IndexVal (IQIndexVal) arguments.

This has at least two major design flaws:

• Unless the tensor variable is const, the compiler will always choose the non-const version of the overload even for read-only access. This can lead to extra work being done and unwanted side-effects.
• It's not possible to overload on return values, so these methods can't be used to access elements of complex tensors as scalars of type Complex.

Of course it has the major advantage of readability and symmetry between getting and setting elements.

Possible fixes:

• Different names (other than operator()) for either the get or set methods, or both.
• operator() could return a temporary object holding the particular element. This object could have overloaded constructors and support automatic conversion to both Real and Complex, giving a warning/error if assigning to Real would drop a non-zero imaginary part. Of course this approach means added complexity and possibly extra confusion to users e.g. who store the returned object to an auto variable yet expect it to be a Real.

The parameters controlling truncation, namely minm_, maxm_ and cut_ within SVDWorker only have indirect control on the final number of states kept through the variable docut. Also the code for computing docut is pretty compliacated. It would be good to simplify this code for readability and maintainability, and to write some tests to ensure that guarantees about minm_ etc. are being enforced.

Make sure to test that the truncate_ flag works correctly (when set false) as well.

For simulations generating a huge number of indices (such as continuum DMRG), collisions are occurring between the total unique reals of groups of indices.

Here is a sample output showing two IndexSets with the same unique real but clearly different indices:

iset =
0 qlink(Link,5420):9 5.4205104068E-01
1 qlink(Link,1662):11 1.6627267101E-01
2 UpDn for site 318(Site,1540):1 1.5402714266E-01
3 Emp for site 319(Site,3385):1 3.3851180119E-01

oset =
0 qlink(Link,2513):33 2.5135507685E-01
1 qlink(Link,7630):8 7.6306582481E-01
2 Dn for site 318(Site,1274):1 1.2744677526E-01
3 UpDn for site 319(Site,589):1 5.8994978620E-02

iset uniqueReal = 1.2008626555E+00
oset uniqueReal = 1.2008626555E+00
uniqueReal diff = 1.4210854715E-13

Currently *= is just the original code for * followed by a self assignment. Much further optimization should be possible.

Currently only non-const versions of ITensor::operator()(IndexVal, ...). It should be possible to use these methods on const ITensor objects for reading out individual components, however.

Currently the IQTensor to ITensor converter expands each block in pieces using expandIndex. This is inefficient because all of the additional zeros outside the block are being allocated over and over. Instead, make the final ITensor with the correct size, loop over the IQTensor blocks, and add their elements to the ITensor's with the correct offsets (similar to expandIndex).

Currently SVDWorker holds the svd, csvd, and denmatDecomp methods, as well as being a struct for setting SVD accuracy parameters and saving results like the untruncated singular values/eigenvalues. A better design could be to have svd, etc. as free methods and pass in a struct with the accuracy params. On return this struct could hold the eigs kept. Instead of setting the showeigs flag, one could just print the resulting struct.

Hello,

More of a question than a issue. Is there a way of removing a site from an MPS/MPO that is more efficient than a simple loop? The loop below copies the MPO tensor to the right across starting at some position (pos) that is to be removed:

for(int i=pos; i<N; ++i)
H.Anc(i) = H.Anc(i+1);

Many thanks

Imaginary part of complex IQTensors not ending up in the right components of the resulting ITensor in the current implementation.

ITensor 2.0 is very exciting! I cloned it immediately after receiving the email notification. It compiles fine with g++ and PLATFORM=macos, but when I tried icpc with PLATFORM=mkl or macos, the compilation process would encounter the same error message (full output attached). Checking the compiler version with "icpc -v" shows "icpc version 16.0.2 (gcc version 4.9.0 compatibility)". Does someone know what's causing the error? Thanks in advance!
error.txt

Although statements like using std::cout help keep the library code less verbose, they can cause problems and conflicts within client code.

Currently we have an SVD method and a density matrix decomposition method in svdworker.h. But neither can be used for doing eigenvalue decompositions even though the underlying code could support this.

In VectorRef & VectorRef::operator *= (Real a) defined in matrixref.cc, is there a reason you declare dscal inside the method, or rather, why you define your own version of dscal in the first place? It causes linkage problem in MSVC, for instance when I try to compile tutorial one.cc.

matrix.lib(matrixref.obj) : error LNK2019: unresolved external symbol "void __cdecl dscal(int,double,double *,int)" (?dscal@@YAXHNPEANH@Z) referenced in function "public: class itensor::VectorRef & __cdecl itensor::VectorRef::operator*=(double)" (??XVectorRef@itensor@@QEAAAEAV01@N@Z)

A quick fix would be to move the declaration outside, perhaps to line 62 together with copyscale and dcopy, or to use blas's dscal_.

The EigenValues method in matrix utility.cc does Z = Z.t() after computing the eigenvectors and storing them into Z. Can we change our conventions in methods calling EigenValues to avoid this copy?

I think there is a scale_.real()* missing in the case that the tensor is diagonal.

Real ITensor::
operator()(const IndexVal& iv1, const IndexVal& iv2) const
    {
    ITENSOR_CHECK_NULL
    if(type_ == Diag)
        {
        if(iv1.i != iv2.i) return 0;
        return r_->v.at(iv1.i-1);
        }
    return scale_.real()*r_->v[_ind2(iv1,iv2)];
    }

Currently conj(...) is Hermitian conjugate. Would be useful to have just complex conjugate, but conj would be an ideal name so this could be a breaking change.

Possible names:

• Non-breaking (Hermitian conj remains just "conj"), complex conj becomes:
  • cconj(...)
  • star(...)
  • cc(...)
• Breaking (complex conj takes over the name "conj"), Hermitian conj becomes:
  • dagger(...)
  • dag(...)
  • hc(...)
• Breaking: rename both, something like cc(...) and hc(...) or cc(...) and dag(...).

My original goal was to see how the LocalOp::diag() function works (defined in itensor/mps/localop.h). Before ITensor 2, it depended on a function called tieIndices(). In ITensor 2, it involves noprime() and diagTensor() functions to tie the indices, and I want to see if it works the same way.

So, in the itensor/eigensolver.h file, to the davidson() function I added a line of
Tensor ddd = A.diag();
which would call the diag() of a LocalMPO (or LocalMPOSet) object and then wrap the function call to the diag() of a LocalOp object. And then I started a usual DMRG run. However, this resulted in
error: no matching function for call to 'diagTensor'
at line 326:31 of itensor/mps/localop.h
auto Diag = nonprime(Op1 * diagTensor(1,toTie,prime(toTie)),toTie);
The complete error message is attached. I think the template argument deduction/substitution failed for the diagTensor() function.

As for my configuration, I used g++ and framework Accelerate to compile ITensor 2 (the default setting of options.mk). And 'g++ --version' shows
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.3.0 (clang-703.0.29)
Target: x86_64-apple-darwin15.4.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

Does anybody know how to fix this? Thank you!
error.txt

Currently, all of the IndexDat read/write is handled in Index. Moving it inside IndexDat should make the encapsulation better and simplify the entire process (e.g. no special case for ReIm type).

Currently the Sweeps class fails silently when using table mode with too few lines of sweeps input.

Currently copying an IQTensor involves copying two large vectors containing many ITensors, Index's etc. It should be possible to pool IQTensor data to save both memory and time spent making unnecessary copies.

I was running lots of DMRG simulations in parallel (using MPI), and one of them just got crashed with following message:

j = 0
iset =
0 qlink(Link,8594):14 | 8.5940636039E-01
1 Up for site67(Site,2368):1 | 2.3686394567E-01
2 Up for site68(Site,3031):1 | 3.0319490407E-01
3 qlink(Link,2891):1 | 2.8917708278E-01

oset = 
0 qlink(Link,2740):5 | 2.7402408078E-01
1 qlink(Link,6117):19 | 6.1170528075E-01
2 Dn for site67(Site,5288):1 | 5.2881146305E-01
3 Dn for site68(Site,2741):1 | 2.7410146833E-01

iset uniqueReal = 1.688642292909474E+00
oset uniqueReal = 1.688642292909474E+00
uniqueReal diff = 0.000000000000000E+00
terminate called after throwing an instance of 'ITError'
[eol-laptop:14374] *** Process received signal ***
[eol-laptop:14374] Signal: Aborted (6)
[eol-laptop:14374] Signal code:  (-6)
[eol-laptop:14374] [ 0] /lib64/libc.so.6(+0x33aa0)[0x7fca13aa7aa0]
[eol-laptop:14374] [ 1] /lib64/libc.so.6(gsignal+0x37)[0x7fca13aa7a27]
[eol-laptop:14374] [ 2] /lib64/libc.so.6(abort+0x16a)[0x7fca13aa8dba]
[eol-laptop:14374] [ 3] /usr/lib/gcc/x86_64-pc-linux-gnu/4.9.2/libstdc++.so.6(_ZN9__gnu_cxx27__verbose_terminate_handlerEv+0x16d)[0x7fca143840bd]
[eol-laptop:14374] [ 4] /usr/lib/gcc/x86_64-pc-linux-gnu/4.9.2/libstdc++.so.6(+0x5ef06)[0x7fca14381f06]
[eol-laptop:14374] [ 5] /usr/lib/gcc/x86_64-pc-linux-gnu/4.9.2/libstdc++.so.6(+0x5ef51)[0x7fca14381f51]
[eol-laptop:14374] [ 6] /usr/lib/gcc/x86_64-pc-linux-gnu/4.9.2/libstdc++.so.6(+0x5f168)[0x7fca14382168]
[eol-laptop:14374] [ 7] ./dmrg(_ZN7itensor7getpermINS_5IndexEEEvRKNS_8IndexSetIT_EERKNS4_7StorageERNS_11PermutationE+0x3d3)[0x4dd253]
[eol-laptop:14374] [ 8] ./dmrg(_ZN7itensor7ITensorpLERKS0_+0xb65)[0x4d34a5]
[eol-laptop:14374] [ 9] ./dmrg(_ZN7itensor8IQTensormLERKS0_+0x931)[0x4eb641]
[eol-laptop:14374] [10] ./dmrg(_ZNK7itensor7LocalOpINS_8IQTensorEE4diagEv+0x110d)[0x4b17a7]
[eol-laptop:14374] [11] ./dmrg(_ZNK7itensor8LocalMPOINS_8IQTensorEE4diagEv+0x27)[0x4ad8c9]
[eol-laptop:14374] [12] ./dmrg(_ZNK7itensor12LocalMPO_MPSINS_8IQTensorEE4diagEv+0x27)[0x4aaad9]
[eol-laptop:14374] [13] ./dmrg(_ZN7itensor15complexDavidsonINS_12LocalMPO_MPSINS_8IQTensorEEES2_EESt6vectorISt7complexIdESaIS6_EERKT_RS4_IT0_SaISC_EERKNS_6OptSetE+0x6b7)[0x4a3782]
[eol-laptop:14374] [14] ./dmrg(_ZN7itensor8davidsonINS_12LocalMPO_MPSINS_8IQTensorEEES2_EESt6vectorIdSaIdEERKT_RS4_IT0_SaISA_EERKNS_6OptSetE+0x8c)[0x49cb4f]
[eol-laptop:14374] [15] ./dmrg(_ZN7itensor8davidsonINS_12LocalMPO_MPSINS_8IQTensorEEES2_EEdRKT_RT0_RKNS_6OptSetE+0xac)[0x4970c2]
[eol-laptop:14374] [16] ./dmrg(_ZN7itensor10DMRGWorkerINS_8IQTensorENS_12LocalMPO_MPSIS1_EEEEdRNS_4MPStIT_EERT0_RKNS_6SweepsERNS_12DMRGObserverIS5_EENS_6OptSetE+0x5ce)[0x4926ce]
[eol-laptop:14374] [17] ./dmrg(_ZN7itensor10DMRGWorkerINS_8IQTensorENS_12LocalMPO_MPSIS1_EEEEdRNS_4MPStIT_EERT0_RKNS_6SweepsERKNS_6OptSetE+0x9e)[0x48c33b]
[eol-laptop:14374] [18] ./dmrg(_ZN7itensor4dmrgINS_8IQTensorEEEdRNS_4MPStIT_EERKNS_4MPOtIS3_EERKSt6vectorIS4_SaIS4_EERKNS_6SweepsERKNS_6OptSetE+0x86)[0x487872]
[eol-laptop:14374] [19] ./dmrg(_Z14run_simulationN4alps6paramsEiSs+0x675)[0x4770db]
[eol-laptop:14374] [20] ./dmrg(_Z5slaveN5boost3mpi12communicatorEN4alps6paramsESs+0x170)[0x477838]
[eol-laptop:14374] [21] ./dmrg(main+0x6db)[0x477f70]
[eol-laptop:14374] [22] /lib64/libc.so.6(__libc_start_main+0xf0)[0x7fca13a93fa0]
[eol-laptop:14374] [23] ./dmrg[0x47699f]
[eol-laptop:14374] *** End of error message ***

To be specific, I was using modified version of Heisenberg hamiltonian (with added random magnetic field), and crash occured somewhere in the middle (so it seem to be hardly reproducible)

Davidson converges much faster in some cases than Lanczos. You get Lanczos if you
set the diagonal part of H to 1 in Davidson's algorithm. That is what we do. It will be
fast and easy to put in the actual diagonal part of H.

It may be possible, and a much better design, to make tensors and sparse tensors the same class/type. They both have nearly the same external interface and behave similarly i.e. have the same purpose. Internally, they both store their data as a Vector of Reals. The only difference is how the data in this Vector is interpreted and used. Therefore, we could simply add a flag that switches a tensor from being sparse to being dense.

For example, a tensor default initialized to zero or initialized by hand (such as setting a few values explicitly or via commaInit) would be sparse (this includes site operators and the tensors in an MPO). Most sparse-sparse products would again be sparse. Sparse-dense products, the results of a Davidson calculation, or the unitaries created by an SVD would be dense.

This could have several advantages, such as leveraging sparsity more throughout the library and making sparse tensors easier to use (e.g. novice users wouldn't have to know about them to use the svd method). Also, it would guarantee that sparse and dense tensors have the same interface and reduce code maintainence and documentation.

The pairTie function uses Kronecker delta type ITensors internally. This means that contracting product loops over many elements that are just set to zero. Replacing these ITensors with explicit sums over the Indices to be tied together should be much more efficient.

Note that an IQTensor with all m=1 sub-indices is nearly identical in content to a SiteOp object. Can think of such an IQTensor as a lazy ITensor constructor. Have model objects return IQTensors instead of SiteOps.

Providing DMRG as a function has many limitations. For one example the argument list is rather rigid and not self-documenting. Reimplementing DMRG as a class makes it easier to use and much more customizable.

Currently SiteSet accepts operator strings of the form "Op1 * Op2". Extend this to "Op1 * Op2 * Op3 * ..."

global.h does not compile in MSVC, because of its dependence on unistd.h and mkdtemp. How about making a change to global.h as the following?

#if defined(_MSC_VER)
#include <process.h>
#include <io.h>
#include <direct.h>

inline char* mkdtemp(char *templat)
    {
    char* retval = _mktemp(templat);
    if (retval) { _mkdir(retval); }
    return retval;
    }
#else
#include <unistd.h>
#endif

Even when using the same minm=maxm=m of the bond between the two tensors fed into SVDWorker::operator(), certain MPOs are getting overly truncated. Probably this is because all of the eigs are put into alleig and sorted, which loses information about which block they came from.

Currently primelevel methods accept optional PrimeType distinct from IndexType to allow e.g. primeBoth. However, simpler & more readable to just add new IndexType All and directly specify the IndexType to be primed. (Using IndexType All in Index constructor is forbidden so improper usage this way not an issue.)

The function putInQNs does imaginary time evolution to ensure that an IQTensor has all possible QN sectors. This is so it can be written to a Vector (using AssignToVec), modified, then re-read (using AssignFromVec). It should be possible to eliminate putInQNs altogether.

Possible approaches include putting in missing blocks explicitly, or writing a custom Davidson code for ITensors and IQTensors, which we may want to do anyway.

I have tried using the HubbardChain Hamiltonian in a simple test code based on the exthubbard.cc sample. It compiles fine but when running I get the error:

"Operator F*Cup name not recognized"

The code is included at the bottom. Changing "HubbardChain" to "ExtendedHubbard" makes everything run fine. As I am just learning to use ITensor I admit that it may well be a user error.

Many thanks

include "core.h"

include "model/spinhalf.h"

include "hams/HubbardChain.h"

using boost::format;
using namespace std;

int main(int argc, char* argv[])
{
//
// Initialize the site degrees of freedom.
//
int N = 50;
int Npart = N; //number of particles, default is N (half filling)

Hubbard model(N);

//
// Create the Hamiltonian matrix product operator.
// Here we use the IQMPO class which is an MPO of 
// IQTensors, tensors whose indices are sorted
// with respect to quantum numbers
//
IQMPO H = HubbardChain(model);

//
// Set the initial wavefunction matrix product state
// to be a Neel state.
//
InitState initState(model);
int p = Npart;
for(int i = N; i >= 1; --i) 
    {
    if(p > i)
        {
        cout << "Doubly occupying site " << i << endl;
        initState.set(i,&Hubbard::UpDn);
        p -= 2;
        }
    else
    if(p > 0)
        {
        cout << "Singly occupying site " << i << endl;
        initState.set(i,(i%2==1 ? &Hubbard::Up : &Hubbard::Dn));
        p -= 1;
        }
    else
        {
        initState.set(i,&Hubbard::Emp);
        }
    }

IQMPS psi(model,initState);

cout << totalQN(psi) << endl;

//
// Set the parameters controlling the accuracy of the DMRG
// calculation for each DMRG sweep. 
// Here less than 5 cutoff values are provided, for example,
// so all remaining sweeps will use the last one given (= 1E-10).
//
Sweeps sweeps(5);
sweeps.maxm() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
sweeps.niter() = 2;
sweeps.noise() = 1E-7,1E-8,0.0;
cout << sweeps;

Real En = dmrg(psi,H,sweeps,Quiet());

//
// Print the final energy reported by DMRG
//
cout << format("\nGround State Energy = %.10f\n")%En;
cout << format("\nUsing psiHphi = %.10f\n") % psiHphi(psi,H,psi);

cout << "\nTotal QN of Ground State = " << totalQN(psi) << "\n";

return 0;

}                           

When arrows are involved, properties of indices can be asymmetric between a pair of tensors. The convention of index_in_common is to take the index from the first tensor. Making it a class method would make its behavior more obvious: A.indexInCommon(B); (the index comes from A).