  We propose an efficient Monte Carlo algorithm for simulating a
``hardly-relaxing" system, in which many replicas with different temperatures
are simultaneously simulated and a virtual process exchanging configurations of
these replica is introduced. This exchange process is expected to let the
system at low temperatures escape from a local minimum. By using this algorithm
the three-dimensional $\pm J$ Ising spin glass model is studied. The ergodicity
time in this method is found much smaller than that of the multi-canonical
method. In particular the time correlation function almost follows an
exponential decay whose relaxation time is comparable to the ergodicity time at
low temperatures. It suggests that the system relaxes very rapidly through the
exchange process even in the low temperature phase.
