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.
$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.