A three-dimensional finite element computer program for thermal stress analyses was developed to deal with the cracking of trigonal class 3m oxide bulk single crystals, such as lithium niobate (LN) and lithium tantalate (LT) during Czochralski growth. A tensor transformation technique was used to obtain the elastic constant matrix and thermal expansion coefficient vector corresponding to an arbitrary pulling direction. The anisotropy of the elastic constants and thermal expansion coefficients as well as their temperature dependence were considered in the program. Using this program, we performed thermal stress analyses of LN bulk single crystals with three types of temperature distribution or the shape of crystal-melt interface. The analyses were performed both for the pulling directions of the a-axis  and c-axis . The results show that the thermal stress is very sensitive to the shape of the crystal-melt interface, and that the thermal stress is lower in the c-axis pulling than in the a-axis pulling. The relation between the thermal stress and the cracking and the quality of a bulk single crystal is also discussed.