I see in the documentation that Intels' Poisson solver routines can be used to solve equations of the form
∇.∇u = f
Can they also (perhaps via some additional steps) be applied to an equation of the form
∇.ε∇u = f
The possible algorithm choices depend on whether ε is (i) a constant, (ii) position dependent, or (iii) a known function of u. In the first two cases, the finite-difference equations are still linear, and can be solved using Pardiso as in the case when ε = 1. In the third case, you may use Kirchoff's transformation, φ = \int ε du, which leads to ∇.∇φ = f; you solve for φ, and recover u from φ at the end.