Models and algorithms for elastoplasticity at finite deformations are well established. They typically entail a local integration of the elastoplastic model within a strain-driven format. The local plastic state is then characterised by a set of internal variables associated with each quadrature point. The classical approach to computational plasticity was devised to be efficient on central processing unit based computing architectures where floating-point operations were the primary bottleneck. By contrast, modern computing architectures, based on, e.g., graphical processing units, can only be effectively exploited by algorithms that minimise memory access. This requires a fundamental change in algorithm design philosophy and provides the primary motivation for matrix-free solvers. However, matrix-free methods are not widely adopted in solid mechanics due to the complex nonlinearities and the widespread use of low-order finite element approximations. Furthermore, the evolution of quantities characterising inelastic processes - including classical plasticity - is nearly always approximated at the local level of the quadrature point, inhibiting the use of matrix-free approaches. A key contribution of this work is the novel multifield formulation for finite strain plasticity well suited for matrix-free approaches. To this end, the balance of linear momentum, the flow relation, and the Karush--Kuhn--Tucker constraints are collectively cast in a variational format. Thus, in addition to the deformation, the plastic strain tensor and the consistency parameter are global degrees of freedom in the resulting spatially discrete problem that follows from a finite element discretisation. We adopt a finite strain plasticity formulation set in logarithmic strain space. The KKT inequality constraints on the evolution of plastic flow are recast as a variant of the nonlinear Fischer--Burmeister complementarity function. This circumvents the need for an active set search at the global level. The multifield approach results in a proliferation of global degrees of freedom. This is addressed by exploiting the block sparse structure of the algebraic system together with a tailored block matrix solver. Details of the extension of the coupled thermoplastic problem with contact for extreme applications such as friction welding are also provided.