Hi Yantao,
This is certainly possible to do in a way that scales exponentially with the number of sites in the region "A" whose entanglement you want to compute.
You can see this previous question as one discussion with some example code for the case of a two site region "A":
http://itensor.org/support/229/evaluate-block-entanglement-blocks-that-extend-the-lattice?show=229#q229
The basic idea is just to explicitly trace rho = psi^2 on all of the sites of region B, which is a trivial operation if the gauge center site is located inside of region A. Then if region A has Na sites, you will get its reduced density matrix as a tensor with 2*Na indices. You can diagonalize this tensor using the diagHermitian function.
There may be a better scaling algorithm in principle, but there is not one which has been published as far as I know. (E.g. one ambitious algorithm would be to try to find the eigenvalues of rho one-by-one using an iterative eigensolver such as Lanczos, where one would multiply an MPS onto rho and construct rho tensor-by-tensor each time to keep the scaling polynomial. But one might need to compute a lot of eigenvalues to get a good estimate of the entropy.)
Best regards,
Miles