Frequentist and likelihood methods of inference based on the multivariate skew-normal model encounter several technical difficulties with this model. In spite of the popularity of this class of densities, there are no broadly satisfactory solutions for estimation and testing problems. A general population Monte Carlo algorithm is proposed which: 1) exploits the latent structure stochastic representation of skew-normal random variables to provide a full Bayesian analysis of the model and 2) accounts for the presence of constraints in the parameter space. The proposed approach can be defined as weakly informative, since the prior distribution approximates the actual reference prior for the shape parameter vector. Results are compared with the existing classical solutions and the practical implementation of the algorithm is illustrated via a simulation study and a real data example. A generalization to the matrix variate regression model with skew-normal error is also presented.