  The Legendre spectral Galerkin method of self-adjoint second order elliptic
equations usually results in a linear system with a dense and ill-conditioned
coefficient matrix. In this paper, the linear system is solved by a
preconditioned conjugate gradient (PCG) method where the preconditioner $M$ is
constructed by approximating the variable coefficients with a ($T$+1)-term
Legendre series in each direction to a desired accuracy. A feature of the
proposed PCG method is that the iteration step increases slightly with the size
of the resulting matrix when reaching a certain approximation accuracy. The
efficiency of the method lies in that the system with the preconditioner $M$ is
approximately solved by a one-step iterative method based on the ILU(0)
factorization. The ILU(0) factorization of $M\in
\mathbb{R}^{(N-1)^d\times(N-1)^d}$ can be computed using $\mathcal{O}(T^{2d}
N^d)$ operations, and the number of nonzeros in the factorization factors is of
$\mathcal{O}(T^{d} N^d)$, $d=1,2,3$. To further speed up the PCG method, an
algorithm is developed for fast matrix-vector multiplications by the resulting
matrix of Legendre-Galerkin spectral discretization, without the need to
explicitly form it. The complexity of the fast matrix-vector multiplications is
of $\mathcal{O}(N^d (\log N)^2)$. As a result, the PCG method has a
$\mathcal{O}(N^d (\log N)^2)$ total complexity for a $d$ dimensional domain
with $(N-1)^d$ unknows, $d=1,2,3$. Numerical examples are given to demonstrate
the efficiency of proposed preconditioners and the algorithm for fast
matrix-vector multiplications.
