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

[wang__bin]

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; 

[Kirill_V_Intel]

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]

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]

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]

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