Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- [Pardiso] high residual for symmetric indefinite system after schur complement

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

wang__bin

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-15-2020
11:30 PM

126 Views

[Pardiso] high residual for symmetric indefinite system after schur complement

I'm trying to use Pardiso to solve a Stokes problem using weak galerkin finite element method.

All **matrix and code are attached** and tested at MKL 2019.5.281, Visual Studio 2019, WIN10

The linear system is a saddle point system which can be solved in Full sysem or Schur complement system. Both linear system are symmetric indefinite.

//2nd order FEM with 384 tets

//**System1. Full system,** where u={u0,ut,v0,vt,w0,wt} and u0,v0,w0 is not shared by neighboring element

[K G * {u = {b_u

G' 0] p} b_p}**system size** = 4080 x 4080

//**System2. Schur complement** **,** where u_bar={ut,vt,wt} as u0,v0,w0 is eliminated

[A B * {u_bar = {b_u

B' C] p} b_p}**system size** = 2640 x 2640

left: Full matrix right: Schur complement matrix

Obviously, System2 can save a lot of memory for a large problme.

**Q: I can use Pardiso to solve the System1 very well (residual = 1e-16). **

** But When I solve the System2, it returns the very high residual (1e-2 -> 1e4). **

** Any suggestions to solve this symmetric indefinite system?**

I already tried the solution of setting following iparm but it is not helps

iparm[9] = 8;

iparm[10] = 1;

iparm[12] = 1;

Link Copied

4 Replies

Kirill_V_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-28-2020
04:47 PM

126 Views

Hello Wang,

I've got a question: isn't your system for the Schur complement which constitutes an equation for the pressure is singular (at least having the constant vector of all ones in the kernel)? If it is, then the solution is defined up to an arbitrary vector from the kernel. So, when you're comparing the obtained solution with the reference solution, you need to filter this component out.

Specifically, if you know that the system is singular and the kernel is the constant vector, then you should subtract from each of the solution vectors their average value and only then compare.

If you don't do the filtering, there might be very well a very large difference between the reference and the obtained solution.

Also, the difference between the solutions is the error, not the residual. It looks to me that your code example measures the error, (x - x_ref) and not the residual which would be Ax - b.

I hope some of this helps.

Best,

Kirill

wang__bin

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-31-2020
09:15 AM

126 Views

Hi Voronin,

I may made an mistake in measuring the residual. But I believe the solution of Pardiso with symmetric solver (mtype=-2 and iparm[9]=8) is wrong.

I can get the right solution by switching it to a unsymmetric solver (mtype=13 and iparm[9]=13) but obviously this will increase the memory usage a lot for a large problem.

Do you have any suggestion to solve this problem using mtype=-2 and iparm[9]=8?

Thanks,

Bin

Voronin, Kirill (Intel) wrote:Hello Wang,

I've got a question: isn't your system for the Schur complement which constitutes an equation for the pressure is singular (at least having the constant vector of all ones in the kernel)? If it is, then the solution is defined up to an arbitrary vector from the kernel. So, when you're comparing the obtained solution with the reference solution, you need to filter this component out.

Specifically, if you know that the system is singular and the kernel is the constant vector, then you should subtract from each of the solution vectors their average value and only then compare.

If you don't do the filtering, there might be very well a very large difference between the reference and the obtained solution.

Also, the difference between the solutions is the error, not the residual. It looks to me that your code example measures the error, (x - x_ref) and not the residual which would be Ax - b.

I hope some of this helps.

Best,

Kirill

Kirill_V_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-31-2020
10:26 AM

126 Views

Hi Bin,

Could you confirm/deny my statement that the system for the Schur complement is singular? If it is singular, I've guessed already what can be the root cause.

Also, how do you switch to the non-symmetric matrix type? For mtype=-2 you need to pass only one triangle of the matrix to the solver. See the doc page, which says "For symmetric matrices, the solver needs only the upper triangular part of the system". So do you use different input data when running with mtype=-2 and mtype=11?

Best,

Kirill

wang__bin

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-01-2020
07:55 PM

126 Views

Hi Voronin,

1. I just check the systems. Both of the system are not singular.

Full system: condition number = 2.9214e+05

Static condensation system: condition number = 1.7266e+05

2. I switched to non-symmetric matrix type in my own code. The testing code uploaded only support symmetric type and the upper triangular part of linear system is offered in my first post.

Do you need the full system for testing purpose?

Thanks,

Bin

Voronin, Kirill (Intel) wrote:Hi Bin,

Could you confirm/deny my statement that the system for the Schur complement is singular? If it is singular, I've guessed already what can be the root cause.

Also, how do you switch to the non-symmetric matrix type? For mtype=-2 you need to pass only one triangle of the matrix to the solver. See the doc page, which says "For symmetric matrices, the solver needs only the upper triangular part of the system". So do you use different input data when running with mtype=-2 and mtype=11?

Best,

Kirill

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.