Subsurface flows, occurring in groundwater movement and production of hydrocarbons in the petroleum industry, are affected by the heterogeneity of the medium varying over large scales. In this thesis, we have used state of the art multiscale methods to solve one such flow model, influenced by high contrast permeability field. The focus is on the elliptic pressure equation which is solved on both fine and coarse scale for comparison purposes. Shared memory parallelism has been achieved for generating basis functions, which is computationally the most expensive portion of the multiscale implementation. Parallel systematic spectral enrichment using the GMsFEM (Generalized Multiscale Finite Element Method) is the key feature of the current work and has been compared with the MsFEM (Multiscale Finite Element Method). Efficiency of algorithmic implementation has been first tabulated for a two-dimensional finite element code using the MATLAB parallel computing toolbox and also has been given a more generalized form using a three-dimensional finite element code written in OpenMP (Open Multiprocessing) and C++. The timing comparison shows a significant decline in the execution time for the algorithms. It indicates that a higher level of enrichment and desired accuracy is achievable for large scale problems. Computational time gain and fewer memory requirements are two key features achieved in this work. Distributed parallel computing can further be implemented to achieve mass parallelism through which one can solve large problems accurately and efficiently when compared to benchmark fine scale solutions where global system solver, memory requirements and execution time become significant issues.