  Linear algebraic expressions are the essence of many computationally
intensive problems, including scientific simulations and machine learning
applications. However, translating high-level formulations of these expressions
to efficient machine-level representations is far from trivial: developers
should be assisted by automatic optimization tools so that they can focus their
attention on high-level problems, rather than low-level details. The
tractability of these optimizations is highly dependent on the choice of the
primitive constructs in terms of which the computations are to be expressed. In
this work we propose to describe operations on multi-dimensional arrays using a
selection of higher-order functions, inspired by functional programming, and we
present rewrite rules for these such that they can be automatically optimized
for modern hierarchical and heterogeneous architectures. Using this formalism
we systematically construct and analyse different subdivisions and permutations
of the dense matrix multiplication problem.
