"Reese" wrote in message <ivkasi$9uf$1@newscl01ah.mathworks.com>...
> "Peter" wrote in message <ivipt2$4hm$1@newscl01ah.mathworks.com>...
> > "Reese" wrote in message <ivhc2a$fq6$1@newscl01ah.mathworks.com>...
> > > "Peter" wrote in message <ivg2km$5rv$1@newscl01ah.mathworks.com>...
> > > > Hello,
> > > >
> > > > I am solving a 4th order PDE using finite differences. This amounts to solving a matrix equation where Ax = b. A is large (1e6 x 1e6) and sparse, but not symmetric (even though the FD stencils are symmetric, my geometry is irregular  circular with holes inside).
> > > >
> > > > Currently, at 1e6x1e6, mldivide works ok. I want to push the limits though, because I am not sure the solution has converged enough. I've gone up to 4e6x4e6, but the program then takes a while, ~ 11.5 hrs.
> > > >
> > > > I tried reverting to iterative methods, but I can't seem to precondition the system well enough. If I just take the diagonal (since the matrix is diagonally dominant) gmres seems promising initially, getting to a rel res ~ 10^2, but then sort of stagnates, and takes many runs to get down to the required tolerance. The other methods are also slow. For preconditioning, ILU doesn't work as there are zeros on the diagonal.
> > > >
> > > > Are there are any direct/iterative methods (and preconditioners) anyone thinks might improve the performance?
> > > >
> > > > Thanks,
> > > > Peer
> > >
> > > If you are just interested in performance I would suggest porting the code to fortran or C. If that doesnt tickle your fancy have you checked that there is nothing in the code you can change to speed it up (not algorithm) such as vectorization? If you haven't already run the profiler and see where your bottlenecks are and focus there.
> >
> > Hi Reese,
> > Thanks for the advice. While improvements probably can be made in the body of the code, there will still always be a matrix inversion step (using either a direct method, such as mldivide, or an iterative method) that is very costly for large N. I'm curious though about your first line regarding C/fortran. Are these languages faster for matrix inversion operations as well? My impression was that MATLAB is optimized for these processes above all other languages. So if I download the appropriate libraries, etc., the code might be faster?
> > Thanks
> > Peter
>
> C/fortran aren't faster for matrix operations they are faster overall due to matlab being an interperted language while they are compilied languages. The difference being as a matlab script is run it compilies each line as it reaches it, and then has to execute the command. The algorithm that is then used may be optimized. In C/fortran the code is compilied once into a binary and then runs directly without any wasted time at runtime.
>
> The downside to using C/fortran is that they are more strict on programing style (declaring variables, pointers etc) which matlab handles for you, not to mention a different language to learn. In addition to that aspect you would have to find the algorithms needed for your solution. While traditional tri/pentdiagonal solvers are common it sounds like you will need a more specialized solver which could be difficult to find (I haven't worked with alot of asymmetric solvers yet).
>
> Another thought to speed this up, depending on your resources you could use GPU acceleration to compute some of this. Supported on 2010(a I think but it might be b) and later you can use the Parellel Computing Toolbox(PCT) to use the computers CUDA gpu to do calculations. Another option is Jacket.
>
> Just curious but you said this is a 4th order PDE so what scheme are you using? Some sort of compact scheme to handle BCs?
Hi Reese,
Thanks again. I looked into GPU processing but unfortunately I do not have NVIDIA.
The scheme is rather simple. I have a square grid, with an embedded circle inside. The circle is the physical domain; however there are various holes in the circle that are excluded from the domain. On all boundaries (exterior of circle + boundaries of holes), the value of the function and its normal derivative are specified.
The 4th order PDE is spatial (no time dependence). For the biharmonic operator, I use a 13 point stencil, that approximates the biharmonic operator as a weighted sum of neighboring points. With that said, each point in the square grid falls under two categories: 1) The point is in the domain and obeys some relation given by the finite difference formula 2) The point lies outside the domain, and I fix its value using the boundary conditions (just with linear extrapolation). Then, all these equations get put into a matrix A*x = b and I use mldivide. So basically, for a point in the domain, the corresponding row of the matrix either has 13 nonzero entries, while a point outside the domain has a row with just a 1 on the diagonal (and the lefthandside, b is adjusted so that the point is set to BC) and a 0 elsewhere.
I don't know if that made sense, but if it sounds overly complicated and there is a simpler way to do this, please let me know. I've tested it in some cases and it seems to be working to some degree of accuracy.
Thanks!
Peter
