A multiscale spectral generalized finite element method (MS-GFEM) is presented for the solution of large two and three dimensional stress analysis problems inside heterogeneous media. It can be employed to solve problems too large to be solved directly with FE techniques and is designed for implementation on massively parallel machines. The method is multiscale in nature and uses an optimal family of spectrally defined local basis functions over a coarse grid. It is proved that the method has an exponential rate of convergence. To fix ideas we describe its implementation for a two dimensional plane strain problem inside a fiber reinforced composite. Here fibers are separated by a minimum distance however no special assumption on the fiber configuration such as periodicity or ergodicity is made. The implementation of MS-GFEM delivers the discrete solution operator using the same order of operations as the number of fibers inside the computational domain. This implementation is optimal in that the number of operations for solution is of the same order as the input data for the problem. The size of the MS-GFEM matrix used to represent the discrete inverse operator is controlled by the scale of the coarse grid and the convergence rate of the spectral basis and can be of order far less than the number of fibers. This strategy is general and can be applied to the solution of very large FE systems associated with the discrete solution of elliptic PDE.