Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7014 Discussions

## how to choose a solver for symmetric positive definite linear problem?

Novice
1,114 Views

Hello,

I am working on a project using Finite Element Method and want to use intel MKL to boost the performance of solving a linear problem `A*x = b`.

- `A` is a real symmetric positive definite matrix, highly sparse, and with NxN dimensions. N is ~4M.

- `b` is a Nx1 vector.

- `x` is the unknown.

In my case, `A` is fixed, and I only need to do factorization once. And `b` changes every time when I call the solver.

I have found that in MKL, the pardiso solver can be used to solve symmetric positive definite matrix (if you set mtype = 2), but there is no C examples about it. Moreover, I have found examples for solving symmetric positive definite system using RCI solver.

Why is there no example for solving symmetric positive definite matrix using Pardiso?

If I use RCI, is it possible to separate the factorization part and the solving part, as in my case `A` is always fixed and only b keeps changing?

If yes, will I expect big performance difference between using pardiso and RCI?

Regards,

Fang

1 Solution
Employee
1,102 Views

Hello Fang,

1) Please have a look at the example pardiso_sym_c.c. It has mtype = -2 but you can easily change it to mtype = 2 and it will work.

2) RCI interface is available for iterative methods while PARDISO is a direct sparse solver. In your case, unless you know a great preconditioner for your problem (possibly, if you know weel your application) or you don't need to solve your systems with low accuracy and since you have many rhs, then it makes more sense to use the direct sparse solver (PARDISO).
So, in general you need to decide whether you want to solve your system with an iterative or a direct method (but I suggest PARDISO). These are two fundamentally different approaches and performance comparison is tricky as it depends on the details.

3) Note (in case you haven't), that using phase parameter you can exactly do what you need with PARDISO API. Do the reordering, symbolic and numerical factorizaton once (phases 12) and then for each new rhs do only phase = 33. Also note, that f you know your rhs all in advance, you can call PARDISO with multiple rhs at once.

Best,
Kirill

4 Replies
Employee
1,103 Views

Hello Fang,

1) Please have a look at the example pardiso_sym_c.c. It has mtype = -2 but you can easily change it to mtype = 2 and it will work.

2) RCI interface is available for iterative methods while PARDISO is a direct sparse solver. In your case, unless you know a great preconditioner for your problem (possibly, if you know weel your application) or you don't need to solve your systems with low accuracy and since you have many rhs, then it makes more sense to use the direct sparse solver (PARDISO).
So, in general you need to decide whether you want to solve your system with an iterative or a direct method (but I suggest PARDISO). These are two fundamentally different approaches and performance comparison is tricky as it depends on the details.

3) Note (in case you haven't), that using phase parameter you can exactly do what you need with PARDISO API. Do the reordering, symbolic and numerical factorizaton once (phases 12) and then for each new rhs do only phase = 33. Also note, that f you know your rhs all in advance, you can call PARDISO with multiple rhs at once.

Best,
Kirill

Novice
1,086 Views

Hi Kirill,

Thank you very much for the explanations. Suppose the multiple rhs has `M` columns, is it faster than running single rhs `M` times, using an intel CPU, for example, an intel CORE i7 8th generation?

Regards,

Fang

Moderator
1,061 Views

Yes, it would be faster, but the level of solution stage scalability depends on a number of factors. You may try and get us know the results.

Moderator
944 Views

The issue is closing and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.