  Models of fermions interacting with classical degrees of freedom are applied
to a large variety of systems in condensed matter physics. For this class of
models, Wei{\ss}e [Phys. Rev. Lett. {\bf 102}, 150604 (2009)] has recently
proposed a very efficient numerical method, called O($N$) Green-Function-Based
Monte Carlo (GFMC) method, where a kernel polynomial expansion technique is
used to avoid the full numerical diagonalization of the fermion Hamiltonian
matrix of size $N$, which usually costs O($N^3$) computational complexity.
Motivated by this background, in this paper we apply the GFMC method to the
double exchange model in three spatial dimensions. We mainly focus on the
implementation of GFMC method using both MPI on a CPU-based cluster and
Nvidia's Compute Unified Device Architecture (CUDA) programming techniques on a
GPU-based (Graphics Processing Unit based) cluster. The time complexity of the
algorithm and the parallel implementation details on the clusters are
discussed. We also show the performance scaling for increasing Hamiltonian
matrix size and increasing number of nodes, respectively. The performance
evaluation indicates that for a $32^3$ Hamiltonian a single GPU shows higher
performance equivalent to more than 30 CPU cores parallelized using MPI.
