PARDISO - huge memory requirement, slow calculation

[jirina]

I implemented PARDISO solver to solve a linear system of equations with a sparse matrix. I tested it using a model that  yields 1 million equations, each of them containing typically 7 nonzero elements.

All PARDISO settings in iparm are default, matrix type is 11 (real and nonsymmetric), one core is used for solving the system, phase is 13. These settings allocate 15 GB of memory during the phase 13, it takes about 9 min to solve the system.

My colleague tested the very same model in OpenFOAM and its matrix solver (I do not know what exactly they use) needed 1 GB of memory and it only took 1 minute to solve the system. His computer has the same CPU and 1 core was used for the solution.

My questions are:

1. Would the solution be significantly faster and significantly more memory efficient if the matrix was symmetrical?
2. Is there any way in MKL to find out whether the matrix is symmetrical and positive definite or indefinite?
3. Does it make sense to try different sparse solvers than PARDISO?

Thank you in advance for any suggestions and tips.

[mecej4]

As far as I can gather from http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2008/TimBehrens/tibeh-report-fin.pdf , OpenFoam uses iterative methods to solve linear equations. If good trial solutions/preconditioners are available, iterative methods can be faster than methods that use matrix factorization, such as Pardiso.

Regarding your remaining questions: it is quite easy to extract just the upper (or lower) triangle of your matrix, and call Pardiso with that triangle, specifying that the matrix is symmetric, and attempting a solution. You may similarly test whether the matrix is positive-definite. Bear in mind that this type of testing tends to be expensive, so it is better to do such testing using smaller matrices.

Checking for symmetry can be done using a few lines of code. I suspect that if your equations came from discretizing the Navier-Stokes equations, your equations will not be symmetric. In fact, they are not even linear, but you may be solving nonlinear equations by solving a sequence of linear equations.

Note, however, that before choosing a method or attempting to solve the equations, it is almost required that you spend some time acquainting yourself with the characteristics of the equation systems that you wish to solve.

[jirina]

I have already investigated the structure of the matrix in the past and doing it again now confirmed it is symmetrical. You are right it comes from discretizing Navier-Stokes equations, which results in a system of linear equations with a symmetrical matrix; NS equations I am solving contain only the diffusion term.

I succeeded creating the upper triangular matrix and calling PARDISO with it; it took two runs to find out the matrix is not positive definite, but indefinite. Everything works well then, results show that in an example with 500 000 equations,

• Upper triangular matrix: solution time 75.56 s, memory 2.6 GB
• Full matrix: time 100.13 s, memory 4.5 GB

Regarding solution time and memory requirements, they are similar no matter whether I use the full matrix or just the upper triangular one.

Finally, I checked with my colleague that OpenFOAM really uses iterative solvers, which means I will be looking at them in MKL as well. I need to be able to calculate systems with millions of equations, and direct solvers are clearly not suitable for that (well, I could still try solving the system in subdomains of the whole domain, but this starts to be a numerical problem instead of programming).

I checked iterative sparse system solvers in MKL and it seems that the best option for systems with a real, indefinite matrix is FGMRES. Would you please comment on this if you are familiar with sparse matrix solvers in MKL?

PS: By the way, I like your profile picture, it is cute. :-)

[mecej4]

If you are solving problems where the "inertia terms" in the N-S equations are negligible, i.e., low Reynolds number flow, you may be able to use the Poisson solver routines that are provided in MKL. Is the problem steady-state or time-dependent?

When using iterative methods, sometimes it is advantageous to use more than one grid. A solution for the coarse grid can be obtained inexpensively, and used with an interpolation scheme to generate a trial solution for a finer grid. Search for "multigrid", if you have not already done so. There has been an immense amount of work published on such problems and techniques, but I have been out of touch with that since many years.