MR susceptometry provides a new approach to enhance contrast in MR imaging and quantify substances such as iron, calcium, blood oxygenation in various organs for the clinical diagnosis of many diseases. Susceptibility is closely related to the magnetic field inhomogeneity in MRI, and calculation of susceptibility from the measured magnetic field is an ill-posed inverse problem. Conjugate gradient least square (CGLS) method with Tikhonov regularization is a powerful tool to solve the inverse problems. However, the estimated quantity in CGLS method is usually in one-dimensional (e.g., a column vector) while the most of MR data are in the form of three dimensions. In this work, the least square problem is modified to enable the calculation of susceptibility directly from three-dimensional MR phase maps, which has the benefit of reducing the size of associated matrices and usage of computer memories. Numerical simulations were used to find the proper regularization parameters and study the influence of different noise levels of the magnetic field on the regularization parameters. Experiments of superparamagnetic iron oxide phantom were also conducted. Both of the results demonstrate the validation and accuracy of this method.