Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Finite Differences- Large Sparse Linear Systems

Subject: Finite Differences- Large Sparse Linear Systems

From: Peter

Date: 11 Jul, 2011 23:57:10

Message: 1 of 9

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, ~ 1-1.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

Subject: Finite Differences- Large Sparse Linear Systems

From: Nasser M. Abbasi

Date: 12 Jul, 2011 04:03:54

Message: 2 of 9

On 7/11/2011 4:57 PM, Peter wrote:
> Hello,
>
>
> Are there are any direct/iterative methods (and preconditioners) anyone thinks might improve the performance?
>
> Thanks,
> Peer

Try multigrid to see if it works for you?

--Nasser

Subject: Finite Differences- Large Sparse Linear Systems

From: Reese

Date: 12 Jul, 2011 11:44:10

Message: 3 of 9

"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, ~ 1-1.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.

Subject: Finite Differences- Large Sparse Linear Systems

From: Peter

Date: 13 Jul, 2011 00:46:26

Message: 4 of 9

"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, ~ 1-1.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

Subject: Finite Differences- Large Sparse Linear Systems

From: Reese

Date: 13 Jul, 2011 14:42:26

Message: 5 of 9

"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, ~ 1-1.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/pent-diagonal 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?

Subject: Finite Differences- Large Sparse Linear Systems

From: Peter

Date: 15 Jul, 2011 19:53:08

Message: 6 of 9

"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, ~ 1-1.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/pent-diagonal 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

Subject: Finite Differences- Large Sparse Linear Systems

From: Piero Lanucara

Date: 16 Jul, 2011 22:22:08

Message: 7 of 9

Are you sure that Matlab is able to solve sparse matrices on GPU? I don't think so....maybe jacket yes...I'm not sure

Piero

> > 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

Subject: Finite Differences- Large Sparse Linear Systems

From: John Melonakos

Date: 17 Jul, 2011 04:04:10

Message: 8 of 9

"Piero Lanucara" <lanucara@caspur.it> wrote in message <ivt2ug$4iu$1@newscl01ah.mathworks.com>...
> Are you sure that Matlab is able to solve sparse matrices on GPU? I don't think so....maybe jacket yes...I'm not sure
>
> Piero

Jacket does contain a Sparse Linear Algebra Library, the only sparse library available for GPU computing in MATLAB. Learn more here, http://wiki.accelereyes.com/wiki/index.php/Jacket_SLA

Learn generally about how Jacket compares, here: http://accelereyes.com/compare

Subject: Finite Differences- Large Sparse Linear Systems

From: Piero Lanucara

Date: 17 Jul, 2011 07:45:09

Message: 9 of 9

that's great news!
do you have some number on Fermi?
thank's
Piero


"John Melonakos" wrote in message <ivtmvq$m62$1@newscl01ah.mathworks.com>...
> "Piero Lanucara" <lanucara@caspur.it> wrote in message <ivt2ug$4iu$1@newscl01ah.mathworks.com>...
> > Are you sure that Matlab is able to solve sparse matrices on GPU? I don't think so....maybe jacket yes...I'm not sure
> >
> > Piero
>
> Jacket does contain a Sparse Linear Algebra Library, the only sparse library available for GPU computing in MATLAB. Learn more here, http://wiki.accelereyes.com/wiki/index.php/Jacket_SLA
>
> Learn generally about how Jacket compares, here: http://accelereyes.com/compare

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us