We report on analytical and numerical studies of the magnetic quantum oscillations---the de Haas--van Alphen (dHvA) effect---in a model of spinless electrons on a two-dimensional (2D) square lattice in a magnetic field $B$ known as the Azbel-Hofstadter (AH) problem. The numerical simulations of the dHvA oscillations are based on the magnetic band structure, which is obtained by numerical diagonalization of the one-electron model Hamiltonian on a square lattice in the magnetic field. At zero magnetic field and for a filling of 0.55 electrons per unit cell, the shape of the Fermi surface (FS) is close to a chessboard in which white and black cells are related to holes and electrons, respectively. The numerical simulations of the dHvA oscillations for this FS clearly display the magnetic breakdown (MB) features in the Fourier spectrum: the combinational and forbidden frequencies, which cannot be explained within the standard Lifshitz-Kosevich theory. Because of the idealized conditions of our numerical simulation (zero temperature, ideal crystal lattice) the only reason for deviation of the dHvA oscillations from the sawtooth magnetization of a 2D electron gas is the fractal magnetic band structure of the AH problem. We have also developed a semiclassical analytical approach to the dHvA oscillations taking account of the MB at the corners of the hole and electron orbits of the FS, the probability of which $0\ensuremath{\le}W(B)=\mathrm{exp}(\ensuremath{-}{B}_{\mathrm{MB}}∕B)\ensuremath{\le}1$ depends on $B$ and the MB field ${B}_{\mathrm{MB}}$. We have shown that the dHvA oscillations in this case are completely determined by the quantization of the closed hole and electron orbits and by special MB factors for these orbits ${I}^{e(h)}(p)$ ($p=1,2,3,\dots{}$, is the harmonic index). The factors ${I}^{e(h)}(p)$ not only strongly modulate the amplitudes of harmonics in the magnetization but act also through the changes in the chemical potential (CP) oscillations. The latter are very important: they control the shape of the magnetization (which can be close either to a direct or an inverted sawtoothlike profile) and add forbidden peaks in the Fourier spectrum. In the case of a fixed CP the forbidden peaks vanish and the magnetization profile is close to the direct sawtooth. As a function of the MB probability $W(B)$, the MB factors ${I}^{e(h)}(p)$ are polynomials of degree $p$ with the following fixed values at the ends: 1 at $W(B)=0$ and 0 at $W(B)=1$. Our analytic expression for the magnetization as a function of magnetic field $B$ agrees pretty well with our numerical simulations using only one fitting parameter---the MB field ${B}_{\mathrm{MB}}$. Small deviations are due to the fractal fine structure of the magnetic bands.