  Let $p(\cdot):\ \mathbb{R}^n\to(0,\infty]$ be a variable exponent function
satisfying the globally log-H\"{o}lder continuous condition and $A$ a general
expansive matrix on $\mathbb{R}^n$. Let $H_A^{p(\cdot)}(\mathbb{R}^n)$ be the
variable anisotropic Hardy space associated with $A$ defined via the
non-tangential grand maximal function. In this article, via the known atomic
characterization of $H_A^{p(\cdot)}(\mathbb{R}^n)$, the author establishes its
molecular characterization with the known best possible decay of molecules. As
an application, the author obtains a criterion on the boundedness of linear
operators on $H_A^{p(\cdot)}(\mathbb{R}^n)$, which is used to prove the
boundedness of anisotropic Calder\'on-Zygmund operators on
$H_A^{p(\cdot)}(\mathbb{R}^n)$. In addition, the boundedness of anisotropic
Calder\'on-Zygmund operators from $H_A^{p(\cdot)}(\mathbb{R}^n)$ to the
variable Lebesgue space $L^{p(\cdot)}(\mathbb{R}^n)$ is also presented. All
these results are new even in the classical isotropic setting.
