Computational multiscale methods such as the FE2 technique (Feyel, 1999) come along with large demands in both CPU time and memory. In order to significantly reduce the computational cost of multiscale methods the authors recently proposed a hybrid computational homogenization method for visco-plastic materials using a reduced basis approach in a mixed variational formulation (Fritzen and Leuschner, 2013). In the present contribution two extensions of the method are introduced: First, the previous proposal is extended by allowing for heterogeneous hardening variables instead of piecewise constant fields. This leads to an improved accuracy of the method. Second, a massively parallel GPU implementation of the algorithm using Nvidia's CUDA framework is presented. The GPU subroutines for the batched linear algebraic operations are integrated into a specialized library in order to facilitate its use. The impact of the heterogeneous hardening states on the accuracy and the performance gains obtained from the dedicated GPU implementation are illustrated by means of numerical examples. An overall speedup in the order of 10(4) with respect to a high performance finite element implementation is achieved while preserving good accuracy of the predicted nonlinear material response. (C) 2014 Elsevier B.V. All rights reserved.