  This paper provides a general framework for Stein's density method for
multivariate continuous distributions. The approach associates to any
probability density function a canonical operator and Stein class, as well as
an infinite collection of operators and classes which we call standardizations.
These in turn spawn an entire family of Stein identities and characterizations
for any continuous distribution on $\mathbb{R}^d$, among which we highlight
those based on the score function and the Stein kernel. A feature of these
operators is that they do not depend on normalizing constants. A new definition
of Stein kernel is introduced and examined; integral formulas are obtained
through a connection with mass transport, as well as ready-to-use explicit
formulas for elliptical distributions. The flexibility of the kernels is used
to compare in Stein discrepancy (and therefore 2-Wasserstein distance) between
two normal distributions, Student and normal distributions, as well as two
normal-gamma distributions. Upper and lower bounds on the 1-Wasserstein
distance between continuous distributions are provided, and computed for a
variety of examples: comparison between different normal distributions
(improving on existing bounds in some regimes), posterior distributions with
different priors in a Bayesian setting (including logistic regression), centred
Azzalini--Dalla Valle distributions. Finally the notion of weak Stein equation
and weak Stein factors is introduced. Bounds for solutions of the weak Stein
equation are obtained for Lipschitz test functions if the distribution admits a
Poincar\'e constant. We use these bounds to compare different copulas on the
unit square in 1-Wasserstein distance.
