The Poisson-Boltzmann (PB) equation is an effective model for the electrostatics analysis of solvated biomolecules. The nonlinearity associated with the PB equation is critical when the underlying electrostatic potential is strong, but is extremely difficult to solve numerically. Recently, we have developed several operator splitting methods to efficiently and stably solve the nonlinear PB equation in a pseudo-transient continuation approach. The operator splitting framework enables an analytical integration of the nonlinear term that suppresses the nonlinear instability. Fully implicit alternating direction implicit (ADI) and locally one-dimensional (LOD) schemes are formulated to reduce 3D linear systems into 1D ones. Central finite difference discretization is conducted in space, which yields a tridiagonal structure for 1D systems. The LOD methods are found to be unconditionally stable in dealing with real proteins with source singularities and nonsmooth solutions. By using a large time increment, the optimized LOD method becomes over 20 times faster in electrostatic free energy analysis, without sacrificing the precision. The resulting PB solver scales linearly with respect to the number of atoms, and is promising for studying large macromolecules.