Efficient memory usage for preconditioned conjugate gradient and interapp

I am solving Laplace's Equation on a very large (400x400x400 voxel^3) grid using an implicit finite-difference scheme (matrix inversion). The size of the problem and symmetry of the matrix (symmetric positive definite) makes the problem attractive to solve using conjugate gradient iterative solver method. The benefit of iterative methods typically require far less memory to solve the matrix equation than direct methods. In order to improve convergence of the solver, I am using an Incomplete Cholesky preconditioner. Convergence is very slow when using the standard conjugate gradient solver without preconditioner. However, the memory usage sky rockets when I use the matrix preconditioner. Looking at the pcg and iterapp (located inside pcg) source code, it appears Matlab invokes the mldivide (a direct solver) operation when applying the preconditioner. Since the point of using the conjugate gradient method is to reduce memory requirements for solving the problem, why does Matlab invoke a direct solver operation in the process that makes the memory requirement skyrocket? Is there a more efficient conjugate gradient solver out there that makes more efficient use of memory than Matlab's pcg function?