Using Pardiso64 in Pardiso.jl for Large Sparse Matrices

[del2341]

I am trying to use Pardiso for sparse matrix diagonalization, but my matrix has more non-zero elements than Int32 can handle. To resolve this, I need to use Pardiso64. However, when I set the solver parameter ps.iparm[35] = 1, I still get an error related to the number of non-zero elements exceeding Int32.

I have confirmed that I am calling the 64-bit MKL version.

LinearAlgebra.BlasInt = Int64, but MklInt remains Int32.

The issue seems to come from Pardiso.jl, where the following condition is checked:

if LinearAlgebra.BLAS.vendor() === :mkl && LinearAlgebra.BlasInt == Int64

However, due to a version change, LinearAlgebra.BLAS.vendor() is now returning :lbt instead of :mkl, even though I am using MKL with 64-bit support.

Attempted Fixes:

I tried manually overriding this by commenting out the above check in Pardiso.jl and defining:

const MklInt = Int64
const PARDISO_FUNC = :pardiso_64

Unfortunately, I still encountered the following error:

[2362860] signal (11.1): Segmentation fault

Request for Help

Is there a correct way to enforce Pardiso64 in Pardiso.jl?

How should I properly handle MklInt to ensure compatibility with large sparse matrices?

Could the segmentation fault be caused by incorrect linking, memory issues, or an internal bug in Pardiso.jl?

Any guidance or suggestions would be greatly appreciated!

[mipals]

I don't have any experience with this, but maybe I could have a look. Do you have a MWE?

In either case the line LinearAlgebra.BLAS.vendor() === :mkl should definitely be changed as the vendor is always :lbt on newer version on Julia that uses libblastrampoline (that is at least how I understand it).

[del2341]

Sorry for the delay. Below is a MWE for testing. It runs perfectly as long as the number of nonzero elements is less than 2^{32} . However, if the count exceeds 2^{32} , an error occurs: “ERROR: InexactError: trunc(Int32, 2147483649)”

using LinearAlgebra
using MKL
using SparseArrays
using Pardiso
using Arpack

function create_tridiagonal_hamiltonian(n)
    diagonal = (2.0 + 0.1im) * ones(ComplexF64, n)
    off_diagonal = (-1.0 + 0.2im) * ones(ComplexF64, n-1)
    return spdiagm(-1 => off_diagonal, 0 => diagonal, 1 => off_diagonal)
end

function calculate_nonzeros_tridiagonal(n)
    return n + 2*(n-1)
end

function calculate_required_dimension(target_nonzeros)
    required_n = ceil(Int, (target_nonzeros + 2)/3)
    return required_n
end

println("2^32 = ", 2^32)

target_nonzeros = 2^32
required_dimension = calculate_required_dimension(target_nonzeros)
println("Required dimension to exceed 2^32 non-zero elements: ", required_dimension)

test_dimension = 1441655766
nonzeros_count = calculate_nonzeros_tridiagonal(test_dimension)
println("\nUsing test dimension n = ", test_dimension)
println("Expected non-zero elements: ", nonzeros_count)


H = create_tridiagonal_hamiltonian(test_dimension)

ps = MKLPardisoSolver()
set_matrixtype!(ps, Pardiso.COMPLEX_HERM_INDEF)
pardisoinit(ps)
fix_iparm!(ps, :N)
H_pardiso = get_matrix(ps, H, :N)
b = rand(ComplexF64, size(H, 1))
set_phase!(ps, Pardiso.ANALYSIS)
pardiso(ps, H_pardiso, b)
set_phase!(ps, Pardiso.NUM_FACT)
pardiso(ps, H_pardiso, b)