Accurate prediction of instabilities in viscoelastic flows presents one of the most challenging problems in non-Newtonian fluid mechanics and has a crucial importance for polymer processing where output quality constraints require that operating conditions should be maintained in the stable flow regime. Therefore, the flow instabilities are the primary constraint for the processing speed in many industrial polymer forming processes such as extrusion, wire coating, blow molding, sheet-formation, etc. These viscoelastic instabilities occur in the creeping motion of non-Newtonian polymeric liquids and are entirely absent in the corresponding motion of Newtonian liquids. It has been suggested that for shear-dominated flows one destabilizing mechanism is a combination of streamline curvature and large elastic stresses along the streamlines, giving rise to an extra hoop stress in a direction normal to the streamlines. During the last decade, significant progress has been made in the development of the numerical schemes for the simulation of viscoelastic fluid flows. However, realistic three-dimensional simulations around complex configurations are still challenging in terms of accuracy, required computer power, stability and convergence. The stability analysis of these fluids is even more challenging due to the often prohibitively large generalized eigenvalue problem (GEVP) resulting from the linear stability analysis. In this work, a parallel adaptive unstructured finite volume method is presented for analysis of the stability of two-dimensional steady Oldroyd-B fluid flows to small amplitude three-dimensional perturbations. A semi-staggered dilation-free finite volume discretization with Newton's method is used to compute steady base flows. The linear stability problem is treated as a GEVP in which the rightmost eigenvalue determines the stability of the base flow. The rightmost eigenvalues associated with the most dangerous eigenfunctions are computed through the use of the shift-invert Arnoldi method. To avoid fine meshing in the regions where the flow variables are changing slowly, a local mesh refinement technique is used in order to increase numerical accuracy with a lower computational cost. The CUBIT mesh generation environment has been used to refine the quadrilateral meshes locally. In order to achieve higher performance on parallel machines, the algebraic systems of equations resulting from the steady problem and the GEVP have been solved by implementing the MUltifrontal Massively Parallel Solver (MUMPS). The proposed method is applied to the linear stability analysis of the flow of an Oldroyd-B fluid past a linear periodic array of circular cylinders in a channel and a linear array of circular half cylinders placed on channel walls. Two different leading eigenfunctions are identified for close and wide cylinder spacing for the periodic array of cylinders. In addition, the numerical method is also extended for the Rolie-Poly model and the numerical results are compared with the experimental results obtained from a Multi-Pass Rheometer (MPR) using the rounded 8:1:8 contraction-expansion slit. The numerical results confirm the instabilities observed in the experiments and identify important parameters within the model, which gives physical insight into the underlying mechanism involved.