tensor.hpp File Reference

A dense, row-major, N-dimensional tensor and Einstein summation. More...

#include <algorithm>
#include <array>
#include <cctype>
#include <cmath>
#include <complex>
#include <cstddef>
#include <deque>
#include <initializer_list>
#include <map>
#include <optional>
#include <stdexcept>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
#include "gemm.hpp"
#include "matrix.hpp"
#include "parallel.hpp"
#include "vector.hpp"

Include dependency graph for tensor.hpp:

Go to the source code of this file.

Classes
class	numerical::einsum_error
	Thrown when an einsum specification is malformed or inconsistent. More...
class	numerical::tensor< T >
	A dense, row-major, N-dimensional tensor of T with value semantics. More...
struct	numerical::detail::contraction_step
	One step of a contraction path: the operands at first and second of the current list, whose result goes to the list's end. More...

Namespaces
namespace	numerical
	The NumericAL public API.
namespace	numerical::detail

Functions
template<typename T >
bool	numerical::approx_equal (const tensor< T > &a, const tensor< T > &b, double tolerance=1e-9)
	Element-wise closeness within a tolerance (for inexact scalars).
template<typename T >
T	numerical::detail::pairwise_sum (const T *first, std::size_t count, std::size_t stride)
	Sum of count values from first at stride stride, by NumPy's pairwise summation: below 8 values in order; up to 128 in eight interleaved partial sums, added pairwise, then the remainder in order; above, the two halves (the first a multiple of 8 long) summed apart.
template<typename T >
void	numerical::detail::pairwise_sum_rows (const T *first, std::size_t count, std::size_t stride, std::size_t width, T *out)
	pairwise_sum of width neighbouring columns at once, into out: column q's sum that of pairwise_sum(first + q, count, stride), to the bit, each step adding a whole row of the columns (contiguous, stride apart) where a column was walked down stride at a time.
template<typename T >
void	numerical::detail::pairwise_sums (const T *data, std::size_t outer, std::size_t count, std::size_t inner, T *out)
	The pairwise sums of outer blocks of count rows of inner elements (row-major), each column of a block summed down its rows, into out (outer × inner): one column at a time where inner is 1, else by rows (pairwise_sum_rows), the same bits.
bool	numerical::detail::is_label (char ch)
	Is ch an ASCII letter (a valid index label)?
std::vector< std::string >	numerical::detail::split_inputs (const std::string &lhs)
	Splits an einsum left-hand side into per-operand subscripts.
template<typename T >
tensor< T >	numerical::detail::tile_gemm (const T *a, const T *b, std::size_t rows, std::size_t inner, std::size_t cols, std::vector< std::size_t > out_shape)
	Matrix product C(M×N) = A(M×K) · B(K×N), row-major.
template<typename T >
tensor< T >	numerical::detail::canonical_product (const T *a, const T *b, std::size_t rows, std::size_t inner, std::size_t cols, std::vector< std::size_t > out_shape)
	The product of the row-major rows × inner a and inner × cols b, shaped out_shape: the canonical two-operand contraction of einsum_impl (a matrix times a vector, a vector times a matrix, a matrix product, a dot product), on any storage that holds the operands contiguously.
template<typename T >
buffer< T >	numerical::detail::vector_times_transposed (const T *a, const T *bt, std::size_t ldbt, std::size_t inner, std::size_t cols)
	canonical_product's vector a (inner) times a matrix given as its transpose bt (cols × inner, its rows ldbt apart, B's columns): each element of the product the sum canonical_product forms, a[k] · B[k][j] added in ascending order of k with add_product's arithmetic (add_one_product), so the same bits; a column of B a run of bt, transposed_columns of them at once.
template<typename T >
std::optional< tensor< T > >	numerical::detail::sum_over_run (const std::string &in, const std::string &output, const std::vector< std::size_t > &shape, const T *data)
	einsum of one row-major operand data (shape, its axes labelled in, each label once) into the labels output, where those it sums are one run of its axes and output the others in order (ij->i, ij->, ij->j, ijk->ik): each output element the sum of evenly spaced elements, by pairwise halving as sum() adds them (pairwise_sums), where einsum's general loop added them one after another. None for any other spec (nothing summed, the summed labels apart, the output reordered).
template<typename T >
buffer< T >	numerical::detail::permuted (const buffer< T > &data, const std::string &from, const std::string &to, const std::map< char, std::size_t > &sizes)
	The elements of data, whose axes carry the labels from (extents in sizes), laid out row-major along the labels to, a permutation of from that differs from it (so at least one label): an odometer over to's outer labels, a strided sweep of its last.
template<typename T >
std::optional< tensor< T > >	numerical::detail::contract_pair (const std::string &sa, const tensor< T > &a, const std::string &sb, const tensor< T > &b, const std::string &output, const std::string &summed, const std::map< char, std::size_t > &sizes)
	A two-operand contraction as matrix products: the labels both operands and the output carry (batch), those one operand and the output carry (free), and those both operands carry alone (summed, in ascending order, the order the general loop sums them in, so each element takes the same terms in the same order) regroup each operand as a stack of matrices, A as batch × free_a × summed and B as batch × summed × free_b, multiplied batch by batch through the packed kernel, batches in parallel; the result is laid out in the output's order; a dot product per batch by dot_product. Nothing for a label repeated within an operand or summed by one operand alone (the general loop takes those).
std::string	numerical::detail::pair_keeps (const std::vector< std::string > &labels, std::size_t i, std::size_t j, const std::string &output)
	The labels of labels[i] and labels[j] that the output or another operand carries, in their order of appearance, once each: what their contraction keeps.
double	numerical::detail::extent_product (const std::string &labels, const std::map< char, std::size_t > &sizes)
	The product of the extents of labels' distinct labels, as a double (a cost: it may pass size_t).
std::vector< std::string >	numerical::detail::after_step (const std::vector< std::string > &labels, contraction_step step, const std::string &output)
	The list after a contraction step: the two operands' labels removed, their result's appended.
void	numerical::detail::cheapest_path (const std::vector< std::string > &labels, const std::string &output, const std::map< char, std::size_t > &sizes, std::vector< contraction_step > &path, double &cost, double &largest)
	The cheapest order of contraction from labels, and its cost and largest intermediate: every order tried, the first found kept among equals.
std::vector< contraction_step >	numerical::detail::contraction_path (std::vector< std::string > labels, const std::string &output, const std::map< char, std::size_t > &sizes)
	The order in which to contract operands labelled labels into output, as opt_einsum chooses it: the cheapest in multiply-adds (the product of the extents a pair's labels span), the smaller largest intermediate among equals, over every order up to exhaustive_path_operands operands; past it, the cheapest pair at each step among those sharing a label (an outer product, cheap in itself, would leave its labels to every later step), the smaller result, then the first pair, among equals. It depends on the shapes alone, so a result is the same on any number of threads.
std::size_t	numerical::detail::ellipsis_in (const std::string &subscripts)
	Evaluates a parsed einsum over its operands.
template<typename T >
tensor< T >	numerical::detail::einsum_impl (const std::string &spec_in, const std::vector< const tensor< T > * > &operands)
template<typename T >
tensor< T >	numerical::detail::einsum_ellipsis (const std::string &spec, const std::vector< const tensor< T > * > &operands)
	einsum_impl of a spec with ellipses, as NumPy reads them: each ... stands for an operand's axes its letters leave, aligned from the right across the operands and broadcast (an extent of 1 repeated), named by letters the spec does not use. An operand's axis of extent 1 under the ellipsis where another's is longer is dropped with its letter, so that operand repeats along it. The output: its ... the broadcast axes (an output without one, where an operand has such axes, is refused, as NumPy refuses it); with no ->, the broadcast axes then the letters that appear once, in order.
template<typename T , typename ... Rest>
tensor< T >	numerical::einsum (const std::string &spec, const tensor< T > &first, const Rest &... rest)
	Einstein summation over one or more tensor operands.
template<typename T >
tensor< T >	numerical::einsum (const std::string &spec, const std::vector< const tensor< T > * > &operands)
	Einstein summation over operands counted at run time: einsum for a caller that holds them in a container (an interpreter evaluating einsum("ij,jk->ik", a, b) from a parse).
template<typename T >
tensor< T >	numerical::matmul (const tensor< T > &a, const tensor< T > &b)
	The matrix product of two tensors, as NumPy's matmul: each operand a stack of matrices along its leading axes, those axes broadcast against the other's (aligned from the right, an extent of 1 repeated), one product per matrix of the broadcast stack through the packed kernel, the products in parallel. A vector first is a row (its axis absent from the result), a vector second a column (likewise); two vectors give their dot product, of rank 0. The broadcast operands are not copied: a repeated matrix is read where it is.

Variables
constexpr std::size_t	numerical::detail::pairwise_columns {512}
	The columns pairwise_sum_rows sums at once: its partial sums of eight rows, 32 KiB of doubles, stay in the L1 cache.
constexpr std::size_t	numerical::detail::vector_chunk {128}
	The columns a thread takes of a vector times a matrix: each of its rows added to them in order.
constexpr std::size_t	numerical::detail::transposed_columns {8}
	The columns of B one task of vector_times_transposed sums together: their sums' chains interleaved, each waiting on its own last addition.
constexpr std::size_t	numerical::detail::exhaustive_path_operands {4}
	The operands at most whose contraction order contraction_path searches exhaustively; more take the cheapest pair at each step.

Detailed Description

A dense, row-major, N-dimensional tensor and Einstein summation.

numerical::tensor is generic over its element type T (the scalar) and over its rank: rank-1 is a vector, rank-2 a matrix, rank-0 a scalar. It is the unifying primitive of the library's vector/tensor algebra. A tensor is built directly from a shape and a row-major buffer, or from a numerical::vector or numerical::matrix (its rank-1 / rank-2 special cases).

numerical::einsum evaluates an Einstein-summation specification such as "ij,jk->ik" (matrix product), "i,i->" (dot product), "ij->ji" (transpose), "ii->" (trace) or "i,j->ij" (outer product), over one or more operands. Each single letter is an index label; a label shared by several axes ties them together; a label absent from the output is summed over. Labels are ASCII letters. Summation runs in a deterministic order (labels sorted ascending), so results are reproducible; a plain dot product ("i,i->") sums in a fixed vectorised order instead, the same on any number of threads (detail::chunked_dot).

Definition in file tensor.hpp.