Tim,
Thanks for all your comments  they've really helped me understand the problem a bit better; and I'll give cholmod+metis a try (I already tried suitesparse, but wasn't using the correct solver).
You're very welcome to add the matrices to your sparse collection; they're probably quite useful as reference, being generic poisson solvers.
Kind regards, and safe travelling
Tom Clark
"Tim Davis" <davis@cise.ufl.edu> wrote in message <gitgrl$lq5$1@fred.mathworks.com>...
> Thanks for the details, Pat. So if you want to use backslash, x=(A)\(b) would be the thing to use. Or x=cholmod(A,b) if I remember the syntax (and be sure to compile CHOLMOD with METIS).
>
> I'm travelling just now, so I might not be able to check the matrix, but if I can I'll let you know what CHOLMOD+METIS can do for it.
>
> Comment to Tom: large 3D problems are posterchildren for PCG and the like. It would be interesting to see what CHOLMOD+METIS can do, though. I'll give it a crack.
> And may I add them to my collection (http://www.cise.ufl.edu/research/sparse)?
>
> "Pat Quillen" <pquillen@mathworks.com> wrote in message <gipfhl$f60$1@fred.mathworks.com>...
> > "Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <giob75$eei$1@fred.mathworks.com>...
> > <SNIPfor brevity>
> > > Can anyone help me to understand why the slight change in sparsity pattern makes so much difference please?
> > >
> > > If anyone can suggest a different method to backslash, which uses less memory during the solution, I'd be very grateful... I really don't want to have to implement a grid based method (trying to keep my code as simple as possible).
> > <SNIPEND>
> >
> > Hi Tom.
> >
> > For starters, it appears that your matrix is not s.p.d. Instead, it looks to be negative definite. So, if you try to solve Ax = b via backslash, MATLAB will not be using the Cholesky factorization as you might hope, but instead, it will use a symmetric indefinite solver which is just a shade slower. If instead, you try to solve Ax = b, then you can use the MATLAB function pcg with A, as that matrix is s.p.d. pcg will use *much* less memory than backslash and you can get just as much accuracy as you like. Try this:
> >
> > % Get coefficient matrix and rhs from poisson3d. Call them A, b respectively.
> > % No need to copy and negate A; use anonymous function to negate A*x.
> > >> f = @(x) (A*x);
> > >> tic; x = pcg(f,b,1e11,1000); toc % See help pcg for more info.
> > pcg converged at iteration 265 to a solution with relative residual 9.8e12
> > Elapsed time is 10.863882 seconds.
> > >> norm(A*xb)
> >
> > ans =
> >
> > 2.5098e09
> >
> > >> tic; y = (A)\(b); toc % Compare to backslash...
> > Elapsed time is 63.965740 seconds.
> >
> > This is exactly the matrix which pcg loves... If you use a decent preconditioner, you can get the solution more quickly, but you have to construct the preconditioner. Under normal circumstances, I'd recommend the zerofill incomplete Cholesky decomposition, but the current one in MATLAB is not as optimal as it could be. Instead, I point you to the ilu function (incomplete LU factorization). You can use this instead. This is sort of cheating because in theory, one should only use s.p.d. preconditioners with pcg, but in practice pcg tends to not notice.
> >
> > >> tic; [L,U] = ilu(A,struct('type','nofill')); toc
> > Elapsed time is 0.059273 seconds.
> > >> % Since we're going to apply pcg to A, need L (or U, but not both).
> > >> tic; x = pcg(f,b,1e11,1000,L,U); toc
> > pcg converged at iteration 94 to a solution with relative residual 8.5e12
> > Elapsed time is 6.227631 seconds.
> > >> % Notice, with this preconditioner, pcg took about 1/3 the iterations and
> > >> % about 2/3 the time. Try the modified ilu:
> > >> tic; [L,U] = ilu(A,struct('type','nofill','milu','row')); toc
> > Elapsed time is 0.194986 seconds.
> > >> tic; x = pcg(f,b,1e11,1000,L, U); toc
> > pcg converged at iteration 64 to a solution with relative residual 7.7e12
> > Elapsed time is 4.116877 seconds.
> >
> > The MILU is a bit more expensive to cook up, but it pays off as pcg with the MILU preconditioner uses a lot fewer iterations and about 40% of the time. Of course, this is about 1/16 the time of \, but with slightly less accuracy. Try asking pcg for a smaller residual norm:
> > >> tic; x = pcg(f,b,1e13,1000,L, U); toc
> > pcg stopped at iteration 91 without converging to the desired tolerance 1e13
> > because the method stagnated.
> > The iterate returned (number 77) has relative residual 2.2e13
> > Elapsed time is 5.829291 seconds.
> > >> norm(A*x  b)
> >
> > ans =
> >
> > 5.5067e11
> >
> > >> norm(A*y  b) % y = (A)\(b);
> >
> > ans =
> >
> > 4.2224e11
> >
> > These two are equivalently good in terms of residual norm, and x was much cheaper in terms of time and memory to come by. The bit about stagnation basically tells us that this is as good as pcg can do. If you would like to see the residual norm histories, try this:
> > >> [x,fl,rr,it,rv] = pcg(f,b,1e13,1000,L, U);
> > >> semilogy(rv);
> >
> > Now, to actually attempt an answer to your question, I believe that AMD, the default method for reordering sparse matrices within MATLAB doesn't give as high quality orderings for matrices arising from 3d geometries as it does for those arising from 2d geometries. Essentially the change in sparsity pattern caused much more fill to occur in the factors of your matrix. Tim Davis is the guy to really answer your question, and I'm sure he'll be along soon. Also, part of the issue is that the symmetric indefinite solver uses more memory than the Cholesky solver. Solving Ax = b required about 400MB less on my machine:
>
> > >> tic; A\b; toc
> > Elapsed time is 118.416031 seconds.
> > >> % Approximate peak memory usage ~2.05 GB
> > >> tic; (A)\(b); toc
> > Elapsed time is 64.046205 seconds.
> > >> % Approximate peak memory usage ~1.65 GB
> >
> > Note that all of these timings are from MATLAB 7.7 (R2008b) on a 64bit Linux machine.
> >
> > Take care.
> > Pat.
