  Variational formulations of time-dependent PDEs in space and time yield
$(d+1)$-dimensional problems to be solved numerically. This increases the
number of unknowns as well as the storage amount. On the other hand, this
approach enables adaptivity in space and time as well as model reduction w.r.t.
both type of variables. In this paper, we show that matrix oriented techniques
can significantly reduce the computational timings for solving the arising
linear systems outperforming both time-stepping schemes and other solvers.
