The new model of the last glacial isostatic adjustment (GIA) is constrained simutaneously by the data sets from geodetic measurements, relative sea levels and seismic shear tomographic data. The viscosities in mantle for this model are allowed to vary both radially and laterally. Firstly, the scalling factor β is tentatively chosen for which the lateral perturbations in viscosity are computed based on the linear relation with seismic shear wave anomalies. These lateral perturbations in viscosity are then supperposed logarithmically on a laterally homogeneous reference viscosity model to give the total viscosities in the mantle layers. Based on this 3-D viscosity model, the GIA predictions are computed using the coupled-Laplace Finite-Element Method. Such processes will not end untill the predictions agree well with the observed data sets. The main results are as follows. Firstly, the laterally heterogeneous viscosity model (RF3L20(β =0. 4)) is determined from which the obvious lateral heterogeneities and the prominouced effects on GIA predictions are found. It is indicated the lateral heterogeneities are attributable not just to themal effects but also to lateral variations in chemistry and other factors. The RF3L20(β = 0. 4) model can be used to the modelling in mantle geodynamics. Secondly, the multiple present GIA rates are predicted that can be served as the important corrections for the monitoring of the tectonic motion, land water storage changes, ocean water mass variations and the non-equilibrium of the ice and snow.