One standard method for viscoelastic (VE) finite element method is the discrete elastic-viscous stress-split method of Guenette and Fortin (J Non-Newt. Fluid Mech., 1995) which uses the Galerkin method for the momentum equation requiring a velocity-pressure space that satisfies the LBB condition, a Petrov-Galerkin method for the differential stress constitutive equation, and the addition of a least-square interpolation of the velocity gradient tensor to a continuous space. Adding the smoothed velocity gradient to the momentum equation acts to stabilize the stress space, allowing for low order stress interpolation. These simulations are computationally demanding since in addition to the standard Navier-Stokes unknowns, we must also solve a symmetric stress tensor and a non-symmetric velocity gradient tensor. Multimode calculations require even more unknowns, adding a stress tensor for each mode. The matrices arising from the discretized problem are poorly conditioned, requiring direct Gaussian elimination methods. Here we investigate methods to speed up multimode viscoelastic flow calculation in both 2D and 3D, including free and moving boundaries represented by an arbitrary-Lagrangian-Eulerian method. Applications of interest include extrusion, coextrusion of layered polymers, and blade coating of polymeric solution on flexible substrates. Results have shown that decoupling the momentum equations from the stress and velocity gradient equations can greatly reduce solution times. Stabilization of the momentum equation using Dorhmann-Bochev stabilization allows for use of Krylov-based iterative solvers in place of direct solvers, even for mixed velocity-pressure formulations. The method is fully parallelized using MPI and a GMRES solver preconditioned with ILUT. Comparison of solution times and scalability for the monolithic and decoupled solve will be discussed.