I implemented a FEM solver in MATLAB for Poisson's equation in 3D, using hexahedron and sparse matrix for the Laplacian. I was using the backslash but now I have to use a few iterative methods (GMRES with ilu for example). With $N_x$, $N_y$ and $N_z$ (in general, $N$) being the number of elements (50, 100 and 200 as test cases) I have some problems. For $N=50$ it takes more than $600$ seconds using cgs with incomplete LU preconditioner. For $N=100$ it takes more than 60 minutes and it is still running. I did not set the tolerance or the maximum number of iterations since it does not seem to affect the timing problem. Is this possible or it should not take that long? Is there a way to speed it up?
Here's the code:
[L,U] = ilu(Kb); [u_h,flag,relres,iter,resvec] = cgs(Kb,fb,,,L,U);
where Kb is the Laplacian and fb the right-hand side.