Abstract：Multiple solutions may be resulted from the inversion of 3D electrical resistivity and the calculation is time-consuming，which restricts the application of resistivity detection. An idea is presented to solve this problem，with the inequality constraints to the inversion equations and the parallel algorithms being applied. Based on the traditional smooth constraints，the relaxed variables were introduced to apply the inequality constraints to the inverse equations，which carry the upper and lower limits information of subsurface media. A new objective function was obtained using the primal dual interior point method. The resistivity values in the function were defined within the inequality constrains range，and the parameters were optimized within the feasible region defined by the constraints. The method theoretically suppressed the multiplicity of inverse solution. The parallel calculation algorithm for partial derivative matrix and the parallel algorithm for Cholesky decomposition to the overall coefficient matrix were designed，which speeds the inversion more than fifty percent. Based on the above research，the parallel 3D electrical resistivity inversion method with inequality constraint based on relax variables was developed，and numerical tests and engineering application were carried out. Results showed that the above method made full use of the inequality constrains and removed the false anomalies，suppressed the multiplicity and improved the accuracy and calculation efficiency.