The neutrinoless double beta decay is of special importance in determining the fundamental properties of neutrinos. The nuclear matrix element of this decay must be evaluated in a sufficient accuracy, and the shell-model calculation can make contributions to this end. This is because the shell-model calculations incorporate basically all correlations into the wave functions of the initial and final states of the decay, and the accuracy of the calculation can be investigated by referring to other observables. I will report results obtained by recent large-scale shell-model calculations on 76Ge and 136Xe as well as their daughter nuclei 76Se and 136Ba. Here the large-scale shell-model calculations mean those by Monte Carlo Shell Model at its most advanced edition. The results are not away from the ranges of earlier studies, but are rather on the edges of smaller values. I will also discuss why such smaller values arise as natural consequences of basic features of the wave functions.