In order to accelerate numerical integration, improve long time simulation accuracy, and sample statistical distributions by dynamics, we develop geometric multiscale integrators. These integrators employ coarse steps that do not resolve the fast scale in the system; nevertheless, they capture the effective contribution of the fast dynamics. They works for a broad class of systems, including stiff ODEs, SDEs and PDEs, by requiring no identification of underlying slow variables or processes. They also numerically preserve intrinsic geometric structures (e.g., symplecticity, conservation laws, and invariant distribution); such preservations lead to not only improved long time accuracy, but also an option of sampling statistical distributions by dynamics.