I am trying to solve the system of linear equations: $AX=B$. For this currently I am using Intel MKL Pardiso solver. It works well when the order of $A$ is around $13500\times13500$ and below. Above that the solver gives wrong answer (perhaps I am giving wrong values of parameters required to be set before calling solver). Dimensions of $A$ and $B$ are needed to be as high as ($15\cdot 10^6\times15\cdot 10^6$) and ($15\cdot 10^6 \times 6$), respectively. My program is in Fortran fixed form and I am using Intel Visual Fortran Compiler.

Reason to use Sparse linear solver:

a) The $A$ matrix is primarily sparse. I have attached morphology of my $A$ matrix for visualization, where yellow pixels represents non-zero values.

b) I was previously using DGELSD LAPACK subroutine, but it is taking so much time. Sparse solvers are about 180 times faster than non-sparse solvers. Also, memory consumption with Pardiso solver is very less as it works on CSR format of the matrix.

My question is: how to set parameters for Pardiso solver, or is there any other solver with less presetting requirements that I can use.

Presently I working on Windows 10. System specification:

Intel Xeon E series processor 32 cores, 64 GB RAM.

Image information:

$A$ matrix morphology, with a total of $135\times135$ elements.

Yellow pixels: 666 in number. Matrix dimensions are $135\times135$ (easiest case). NOTE THAT MATRIX STRUCTURE REMAINS EXACTLY IN THE SAME FASHION WITH INCREASE IN ITS DIMENSIONS, also the A matrix is non symmetric. enter image description here


  • 1
    $\begingroup$ A couple of questions right away: 1) why use DGELSD to solve AX=B? This one is for the least-squares solution, and I doubt that is what you need. 2) Do you have to use a sparse direct solver, or you are open to iterative solver options? 3) It would be helpful if you list you sample PARDISO timing report, to understand what options are you trying, and which steps take most of the time. $\endgroup$
    – Anton Menshov
    Apr 9, 2017 at 11:53
  • 3
    $\begingroup$ Specifically, what is your question? Is there some aspect of the Pardiso documentation that you do not understand? Since you are getting correct answers for smaller matrices, it seems likely you have a bug in your program that you will have to track down. If you search on this site you will find many discussions about different direct and iterative sparse solvers. $\endgroup$ Apr 9, 2017 at 11:53
  • $\begingroup$ In what language should the routine be written? Must it be in Fortran? $\endgroup$
    – DanielRch
    Apr 9, 2017 at 17:55
  • $\begingroup$ DanielRch>>> It can be in Fortran or any other compile level language................ currently my code is in Fortran $\endgroup$ Apr 10, 2017 at 6:47

1 Answer 1


A matrix of size 15M x 15M is likely too big for a (sparse) direct solver on a single machine -- it is going to take too much time and memory. If you wanted to use a direct solver, you could try parallel sparse direct solvers such as MUMPS or SuperLU-dist, both of which are conveniently called via PETSc.

There is also the option of using iterative solvers. This is particularly true if you have knowledge where these linear systems come from, and can come up with a good preconditioner. My thoughts on how to choose solvers and preconditioners can -- in significant detail -- be found in lectures 34-38 here: http://www.math.colostate.edu/~bangerth/videos.html .


Not the answer you're looking for? Browse other questions tagged or ask your own question.