A computational bottleneck during the solution to multiphase formulations of the incompressible Navier-Stokes equations is often during the implicit solution of the pressure-correction equation that results from operator-splitting methods. Since density is a coefficient in the pressure-correction equation, large variations or discontinuities among the phase densities greatly increase the condition number of the pressure-correction matrix and impede the convergence of iterative methods employed in its solution. To alleviate this shortcoming, the open-source multiphase code MFiX is interfaced with the linear solver library PETSc. Through an appropriate mapping of matrix and vector data structures between the two software, the access to a suite of robust, scalable, solver options in PETSc is obtained. Verification of the implementation of MFiX-PETSc is demonstrated through predictions that are identical to those obtained from MFiX's native solvers for a simple heat conduction case with a well-known solution. After verifying the framework, several cases were tested with MFiX-PETSc to analyze the performance of various solver and preconditioner combinations. For a low Reynolds number, flow over a cylinder case, applying right-side Block Jacobi preconditioning to the BiCGSTAB iterative solver in MFiX-PETSc was 28-40% faster than MFiX's native solver at the finest mesh resolution. Similarly, the left-side Block Jacobi preconditioner in MFiX-PETSc was 27--46% faster for the same fine meshing. Further assessments of these preconditioning options were then made for a fluidized bed problem involving different bed geometries, convergence tolerances, material densities, and inlet velocities. For a three-dimensional geometry with uniform meshing, native MFiX was faster than MFiX-PETSc for each simulation. The difference in speed was minimized when a low density fluidization material (polypropylene) was used along with a higher order discretization scheme. With these settings, MFiX-PETSc was only 2-6% slower than native MFiX when right-side Block Jacobi preconditioning was employed. The fluidized bed was then represented by a two-dimensional geometry with fine meshing towards the center. When this bed was filled with glass beads, right-side Block Jacobi was 28% faster than MFiX's native solver, which was the largest speedup encountered throughout this 2D case.