Linear inverse Gaussian problems are traditionally solved using least squares-based inversion. The center of the posterior Gaussian probability distribution is often chosen as the solution to such problems, while the solution is in fact the posterior Gaussian probability distribution itself. We present an algorithm, based on direct sequential simulation, which can be used to efficiently draw samples of the posterior probability distribution for linear inverse problems. There is no Gaussian restriction on the distribution in the model parameter space, as inherent in traditional least squares-based algorithms. As data for linear inverse problems can be seen as weighed linear averages over some volume, block kriging can be used to perform both estimation (i.e. finding the center of the posterior Gaussian pdf) and simulation (drawing samples of the posterior Gaussian pdf). We present the kriging system which we use to implement a flexible GSLIB-based algorithm for solving linear inverse problems. We show how we implement such a simulation program conditioned to linear average data. The program is called VISIM as an acronym for Volume average Integration SIMulation. An effort has been made to make the program efficient, even for larger scale problems, and the computational efficiency and accuracy of the code is investigated. Using a synthetic cross-borehole tomography case study, we show how the program can be used to generate realizations of the a posteriori distributions (i.e. solutions) from a linear tomography problem. Both Gaussian and non-Gaussian a priori model parameter distributions are considered.