Why lsqr is slow?

Hello everyone, I am trying to solve a large linear equations Ax=b in which the size of coefficient matrix A should be in the range 300,000 by 300,000 or even larger. The equation is acquired from the FEM so it is sparse. Because the original PDE is not a 2nd order elliptic equation, so the matrix A is not symmetric. I tried several different built-in functions in Matlab (2012 a) to solve it, \ dosen't work and it took really long time and got a wrong answer. cgs.m dosen't converge for 1000 iterations. If I pre-processed A using ilu.m, I can't wait. lsqr.m can only converge with relative large error 10^(-3) and it took over 10 minutes on my laptop (Win 7, Core i7 processor with 2.1GHz, 6G RAM). But if I used smaller error bond, it cannot converge within 1000 iterations. So any suggestions on this problem? Thanks!

Best, Feixiao