In order to smooth point-based surface with noise, Bezier surface fitting method was applied to compute principal curvature and normal vector at each point. We construct an undirected weighted graph from scattered points, and a heat diffusion partial differential equation was defined on this graph. Considering principal curvature as heat power, the smoothing process is realized by diffusing principal curvature along graph structure. We adapt the position of each point along normal direction according to the difference of curvature. Computing heat diffusion equation can be boiled down to the spectral decomposition of Laplacian matrix. We use parallel computing method to solve the spectral decomposition of Laplacian Matrix on GPU to improve the efficiency. Some examples show that our method works well on complex models with large scale point sets.