Implicit integration of the viscous term can significantly improve performance in computational fluid dynamics for highly viscous fluids such as lava. We show improvements over our previous proposal for semi-implicit viscous integration in Smoothed Particle Hydrodynamics, extending it to support a wider range of boundary models. Due to the resulting loss of matrix symmetry, a key advancement is a more robust version of the biconjugate gradient stabilized method to solve the linear systems, that is also better suited for parallelization in both shared-memory and distributed-memory systems. The advantages of the new solver are demostrated in applications with both Newtonian and non-Newtonian fluids, covering both the numerical aspect (improved convergence thanks to the possibility to use more accurate boundary model) and the computing aspect (with excellent strong scaling and satisfactory weak scaling).