Hi Aldo,

Thank you very much for your reply!

Yes, PARDISO uses palletization and I used call mkl_set_num_threads(numthred) to set the number of threads (numthred) and by setting  msglvl = 0 to print statistical information from which I can confirm how many threads are actually using.  By varying numthred from 1 to 8 using the same A and b matrix,   i got wall time for the first try: 18.5s, 15.1s, 13.6s, 10.4s, 11.6s and for the second try 20.1s, 22.9s, 11.7s, 9.5s, 34.6s.

I appears time deduction as numthred increases (but not pronounced as i expected i even thought is fluctuaion may not comes from the contribution of more threads) if we take into account of the fluctuations of wall time at different runs of the same code. 

My desktop is 8-core and I will further test it on a more powerful computer once I can access it. 

So do you have any suggestions for how many threads is suitable for the problem of this size, or is there any better linear solver  for this size of problem?

Thank you very much and have a nice day!

Best,

Xiang

Hi Aldo,

Than you for your reply!

You are right that solvers are often dedicated to a specific matrix structure and  the matrix is randomly generated non-symmetric sparse matrix  with density 0.6 so it is not block tridational matrix.  I use MATLAB function sprand to do it. 

PARDISO should be good at sparse matrix according to the manual, at least better than DGESV I think. Matlab backslash is actually a very powerful function that take different factorization schemes by analyse the matrix structure first and they also consider multi-threading too. So my goal is actually to see PARDISO faster than DGESV at least and hopefully can approach MATLAB's speed.

I am not sure whether have I made anything wrong , but I did check that the solutions from both DGESV and PARDISO are correct and accurate (L2 norm of the error between the exact solution and the solution is at O(-3)).

Will report here if I have some new findings! thank you again Aldo!

Best,

Xiang

Hi Aldo,

Apologises  for my late response as I was focusing of another code earlier today.

Your suggestions actually point me to a better direction toward understanding the problem I think:

1) You are right that PARDISO conduct steps including analysis, numerical factorization, solve and  which can be controlled by the flag you provide. A lot of physical problems have a constant A but different b and factorization is needed once. The following output of one run of my code gives the time consumed at different steps:

Summary: ( starting phase is reordering, ending phase is solution )

================

Times:

======

Time spent in calculations of symmetric matrix portrait (fulladj): 1.300993 s

Time spent in reordering of the initial matrix (reorder) : 1.639168 s

Time spent in symbolic factorization (symbfct) : 0.355624 s

Time spent in data preparations for factorization (parlist) : 0.000218 s

Time spent in copying matrix to internal data structure (A to LU): 0.000000 s

Time spent in factorization step (numfct) : 3.807952 s

Time spent in direct solver at solve step (solve) : 1.310971 s

Time spent in allocation of internal data structures (malloc) : 6.229885 s

Time spent in additional calculations : 1.135247 s

Total time spent : 15.780058 s

Which basically says if I were to solve the system 10 times, only another Time spent in direct solver at solve step*9 is needed.

Also, Time spent in allocation of internal data structures is interesting to me as I am trying to see whether can I decrease it.

2) My original Matlab code to calculate MATLAB speed is incorrect. As Matlab did something similar  as PARDISO and I generate the code and do the calculation within the same code. That means Matlab already know the structure of the matrix before solving the system which may indicate it may save the time of analyse the matrix. By using an independent code to read the matrix and solve the system, Matalb need  3.x seconds, that is one seconds longer than previous.

3) I will do a thorough reading of the PARDISO manual later this week and hopefully find something useful.

I appreciate your feedback and will report to you once I got some new findings.

Best,

Xiang

Hi Peter,

Thank you very much for you comment and suggestions. Will definitely investigate it further along the lines you mentioned and report back here. Need time for another code and may only be able to give the feedback one or two days later. Thank you again and have a nice day!

Best,

Xiang

Matrix size of 6000 to 20000 is quite small for modern computers, especially when it's real (not complex-valued). For my problems, PARDISO performed much faster than Matlab backslash. Scalability of matrix factorization with PARDISO was close to linear using 1-8 cores.

As Peter has noticed, 6s spent in malloc is very strange.

Hi Vladimir,

Thank you very much for your comments.

I will check the malloc as you and Peter suggested tonight and tomorrow. And also what is the problem size of yours that demonstrate faster speed and good scalability? Or do you have any experience what solver is appropriate for a problem size in my case?

Thank you very much!

Best,

Xiang 

I think that the main problem is in high density (density between 0.6 and 0.8) , so the matrix is almost full. PARDISO (and other sparse libraries) are aimed at max. speed, but also on minimization of fills (keep sparsity of output as low as possible). PARDISO assumes only few nonzero elements per row. This is reason for long allocation and long proccesing time. Matlab is the fastest because it solves the matrix as a dense. Maybe you could try SuperLU or dense LAPACK implementation.

Hi Ivan,

Thank you very much for your comment!

I think you made a very good point here that PARDISO performs well at higher sparsity. When I looked at examples with good performance other people have with PARDISO, they generally have very sparse matrix (Density even less than 0.01 percent ).  For example,

https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/611238

and

https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/509122

On the other hand, there should be  a better way to do the Malloc as suggested by Peter and Vladimir. However, as I was trying code examples of mkl_malloc from mkl manual  with c preprocessor and have not yet got comfortable/familiar enough to do it by myself ahead of PARDISO call.

The problem size can be another possible factor too as suggested by Vladimir, this can also be seen from the above two link which have much larger problem size.

Currently, for my problem the LAPACK routines DGESV (falls into the category of dense solver you mentioned) is faster than PARDISO. I will compare them again after I can do my own Malloc and also a try on SuperLU.

Hi Xiang,

Regarding you question about the size, we have recently tried quite large complex matrices in our paper "Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems". There is a post-print in my profile here: http://www.researchgate.net/publication/291420521_Evaluation_of_parallel_direct_sparse_linear_solvers_in_electromagnetic_geophysical_problems

In particular, it has some numbers and figures for PARDISO applied to matrices with 1-4 millions of complex unknowns.

By "density between 0.6 and 0.8" you mean that 60-80% elements of our matrix are not zeros, right? It is a dense matrix in that case.

Article Evaluation of parallel direct sparse linear solvers in elect...

Hi Vladimir,

Thank you very much for pointing me to your very interesting paper. Very informative and comprehensive comparison of several currently popular linear solvers.

Yes, that is what I mean by density and the density of four matrices in your problem are 0.00003098122, 0.00000592369, 0.00000668426 and 0.00000160117. Besides, they have a much larger DOF. I think these two factors as pointed out also by Ivan are  why they perform and scale so well in your examples.

So currently, I tend to focus more on DGESV other than sparse solver for my type of problem.  I appreciate your time and all these I learned from you (and also from above repliers) !

Best,

Xiang

Hi Xiang,

Then I guess Ivan is right in his comment above. Your main problem is in high matrix density.

Hello,

I have a similar problem as Xiang, that PARDISO library scales badly with the number of cores, sometimes is even slower with growing number of cores. I wonder why is it so? (My matrix is indeed sparse)

In one post Aldo asked is the matrix block-tridiagonal and my matrix is. Maybe it is an unfavorable pattern for the solver? I use PARDISO (5.0.0) because it has a selective inversion algorithm which is claimed to be also parallelized.

Thank you for help,

Alex.

Hi Alexander,

Thank you for your post. It might be the reason of unfavorable pattern for the solve. Another thing you need to pay attention is the size and sparsity of your matrices. As commented by Ivan and Vladimir, these two are also important factors that might affect the performance of PARDISO. Thank you!

Xiang

Thanks, Xiang

My matrix is sparse, I checked. Maybe it is a block+tridiagonal pattern then...

I am also using the original PARDISO library not the one from intel. I wrote a python wrapper for it. So, I thought maybe it is an issue of compilation....

Alex.

Hi Alexander,

Do you mean that numerical factorization scales badly with the number of cores? What is the size and nnz of your typical matrix?

You may try using different reordering. In PARDISO 5.0.0 this is IPARM(2). If IPARM(2) is 0 then the minimum degree algorithm is applied, if IPARM(2) is 2 the solver uses the nested dissection algorithm from the METIS package (default).

Vladimir

Hi, Vladimir

I am solving a linear system (symbolic factorization, numerical factorization, solving), computing the determinant and computing selective inversion. Currently, I do these 3 operations in a row, so I am not quite sure what is the slowest one.

I run on my laptop and also checked on a desktop. I am varying the number of cores by the environment variable OMP_NUM_THREADS. I tried only 3 different values: 1,2 and 4. Then I monitor the actually working cores by looking at the output of htop command. Also, I measure the running time. (By the way, when I  do multiplication of large matrices and check how it scales with the number of working cores, I observe a clear speed up.)

So, the result is basically that increasing the number of working cores does not lead any speed up even to slowing down. In htop sometimes I see that several cores are working, but most of the time only one irrespective of the value in  OMP_NUM_THREADS.

I have tried to set IPARM(2)=0, it has no effect. I have block-tridiagonal matrix with 50000 blocks (on diagonal) and size of the block is 16. So, the matrix is 800000 * 800000. The number of nonzeros is 2649969, and percentage of nonzeros is 0.000414.

Thank you for help,

Alex.

Hi Alexander,

This is strange indeed. Perhaps you have an unfavorable pattern of the matrix as Xiang suggested. You may check if the same code scales well on another matrix with a different pattern.

Vladimir