A Hamiltonian decomposition of a regular graph is a partition of its edge set into Hamiltonian cycles. The problem of finding edge-disjoint Hamiltonian cycles in a given regular graph has many applications in combinatorial optimization and operations research. Our motivation for this problem comes from the field of polyhedral combinatorics, as a sufficient condition for vertex nonadjacency in the 1-skeleton of the traveling salesperson polytope can be formulated as the Hamiltonian decomposition problem in a 4-regular multigraph with one forbidden decomposition. In our approach, the algorithm starts by solving the relaxed 2-matching problem, then iteratively generates subtour elimination constraints for all subtours in the solution and solves the corresponding ILP-model to optimality. The procedure is enhanced by the local search heuristic based on chain edge fixing and cycle merging operations. In the computational experiments, the iterative ILP algorithm showed comparable results with the previously known heuristics on undirected multigraphs and significantly better performance on directed multigraphs.