  We describe an algorithm to decompose deep images of Active Galactic Nuclei
into host galaxy and nuclear components. Currently supported are three galaxy
models: A de-Vaucouleurs spheroidal, an exponential disc, and a two-component
disc+bulge model. Key features of the method are: (semi-)analytic
representation of a possibly spatially variable point-spread function; full
two-dimensional convolution of the model galaxy using gradient-controlled
adaptive subpixelling; multiple iteration scheme. The code is computationally
efficient and versatile for a wide range of applications. The quantitative
performance is measured by analysing simulated imaging data. We also present
examples of the application of the method to small test samples of nearby
Seyfert 1 galaxies and quasars at redshifts z < 0.35.
