Multidimensional arrays

TRIQS comes with a library of multidimensional arrays. This library, among others, allows for easy slicing, archiving and algebraic manipulations of multidimensional arrays. Here are a couple of simple examples showing the basic use of this class.

Declaring and printing an array

#include <nda/nda.hpp>
using nda::array;
int main() {
  array<double, 1> A(20);
  std::cout << "A = " << A << std::endl; // arrays are not init by default: this is random
  A() = 2;                               //  assign 2 to a complete view of A.
  std::cout << "A = " << A << std::endl;
  A(4) = 5;
  std::cout << "A = " << A << std::endl;
}

Simple operations

#include <nda/nda.hpp>
using nda::array;
int main() {
  array<double, 1> A(10), B(10);
  A()                = 2;
  B()                = 3;
  array<double, 1> C = A + B;
  std::cout << "C = " << C << std::endl;
}

HDF5 Archiving

Archiving an array into an HDF5 file is easy:

#include <nda/nda.hpp>
#include <nda/h5.hpp>
#include <h5/h5.hpp>

using nda::array;
int main() {

  array<double, 2> A(2, 2);
  A() = 3; // declare and init

  h5::file file("store_A.h5", 'w'); // open the file
  h5::write(file, "A", A);          // write the array as 'A' into the file

  array<double, 2> B; // read the file into B
  h5::read(file, "A", B);
  std::cout << "B = " << B << std::endl;
}

Views: ranges and slices

One can easily take a slice of an array to view and modify only part of the underlying data.

#include <nda/nda.hpp>
using itertools::range;
using nda::array;
using nda::array_view;
int main() {
  array<double, 2> A(3, 3);
  A() = 2.5;
  std::cout << A << std::endl;

  array_view<double, 1> B = A(1, range::all); // select the first line of the matrix
  std::cout << "B = " << B << std::endl;
  B(0) = 1;

  std::cout << "A = " << A << std::endl;
}

Matrices and vectors

Arrays must be distinguished from vectors and matrices, which have an algebra of their own.

#include <nda/nda.hpp>

using nda::array;
using nda::matrix;
using nda::vector;
int main() {
  array<double, 2> A(2, 2), B(2, 2), C;

  A() = 3;
  B() = 1;
  C   = A * B;
  std::cout << "A*B = " << C << std::endl;

  matrix<double> D(2, 2), E(2, 2), F;
  E() = 3;
  E() = 1;
  F   = D * E;
  std::cout << "C*D = " << F << std::endl;

  vector<double> u(2), v(2), w;
  u() = 1;
  v() = 2;
  w   = u + v;

  std::cout << "u+v = " << w << std::endl;
}

Defining through a lazy expression

#include <nda/nda.hpp>
using nda::array;
namespace tql = nda::clef;

int main() {
  tql::placeholder<0> i_;
  tql::placeholder<1> j_;
  array<double, 2> A(2, 2);
  A(i_, j_) << i_ + j_;
  std::cout << "A = " << A << std::endl;
}

Linear algebra

#include <nda/nda.hpp>
#include <nda/linalg.hpp>

using nda::array;
using nda::matrix;
using nda::clef::placeholder;
int main() {
  placeholder<0> i_;
  placeholder<1> j_;
  matrix<double> A(2, 2);
  A(i_, j_) << i_ + j_;

  matrix<double> B = inverse(A);
  double C         = determinant(A);

  std::cout << "A^(-1) = " << B << std::endl;
  std::cout << "det(A) = " << C << std::endl;
}

Map and fold

#include <nda/nda.hpp>
#include <nda/map.hpp>
using nda::array;

double f(int i) { return i * 10; }

int main() {
  auto F = nda::map(std::function<double(int)>(f));
  array<int, 2> A(2, 2);
  A() = 2;

  array<double, 2> B, C;

  A() = 2;
  B   = F(A);
  C   = F(2 * A); // works also with expressions of course

  std::cout << "A = " << A << std::endl;
  std::cout << "F(A) = " << B << std::endl;
  std::cout << "F(2*A) = " << C << std::endl;
}

The full reference of the array library can be found here.