Solving sparse linear system with eigen/Intel MKL

I want to solve the Ax=b equation where A is a sparse matrix (1,964,568 x 1,964,568 nnz=75256446) and b is also sparse (1,964,568 x 1,964,568 nnz= 25354926) by using Eigen library in C++.

At first I was trying to use the Eigen Sparse LU to solve my problem, after a few hours I ran out of memory (I have 128GB RAM). After this I included the INTEL MKL library with Pardiso solver. Even with this I cant solve my problem. Maybe anyone has some tips to solve my problem?

#define EIGEN_USE_MKL_ALL

#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
#include <iostream>
#include <Eigen/OrderingMethods>
#include <Eigen/PardisoSupport>


typedef Eigen::SparseMatrix<double>SpMat; 
typedef Eigen::COLAMDOrdering<int>Order;

int main()
{

    SpMat A;
    SpMat B;
    SpMat X;

    Eigen::loadMarket(A, "MatK.mtx");
    Eigen::loadMarket(B, "MatM.txt");

    A.makeCompressed();
    B.makeCompressed();

    Eigen::PardisoLU<SpMat>solver;

    solver.analyzePattern(A);

    solver.factorize(A);

    X = solver.solve(B);


}

I can compile my code and run it. I just need better performance and less memory.

Since I haven't found any method to solve this. I tried to change my RHS matrix, because the huge amount of columns was a big problem. I saw in my algorithm that my RHS gets multiplied by another matrix, which reduces the columns to 10. With this it was possible to solve with Intel Pardiso LDLT solver instead of LU. The iterative solver took too many iterations to converge.