The Helmholtz equation is the simplest possible model for the wave propagation. Perhaps this is the reason, despite denying traditional iterative methods like Krylov sub-space methods, Multigrids, etcetera, numerical solution of the Helmholtz equation has been an interesting and abundant problem to researchers since years. The work in this dissertation is also classified as an attempt to develop fast and robust iterative methods for the solution of the Helmholtz equation. This works is specified for applications in seismic imaging-Geophysics, where usually high frequency are used. Thus we will be targeting large wavenumber Helmholtz problems. The finite difference discretization of the Helmholtz equation with typically given Absorbing (Sommerfeld) boundary conditions gives rise to symmetric, non-Hermitian, indefinite linear systems. Resolution of large wavenumber requires larger number of grid points, thus large linear systems. Many (sparse) direct solvers and hybrid (direct and iterative) solvers have been proposed, but it is quite obvious for very large problems that (sparse) direct solvers have been too much depending upon memory, which makes them less acceptable. Quite a lot of work has been invested in researching iterative solution methods for the Helmholtz equation since many decades. The indefiniteness, which increases with respect to an increase in the wavenumber, poses extra problems for iterative solvers and robust solution of indefinite (large) linear system forms an important research activity. Many iterative techniques like domain decomposition methods, multigrid methods and preconditioners for Krylov subspace methods have been proposed but non of them has been quite robust. For multigrid methods, indefiniteness arises difficulties in having both good smoothing property and constructing appropriate coarse-grid approximations of the problem, which are responsible for further reduction of low frequency errors. Many attempts have been spent in algebraic variants of multigrid methods. Some of them works well with limitation of homogeneity. Most of them fails to show satisfactory convergence. The same holds for Krylov subspace methods. One of the difficulties for Krylov methods is to find a cheap and performing preconditioner for the indefinite Helmholtz equation. An overview of preconditioners, ranging from classical to matrix based, for indefinite Helmholtz linear system has been give in this thesis. A matrix-based complex shifted Laplace preconditioner (CSLP) has been seen as best in the available ones. However, with increasing wavenumbers CSLP shows a slow convergence behavior. We address this issue continuing using CSLP while taking care of its requirement of specific complex shifts. The projection-type preconditioners have been widely investigated by researchers in numerical analysis community. We propose the projection-type deflation preconditioner to tackle the near-singular nodes, which are the cause of the decay the convergence of, this otherwise well performing, CSLP. Like multigrid, this deflation pre-conditioner, named as ADEF1, requires to solve coarse problems at different coarser levels. An optimized algorithm has been tested and proposed suggesting iterative solution of coarse problems at different levels. This finalizes as a multilevel preconditioner. The re-discretization coarsening strategy that we propose and investigate in this thesis is aimed at reducing the memory size and maintaining stencil size. The multilevel Krylov method (MLKM) has also been investigated and compared with its counterpart ADEF1. The rigorous Fourier analysis (RFA) to investigate the convergence of iterative methods forms a separate research theme, which is included in the thesis. We analyse the proposed multilevel preconditioners ADEF1 and MLKM for two-levels. Analysis shows spectral behavior of the preconditioner, which can be taken as favorable for Krylov methods. RFA points out near-singular modes and highlights their contribution in prevailing stagnation. Further the convergence can be enhanced by adapting coarse grid operator at different levels. The proposed preconditioners have been tested on academic as well as the bench mark Marmousi problem. A huge reduction in number of iterations can be noticed. A comparison in the amount of iterations and solve time, specially for three-dimensional problem, shows that the invested work has paid-off. Proposed preconditioners has been uniformly performing for one- to three-dimensions as well as for heterogeneous medium problems.