README.md

Each coefficient in the result is equal to the corresponding coefficient in the
'then' tensor if the corresponding value in the 'if' tensor is true. If not, the
resulting coefficient will come from the 'else' tensor.


## Contraction

Tensor *contractions* are a generalization of the matrix product to the
multidimensional case.

    // Create 2 matrices using tensors of rank 2
    Eigen::Tensor<int, 2> a(2, 3);
    a.setValues({{1, 2, 3}, {6, 5, 4}});
    Eigen::Tensor<int, 2> b(3, 2);
    a.setValues({{1, 2}, {4, 5}, {5, 6}});

    // Compute the traditional matrix product
    array<IndexPair<int>, 1> product_dims = { IndexPair(1, 0) };
    Eigen::Tensor<int, 2> AB = a.contract(b, product_dims);

    // Compute the product of the transpose of the matrices
    array<IndexPair<int>, 1> transpose_product_dims = { IndexPair(0, 1) };
    Eigen::Tensor<int, 2> AtBt = a.contract(b, transposed_product_dims);


## Reduction Operations

A *Reduction* operation returns a tensor with fewer dimensions than the
original tensor.  The values in the returned tensor are computed by applying a
*reduction operator* to slices of values from the original tensor.  You specify
the dimensions along which the slices are made.

The Eigen Tensor library provides a set of predefined reduction operators such
as ```maximum()``` and ```sum()``` and lets you define additional operators by
implementing a few methods from a reductor template.

### Reduction Dimensions

All reduction operations take a single parameter of type
```<TensorType>::Dimensions``` which can always be specified as an array of
ints.  These are called the "reduction dimensions."  The values are the indices
of the dimensions of the input tensor over which the reduction is done.  The
parameter can have at most as many element as the rank of the input tensor;
each element must be less than the tensor rank, as it indicates one of the
dimensions to reduce.

Each dimension of the input tensor should occur at most once in the reduction
dimensions as the implementation does not remove duplicates.

The order of the values in the reduction dimensions does not affect the
results, but the code may execute faster if you list the dimensions in
increasing order.

Example: Reduction along one dimension.

    // Create a tensor of 2 dimensions
    Eigen::Tensor<int, 2> a(2, 3);
    a.setValues({{1, 2, 3}, {6, 5, 4}});
    // Reduce it along the second dimension (1)...
    Eigen::array<int, 1> dims({1 /* dimension to reduce */});
    // ...using the "maximum" operator.
    // The result is a tensor with one dimension.  The size of
    // that dimension is the same as the first (non-reduced) dimension of a.
    Eigen::Tensor<int, 1> b = a.maximum(dims);
    cout << "a" << endl << a << endl << endl;
    cout << "b" << endl << b << endl << endl;
    =>
    a
    1 2 3
    6 5 4

    b
    3
    6

Example: Reduction along two dimensions.

    Eigen::Tensor<float, 3, Eigen::ColMajor> a(2, 3, 4);
    a.setValues({{{0.0f, 1.0f, 2.0f, 3.0f},
                  {7.0f, 6.0f, 5.0f, 4.0f},
                  {8.0f, 9.0f, 10.0f, 11.0f}},
                 {{12.0f, 13.0f, 14.0f, 15.0f},
                  {19.0f, 18.0f, 17.0f, 16.0f},
                  {20.0f, 21.0f, 22.0f, 23.0f}}});
    // The tensor a has 3 dimensions.  We reduce along the
    // first 2, resulting in a tensor with a single dimension
    // of size 4 (the last dimension of a.)
    // Note that we pass the array of reduction dimensions
    // directly to the maximum() call.
    Eigen::Tensor<float, 1, Eigen::ColMajor> b =
        a.maximum(Eigen::array<int, 2>({0, 1}));
    cout << "b" << endl << b << endl << endl;
    =>
    b
    20
    21
    22
    23

#### Reduction along all dimensions

As a special case, if you pass no parameter to a reduction operation the
original tensor is reduced along *all* its dimensions.  The result is a
scalar, represented as a zero-dimension tensor.

    Eigen::Tensor<float, 3> a(2, 3, 4);
    a.setValues({{{0.0f, 1.0f, 2.0f, 3.0f},
                  {7.0f, 6.0f, 5.0f, 4.0f},
                  {8.0f, 9.0f, 10.0f, 11.0f}},
                 {{12.0f, 13.0f, 14.0f, 15.0f},
                  {19.0f, 18.0f, 17.0f, 16.0f},
                  {20.0f, 21.0f, 22.0f, 23.0f}}});
    // Reduce along all dimensions using the sum() operator.
    Eigen::Tensor<float, 0> b = a.sum();
    cout << "b" << endl << b << endl << endl;
    =>
    b
    276


### <Operation> sum(const Dimensions& new_dims)
### <Operation> sum()

Reduce a tensor using the sum() operator.  The resulting values
are the sum of the reduced values.

### <Operation> mean(const Dimensions& new_dims)
### <Operation> mean()

Reduce a tensor using the mean() operator.  The resulting values
are the mean of the reduced values.

### <Operation> maximum(const Dimensions& new_dims)
### <Operation> maximum()

Reduce a tensor using the maximum() operator.  The resulting values are the
largest of the reduced values.

### <Operation> minimum(const Dimensions& new_dims)
### <Operation> minimum()

Reduce a tensor using the minimum() operator.  The resulting values
are the smallest of the reduced values.

### <Operation> prod(const Dimensions& new_dims)
### <Operation> prod()

Reduce a tensor using the prod() operator.  The resulting values
are the product of the reduced values.

### <Operation> all(const Dimensions& new_dims)
### <Operation> all()
Reduce a tensor using the all() operator.  Casts tensor to bool and then checks
whether all elements are true.  Runs through all elements rather than
short-circuiting, so may be significantly inefficient.

### <Operation> any(const Dimensions& new_dims)
### <Operation> any()
Reduce a tensor using the any() operator.  Casts tensor to bool and then checks
whether any element is true.  Runs through all elements rather than
short-circuiting, so may be significantly inefficient.


### <Operation> reduce(const Dimensions& new_dims, const Reducer& reducer)

Reduce a tensor using a user-defined reduction operator.  See ```SumReducer```
in TensorFunctors.h for information on how to implement a reduction operator.


## Scan Operations

A *Scan* operation returns a tensor with the same dimensions as the original
tensor. The operation performs an inclusive scan along the specified
axis, which means it computes a running total along the axis for a given
reduction operation.
If the reduction operation corresponds to summation, then this computes the
prefix sum of the tensor along the given axis.

Example:
dd a comment to this line

    // Create a tensor of 2 dimensions
    Eigen::Tensor<int, 2> a(2, 3);
    a.setValues({{1, 2, 3}, {4, 5, 6}});
    // Scan it along the second dimension (1) using summation
    Eigen::Tensor<int, 2> b = a.cumsum(1);
    // The result is a tensor with the same size as the input
    cout << "a" << endl << a << endl << endl;
    cout << "b" << endl << b << endl << endl;
    =>
    a
    1 2 3
    6 5 4

    b
    1  3  6
    4  9 15

### <Operation> cumsum(const Index& axis)

Perform a scan by summing consecutive entries.

### <Operation> cumprod(const Index& axis)

Perform a scan by multiplying consecutive entries.


## Convolutions

### <Operation> convolve(const Kernel& kernel, const Dimensions& dims)

Returns a tensor that is the output of the convolution of the input tensor with the kernel,
along the specified dimensions of the input tensor. The dimension size for dimensions of the output tensor
which were part of the convolution will be reduced by the formula:
output_dim_size = input_dim_size - kernel_dim_size + 1 (requires: input_dim_size >= kernel_dim_size).
The dimension sizes for dimensions that were not part of the convolution will remain the same.
Performance of the convolution can depend on the length of the stride(s) of the input tensor dimension(s) along which the
convolution is computed (the first dimension has the shortest stride for ColMajor, whereas RowMajor's shortest stride is
for the last dimension).

    // Compute convolution along the second and third dimension.
    Tensor<float, 4, DataLayout> input(3, 3, 7, 11);
    Tensor<float, 2, DataLayout> kernel(2, 2);
    Tensor<float, 4, DataLayout> output(3, 2, 6, 11);
    input.setRandom();
    kernel.setRandom();

    Eigen::array<ptrdiff_t, 2> dims({1, 2});  // Specify second and third dimension for convolution.
    output = input.convolve(kernel, dims);

    for (int i = 0; i < 3; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 6; ++k) {
          for (int l = 0; l < 11; ++l) {
            const float result = output(i,j,k,l);
            const float expected = input(i,j+0,k+0,l) * kernel(0,0) +
                                   input(i,j+1,k+0,l) * kernel(1,0) +
                                   input(i,j+0,k+1,l) * kernel(0,1) +
                                   input(i,j+1,k+1,l) * kernel(1,1);
            VERIFY_IS_APPROX(result, expected);
          }
        }
      }
    }


## Geometrical Operations

These operations return a Tensor with different dimensions than the original
Tensor.  They can be used to access slices of tensors, see them with different
dimensions, or pad tensors with additional data.

### <Operation> reshape(const Dimensions& new_dims)

Returns a view of the input tensor that has been reshaped to the specified
new dimensions.  The argument new_dims is an array of Index values.  The
rank of the resulting tensor is equal to the number of elements in new_dims.

The product of all the sizes in the new dimension array must be equal to
the number of elements in the input tensor.

    // Increase the rank of the input tensor by introducing a new dimension
    // of size 1.
    Tensor<float, 2> input(7, 11);
    array<int, 3> three_dims{{7, 11, 1}};
    Tensor<float, 3> result = input.reshape(three_dims);

    // Decrease the rank of the input tensor by merging 2 dimensions;
    array<int, 1> one_dim{{7 * 11}};
    Tensor<float, 1> result = input.reshape(one_dim);

This operation does not move any data in the input tensor, so the resulting
contents of a reshaped Tensor depend on the data layout of the original Tensor.

For example this is what happens when you ```reshape()``` a 2D ColMajor tensor
to one dimension:

    Eigen::Tensor<float, 2, Eigen::ColMajor> a(2, 3);
    a.setValues({{0.0f, 100.0f, 200.0f}, {300.0f, 400.0f, 500.0f}});
    Eigen::array<Eigen::DenseIndex, 1> one_dim({3 * 2});
    Eigen::Tensor<float, 1, Eigen::ColMajor> b = a.reshape(one_dim);
    cout << "b" << endl << b << endl;
    =>
    b
      0
    300
    100
    400
    200
    500

This is what happens when the 2D Tensor is RowMajor:

    Eigen::Tensor<float, 2, Eigen::RowMajor> a(2, 3);
    a.setValues({{0.0f, 100.0f, 200.0f}, {300.0f, 400.0f, 500.0f}});
    Eigen::array<Eigen::DenseIndex, 1> one_dim({3 * 2});
    Eigen::Tensor<float, 1, Eigen::RowMajor> b = a.reshape(one_dim);
    cout << "b" << endl << b << endl;
    =>
    b
      0
    100
    200
    300
    400
    500

The reshape operation is a lvalue. In other words, it can be used on the left
side of the assignment operator.

The previous example can be rewritten as follow:

    Eigen::Tensor<float, 2, Eigen::ColMajor> a(2, 3);
    a.setValues({{0.0f, 100.0f, 200.0f}, {300.0f, 400.0f, 500.0f}});
    Eigen::array<Eigen::DenseIndex, 2> two_dim({2, 3});
    Eigen::Tensor<float, 1, Eigen::ColMajor> b;
    b.reshape(two_dim) = a;
    cout << "b" << endl << b << endl;
    =>
    b
      0
    300
    100
    400
    200
    500

Note that "b" itself was not reshaped but that instead the assignment is done to
the reshape view of b.


### <Operation> shuffle(const Shuffle& shuffle)

Returns a copy of the input tensor whose dimensions have been
reordered according to the specified permutation. The argument shuffle
is an array of Index values. Its size is the rank of the input
tensor. It must contain a permutation of 0, 1, ..., rank - 1. The i-th
dimension of the output tensor equals to the size of the shuffle[i]-th
dimension of the input tensor. For example:

    // Shuffle all dimensions to the left by 1.
    Tensor<float, 3> input(20, 30, 50);
    // ... set some values in input.
    Tensor<float, 3> output = input.shuffle({1, 2, 0})

    eigen_assert(output.dimension(0) == 30);
    eigen_assert(output.dimension(1) == 50);
    eigen_assert(output.dimension(2) == 20);

Indices into the output tensor are shuffled accordingly to formulate
indices into the input tensor. For example, one can assert in the above
code snippet that:

    eigen_assert(output(3, 7, 11) == input(11, 3, 7));

In general, one can assert that

    eigen_assert(output(..., indices[shuffle[i]], ...) ==
                 input(..., indices[i], ...))

The shuffle operation results in a lvalue, which means that it can be assigned
to. In other words, it can be used on the left side of the assignment operator.

Let's rewrite the previous example to take advantage of this feature:

    // Shuffle all dimensions to the left by 1.
    Tensor<float, 3> input(20, 30, 50);
    // ... set some values in input.
    Tensor<float, 3> output(30, 50, 20);
    output.shuffle({2, 0, 1}) = input;


### <Operation> stride(const Strides& strides)

Returns a view of the input tensor that strides (skips stride-1
elements) along each of the dimensions.  The argument strides is an
array of Index values.  The dimensions of the resulting tensor are
ceil(input_dimensions[i] / strides[i]).

For example this is what happens when you ```stride()``` a 2D tensor:

    Eigen::Tensor<int, 2> a(4, 3);
    a.setValues({{0, 100, 200}, {300, 400, 500}, {600, 700, 800}, {900, 1000, 1100}});
    Eigen::array<Eigen::DenseIndex, 2> strides({3, 2});
    Eigen::Tensor<int, 2> b = a.stride(strides);
    cout << "b" << endl << b << endl;
    =>
    b
       0   200
     900  1100

It is possible to assign a tensor to a stride:
    Tensor<float, 3> input(20, 30, 50);
    // ... set some values in input.
    Tensor<float, 3> output(40, 90, 200);
    output.stride({2, 3, 4}) = input;


### <Operation> slice(const StartIndices& offsets, const Sizes& extents)

Returns a sub-tensor of the given tensor. For each dimension i, the slice is
made of the coefficients stored between offset[i] and offset[i] + extents[i] in
the input tensor.

    Eigen::Tensor<int, 2> a(4, 3);
    a.setValues({{0, 100, 200}, {300, 400, 500},
                 {600, 700, 800}, {900, 1000, 1100}});
    Eigen::array<int, 2> offsets = {1, 0};
    Eigen::array<int, 2> extents = {2, 2};
    Eigen::Tensor<int, 1> slice = a.slice(offsets, extents);
    cout << "a" << endl << a << endl;
    =>
    a
       0   100   200
     300   400   500
     600   700   800
     900  1000  1100
    cout << "slice" << endl << slice << endl;
    =>
    slice
     300   400
     600   700


### <Operation> chip(const Index offset, const Index dim)

A chip is a special kind of slice. It is the subtensor at the given offset in
the dimension dim. The returned tensor has one fewer dimension than the input
tensor: the dimension dim is removed.

For example, a matrix chip would be either a row or a column of the input
matrix.

    Eigen::Tensor<int, 2> a(4, 3);
    a.setValues({{0, 100, 200}, {300, 400, 500},
                 {600, 700, 800}, {900, 1000, 1100}});
    Eigen::Tensor<int, 1> row_3 = a.chip(2, 0);
    Eigen::Tensor<int, 1> col_2 = a.chip(1, 1);
    cout << "a" << endl << a << endl;
    =>
    a
       0   100   200
     300   400   500
     600   700   800
     900  1000  1100
    cout << "row_3" << endl << row_3 << endl;
    =>
    row_3
       600   700   800
    cout << "col_2" << endl << col_2 << endl;
    =>
    col_2
       100   400   700    1000

It is possible to assign values to a tensor chip since the chip operation is a
lvalue. For example:

    Eigen::Tensor<int, 1> a(3);
    a.setValues({{100, 200, 300}});
    Eigen::Tensor<int, 2> b(2, 3);
    b.setZero();
    b.chip(0, 0) = a;
    cout << "a" << endl << a << endl;
    =>
    a
     100
     200
     300
    cout << "b" << endl << b << endl;
    =>
    b
       100   200   300
         0     0     0


### <Operation> reverse(const ReverseDimensions& reverse)

Returns a view of the input tensor that reverses the order of the coefficients
along a subset of the dimensions.  The argument reverse is an array of boolean
values that indicates whether or not the order of the coefficients should be
reversed along each of the dimensions.  This operation preserves the dimensions
of the input tensor.

For example this is what happens when you ```reverse()``` the first dimension
of a 2D tensor:

    Eigen::Tensor<int, 2> a(4, 3);
    a.setValues({{0, 100, 200}, {300, 400, 500},
                {600, 700, 800}, {900, 1000, 1100}});
    Eigen::array<bool, 2> reverse({true, false});
    Eigen::Tensor<int, 2> b = a.reverse(reverse);
    cout << "a" << endl << a << endl << "b" << endl << b << endl;
    =>
    a
       0   100   200
     300   400   500
     600   700   800
     900  1000  1100
    b
     900  1000  1100
     600   700   800
     300   400   500
       0   100   200


### <Operation> broadcast(const Broadcast& broadcast)

Returns a view of the input tensor in which the input is replicated one to many
times.
The broadcast argument specifies how many copies of the input tensor need to be
made in each of the dimensions.

    Eigen::Tensor<int, 2> a(2, 3);
    a.setValues({{0, 100, 200}, {300, 400, 500}});
    Eigen::array<int, 2> bcast({3, 2});
    Eigen::Tensor<int, 2> b = a.broadcast(bcast);
    cout << "a" << endl << a << endl << "b" << endl << b << endl;
    =>
    a
       0   100   200
     300   400   500
    b
       0   100   200    0   100   200
     300   400   500  300   400   500
       0   100   200    0   100   200
     300   400   500  300   400   500
       0   100   200    0   100   200
     300   400   500  300   400   500

### <Operation> concatenate(const OtherDerived& other, Axis axis)

TODO

### <Operation>  pad(const PaddingDimensions& padding)

Returns a view of the input tensor in which the input is padded with zeros.

    Eigen::Tensor<int, 2> a(2, 3);
    a.setValues({{0, 100, 200}, {300, 400, 500}});
    Eigen::array<pair<int, int>, 2> paddings;
    paddings[0] = make_pair(0, 1);
    paddings[1] = make_pair(2, 3);
    Eigen::Tensor<int, 2> b = a.pad(paddings);
    cout << "a" << endl << a << endl << "b" << endl << b << endl;
    =>
    a
       0   100   200
     300   400   500
    b
       0     0     0    0
       0     0     0    0
       0   100   200    0
     300   400   500    0
       0     0     0    0
       0     0     0    0
       0     0     0    0


### <Operation>  extract_patches(const PatchDims& patch_dims)

Returns a tensor of coefficient patches extracted from the input tensor, where
each patch is of dimension specified by 'patch_dims'. The returned tensor has
one greater dimension than the input tensor, which is used to index each patch.
The patch index in the output tensor depends on the data layout of the input
tensor: the patch index is the last dimension ColMajor layout, and the first
dimension in RowMajor layout.

For example, given the following input tensor:

  Eigen::Tensor<float, 2, DataLayout> tensor(3,4);
  tensor.setValues({{0.0f, 1.0f, 2.0f, 3.0f},
                    {4.0f, 5.0f, 6.0f, 7.0f},
                    {8.0f, 9.0f, 10.0f, 11.0f}});

  cout << "tensor: " << endl << tensor << endl;
=>
tensor:
 0   1   2   3
 4   5   6   7
 8   9  10  11

Six 2x2 patches can be extracted and indexed using the following code:

  Eigen::Tensor<float, 3, DataLayout> patch;
  Eigen::array<ptrdiff_t, 2> patch_dims;
  patch_dims[0] = 2;
  patch_dims[1] = 2;
  patch = tensor.extract_patches(patch_dims);
  for (int k = 0; k < 6; ++k) {
    cout << "patch index: " << k << endl;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        if (DataLayout == ColMajor) {
          cout << patch(i, j, k) << " ";
        } else {
          cout << patch(k, i, j) << " ";
        }
      }
      cout << endl;
    }
  }

This code results in the following output when the data layout is ColMajor:

patch index: 0
0 1
4 5
patch index: 1
4 5
8 9
patch index: 2
1 2
5 6
patch index: 3
5 6
9 10
patch index: 4
2 3
6 7
patch index: 5
6 7
10 11

This code results in the following output when the data layout is RowMajor:
(NOTE: the set of patches is the same as in ColMajor, but are indexed differently).

patch index: 0
0 1
4 5
patch index: 1
1 2
5 6
patch index: 2
2 3
6 7
patch index: 3
4 5
8 9
patch index: 4
5 6
9 10
patch index: 5
6 7
10 11

### <Operation>  extract_image_patches(const Index patch_rows, const Index patch_cols,
                          const Index row_stride, const Index col_stride,
                          const PaddingType padding_type)

Returns a tensor of coefficient image patches extracted from the input tensor,
which is expected to have dimensions ordered as follows (depending on the data
layout of the input tensor, and the number of additional dimensions 'N'):

*) ColMajor
1st dimension: channels (of size d)
2nd dimension: rows (of size r)
3rd dimension: columns (of size c)
4th-Nth dimension: time (for video) or batch (for bulk processing).

*) RowMajor (reverse order of ColMajor)
1st-Nth dimension: time (for video) or batch (for bulk processing).
N+1'th dimension: columns (of size c)
N+2'th dimension: rows (of size r)
N+3'th dimension: channels (of size d)

The returned tensor has one greater dimension than the input tensor, which is
used to index each patch. The patch index in the output tensor depends on the
data layout of the input tensor: the patch index is the 4'th dimension in
ColMajor layout, and the 4'th from the last dimension in RowMajor layout.

For example, given the following input tensor with the following dimension
sizes:
 *) depth:   2
 *) rows:    3
 *) columns: 5
 *) batch:   7

  Tensor<float, 4> tensor(2,3,5,7);
  Tensor<float, 4, RowMajor> tensor_row_major = tensor.swap_layout();

2x2 image patches can be extracted and indexed using the following code:

*) 2D patch: ColMajor (patch indexed by second-to-last dimension)
  Tensor<float, 5> twod_patch;
  twod_patch = tensor.extract_image_patches<2, 2>();
  // twod_patch.dimension(0) == 2
  // twod_patch.dimension(1) == 2
  // twod_patch.dimension(2) == 2
  // twod_patch.dimension(3) == 3*5
  // twod_patch.dimension(4) == 7

*) 2D patch: RowMajor (patch indexed by the second dimension)
  Tensor<float, 5, RowMajor> twod_patch_row_major;
  twod_patch_row_major = tensor_row_major.extract_image_patches<2, 2>();
  // twod_patch_row_major.dimension(0) == 7
  // twod_patch_row_major.dimension(1) == 3*5
  // twod_patch_row_major.dimension(2) == 2
  // twod_patch_row_major.dimension(3) == 2
  // twod_patch_row_major.dimension(4) == 2

## Special Operations

### <Operation> cast<T>()

Returns a tensor of type T with the same dimensions as the original tensor.
The returned tensor contains the values of the original tensor converted to
type T.

    Eigen::Tensor<float, 2> a(2, 3);
    Eigen::Tensor<int, 2> b = a.cast<int>();

This can be useful for example if you need to do element-wise division of
Tensors of integers.  This is not currently supported by the Tensor library
but you can easily cast the tensors to floats to do the division:

    Eigen::Tensor<int, 2> a(2, 3);
    a.setValues({{0, 1, 2}, {3, 4, 5}});
    Eigen::Tensor<int, 2> b =
        (a.cast<float>() / a.constant(2).cast<float>()).cast<int>();
    cout << "a" << endl << a << endl << endl;
    cout << "b" << endl << b << endl << endl;
    =>
    a
    0 1 2
    3 4 5

    b
    0 0 1
    1 2 2


### <Operation>     eval()

TODO


## Representation of scalar values

Scalar values are often represented by tensors of size 1 and rank 1. It would be
more logical and user friendly to use tensors of rank 0 instead. For example
Tensor<T, N>::maximum() currently returns a Tensor<T, 1>. Similarly, the inner
product of 2 1d tensors (through contractions) returns a 1d tensor. In the
future these operations might be updated to return 0d tensors instead.

## Limitations

*   The number of tensor dimensions is currently limited to 250 when using a
    compiler that supports cxx11. It is limited to only 5 for older compilers.
*   The IndexList class requires a cxx11 compliant compiler. You can use an
    array of indices instead if you don't have access to a modern compiler.
*   On GPUs only floating point values are properly tested and optimized for.
*   Complex and integer values are known to be broken on GPUs. If you try to use
    them you'll most likely end up triggering a static assertion failure such as
    EIGEN_STATIC_ASSERT(packetSize > 1, YOU_MADE_A_PROGRAMMING_MISTAKE)