Nonnegative matrix factorization (NMF) is a powerful technique for dimension reduction, extracting latent factors and learning part-based representation. For large datasets, NMF performance depends on some major issues: fast algorithms, fully parallel distributed feasibility and limited internal memory. This research aims to design a fast fully parallel and distributed algorithm using limited internal memory to reach high NMF performance for large datasets. In particular, we propose a flexible accelerated algorithm for NMF with all its $L_1$ $L_2$ regularized variants based on full decomposition, which is a combination of an anti-lopsided algorithm and a fast block coordinate descent algorithm. The proposed algorithm takes advantages of both these algorithms to achieve a linear convergence rate of $\mathcal{O}(1-\frac{1}{||Q||_2})^k$ in optimizing each factor matrix when fixing the other factor one in the sub-space of passive variables, where $r$ is the number of latent components; where $\sqrt{r} \leq ||Q||_2 \leq r$. In addition, the algorithm can exploit the data sparseness to run on large datasets with limited internal memory of machines. Furthermore, our experimental results are highly competitive with 7 state-of-the-art methods about three significant aspects of convergence, optimality and average of the iteration number. Therefore, the proposed algorithm is superior to fast block coordinate descent methods and accelerated methods.