#5 Try out numerical libraries

Before we start with the serious coding, we need to decide on the mathematics libraries to use. Two components are needed:

1. A library for matrix handling. Performance has a lower priority than ease of use. The library needs to be able to:

   • Represent vectors and matrices
   • Allow the basic maths: Scalar products between vectors, elementwise multiplications, and matrix multiplications
   • Represent higher-dimensional objects, such as a "vector" A_{ijk} with i,j,k being indices.
   • Allow easy manipulations and e.g. dimension conversions, something like Matlab or python's Numlib.

2. A library for doing the Fast Fourier Transform

   • It should be preferably easy to use.
   • It must be in some way compatible, or easily be made compatible, with the objects appearing the matrix handling library.

Your task: Look for some libraries, read the documentation/play around with them, and post a summary to the mailing list. Preferably, assemble a small Hello World project that demonstrates the use of the libraries.

[Anonymous]

After just looking at freecode and filtering out projects without releases for years, I got for the matrix libraries:

Eigen
IT++
Seldon
Aramdillo

[Anonymous]

All the matrix libraries seem to deal only with matrices, which is bad; we would like to operate with tensors and would have to wrap the matrices up in some way, which is complicated.

TODO: another search for tensor libraries.

[Anonymous]

After another search for tensor libraries, I ended up with the following candidates:

Tensor C++ library; looks like a smaller library
ITensor; looks cute, not too big, and pretty much exactly what I was thinking about (also used for electronic structure problems), no idea how good an FFT conncetion could be done.
Portable science toolkit; looks like overkill, but maybe it is nice?
GSL tensor extensions; GSL is a pretty standard library, so maybe this is worth a try.

[Anonymous]

Sort summary after a bit more digging:

1. GSL looks bad. No documentation, clunky C interface and no releases for two years.

2. Portable science toolkit looks a bit like overkill. More on the side of massively parallel solvers that can also be run on a single machine. Many dependencies and such. On the other hand, there is an interface to MPI and GPU calculations... Will try it out if time and interest remains.

3. ITensor and Tensor C++ library look good. Clean interface, easy usage, even unit testing.

Now I need to develop some test cases...

[Anonymous]

Suggestion for the test case:

1. A two-dimensional problem, so that the wave function is a tensor of rank 2. Let us call the degrees of freedom "x" and "y".

2. Create some testing wave function X(x,y); random numbers are fine.

3. Typical use-case: Apply an operator A that acts only along x. Create a 2D tensor A(x,x') that corresponds to the operator, and apply it in a generic way to the wave function X'{ij} = sum_k A{ik} X_{kj}

4. Typical usage: Apply a potential energy operator B (element-wise multiplication). Create a tensor B(x,y) and create X'{ij} = B{ij} * X_{ij}

5. Typical usage: Calculate an expectation value. Needs a 1-D array x, assume a homogeneous weight h along y and a weight array w along x, and calculate <x> = sum_{ij} conj(X_{ij}) * x_i * X_{ij} * w_i * h.

6. Typical usage: Standard functions (exponentials). Take a tensor T and calculate exp(T), elementwise multiplication of two tensors, slicing (setting/getting a slice), tensor product (get a 2D-tensor as the dyadic product of two 1-D tensors).

7. Check connection to FFT algorithms.

8. Calculate eigenvalues of a Hamiltonian. Set up a Hamiltonian H(x,x') (use 1-D problem for simplicity), and calculate eigenvalues; how is connection to LAPACK?

9. Try out sparse tensors.

[Anonymous]

Issues with the Tensor C++ library:

1. The checked out sources (at least the git head) were broken in several places. What took me longest was an error where an include file had another function declaration than the corresponding code file, leading to a library with an unresolved dependency.
2. The documentation is severely lacking. There is a good tutorial, and you can pretty much look up things in the code, but many functions are simply not documented. So you have to look into the code to figure out what, e.g., foldc does.
3. Many functions are simply awkward to use. Not at last because they mimic Matlab's functions. For example, getting a Row vector of 10 ones is not done via RTensor::ones(10), because that produces (Matlab style) a 10x10 matrix. Also there are many cases just lacking. For example, I wanted to foldin() a CTensor with an RTensor. This does not work; you either need two CTensors or two RTensors, mixed cases are not considered.
4. There are some assumptions hidden in the code. For example, it is always assumed that the second dimension has more points than the first one. Either that or the code is bugged.

Edit: A new checkout compiled and ran fine; no idea what was wrong. The last bug was in my testing script (never use nx and ny, they look too similar ...)

[Anonymous]

Next: Itensor.

Issues encountered so far:

1. The library is surprisingly huge; sloccount gives me 35.000 lines of code instead of the 3.500 for Tensor C++. Not quite sure where they all come from.
2. Part of the first issue: The tensor code is intertwined with DMRG routines (electronic structure method); not sure how and if this can be disentangled, but probably some work.
3. Installation is a pain. The library requires either Intel MathKernel libraries (expensive), AMD's AMCL libraries (free, but restricted license, also non-standard BLAS interfaces), or a Mac, which is rather strange.

To be continued...

[Anonymous]

Since I could not easily install ITensor, I decided to instead have a look at the code, especially at the tests........Do all C++ developers have to write such ugly code?

Looking at some headers, I can judge the following:

1. Inconsistent interface. Assume 2 Tensors A,B, and a real number r. Then the tensor contraction is given by AB, the tensor with all elements multiplied by a real number by rA, tensor with all elements divided by a real number A/r, and the elementwise product of two tensors by A/B (!!!!). I understand the rationale, but this is pretty unintuitive
2. Some functions I do not like. For example, the rank of a tensor is given by A.r(). Not self-explaining at all.
3. A lot of bloat. There is for example the multiplication with a LogNumber (very small or very large number). Normally you would just say "don't do this.". As a consequence, the definition of the tensor itself is already > 400 Lines of Code.

[Anonymous]

After finally compiling the library on a cluster at work and trying it out, two more things:

1. Either I am missing the obvious or a simple matrix multiplication gives me completely wrong results. Might be a development version or whatever, but this is pretty bad.
2. The tensors are surprisingly awkward to use. The problem is that you do not access a tensor via A(i,j,k), but always have to supply and memorize indices, i.e., A(index1(i), index2(j), index3(k)). I found no easy way to set ranges (multiple values at once); there seems to be a complicated one that involves conversion to a vector and setting a range there or so.

To say also something good: There are a couple of interesting functions, for example tieing indices together (taking a 3D-tensor A_ijk and making a 2D-Tensor A_xk from it; that is something that I was missing in Matlab once).

Final verdict: The library is probably ok for Density Matrix Renormalization Groups, where you create high-dimensional tensors with few entries in each direction (spins, i.e., 2 in the examples) and want to contract them here and there. But I am not quite sure how far it matches our use-case.

[Anonymous]

Selected Tensor C++ library.

This project has moved to github now:
https://github.com/juanjosegarciaripoll/tensor

[Anonymous]

The FFT part will also be done later. In principle, you can do a FT with a matrix multiplication; this is not fast, but it works. So the FFT part should be deferred in favor of getting ahead.