// Declare some random vector and matrices.
satyr::vector<double> v(5);
satyr::matrix<double> A(5, 5), B(5, 5);
satyr::symmetric_matrix<double> S(5);
satyr::lower_triangular_matrix<double> L(5);

// The standard arithmetic operators and mathematical functions can be used to
// execute expression templates.
A = B + 2.0 * cos(S);

// Expressions involving only structural matrices are computed in an efficient
// manner that avoids unnecessary work.
S = sqrt(abs(S));  // computes only over a triangular portion of the matrix.

// Additionally execution policies can be applied to parallelize or vectorize
// the computation.
A = sqrt(abs(L)) << satyr::parallel_v << satyr::simd_v;

// For parallelization, you can also specify a grainsize if the cost of
// managing tasks could potentially be more expensive than the computation
// itself.
A += cos(B) - as_diagonal_matrix(v) << satyr::grainsize{
         10};  // Don't create tasks with fewer than 10 iterations.

// Additionally, you can declare numerical arrays of arbitrary dimension.
satyr::n_array<float, 3> H(5, 2, 6);

// And all array-like objects support indexing and slicing.
A(0, 0) -= 5;
A += H(satyr::all_v, 1, satyr::range{1, 6});

// There is a wrapper for many BLAS-LAPACK functions.
auto C = product(A, B);     // calls gemm.
auto w = left_solve(L, v);  // calls trsv.