Why is Intel MKL failing to solve this eigenproblem while SciPy has no problems

[csouto]

First of all, this is my first time on the Intel Community forums, so I apologize if this is posted in the wrong section.

---

I am trying to solve a generalized eigenvalue problem (I want both the eigenvalues and eigenvectors):

[A]{x} = lambda[B]{x}

Or equivalent (Finite Element Method):

[M]{x} = (1/w^2)[K]{x}

Where [M]=[A] and [K]=[B] (mass and stiffness matrices, respectively).

I have both [M] and [K] in plain text files. These matrices are stored in COO format, 0-based indexing, not sorted.

My files (attached, or download from Google Drive

kc: where K columns are stored
kr: where K rows are stored
kv: where K values are stored
mc, mr, mv: analogous to K, but for the M matrix

Using the data from the files, I can solve the problem with SciPy (Python). Since this is a FEA problem, I've also solved it in Abaqus CAE, so I know for a fact that both eigenvalues and eigenvectors are correctly calculated in SciPy:

MODE ABAQUS (w^2) SciPy (w^2)
1 +6.07235E+06 +6.50440E+06
2 +2.28087E+08 +2.44463E+08
3 +1.67357E+09 +1.79572E+09
4 +3.36973E+09 +3.37316E+09
5 +5.88655E+09 +6.32761E+09
(value of 1/lambda = w^2)

My Python code:

from scipy.sparse import coo_matrix
from scipy.sparse.linalg import eigs
from numpy import real

rows = 610
cols = 610
a_nnz = 4628
b_nnz = 9266

# read a
with open('./coo/mr') as f: a_row_indx = [int(i) for i in f]
with open('./coo/mc') as f: a_col_indx = [int(i) for i in f]
with open('./coo/mv') as f: a_val_indx = [float(i) for i in f]

# read b
with open('./coo/kr') as f: b_row_indx = [int(i) for i in f]
with open('./coo/kc') as f: b_col_indx = [int(i) for i in f]
with open('./coo/kv') as f: b_val_indx = [float(i) for i in f]

a = coo_matrix((a_val_indx, (a_row_indx, a_col_indx)), shape=(rows, cols))
b = coo_matrix((b_val_indx, (b_row_indx, b_col_indx)), shape=(rows, cols))

eigenvalues, eigenvectors = eigs(A=a, M=b, k=10)

val = real(eigenvalues) # to remove "+0j", eigenvalues are real
v = 1.0/val

My question is, why is Intel MKL failing to solve this problem? I've tried tinkering with the code, but I always get errors or exceptions:

pm(3) = 1 ! -> Exception thrown at 0x00007FFCE77556EC (mkl_core.dll) in Console1.exe: 0xC0000005: Access violation accessing location 0x0000000000000000.
pm(3) = 2 ! -> STAT = 2 SPARSE_STATUS_ALLOC_FAILED

My Fortran code (Visual Studio):

include 'mkl_solvers_ee'
include 'mkl_spblas'

program main

	use mkl_solvers_ee
	use mkl_spblas
	use, intrinsic :: iso_c_binding , only : c_int, c_double
	
	implicit none
	
	type(sparse_matrix_t) :: a
	type(sparse_matrix_t) :: b
	integer(c_int), parameter :: rows = 610
	integer(c_int), parameter :: cols = 610
	integer(c_int), parameter :: a_nnz = 4628
	integer(c_int), parameter :: b_nnz = 9266
	integer(c_int) :: a_row_indx(a_nnz)
	integer(c_int) :: a_col_indx(a_nnz)
	real(c_double) :: a_values(a_nnz)
	integer(c_int) :: b_row_indx(b_nnz)
	integer(c_int) :: b_col_indx(b_nnz)
	real(c_double) :: b_values(b_nnz)
	integer :: stat
	
	character, parameter :: which = 'S'
	integer(c_int) :: pm(128)
	type(matrix_descr), parameter :: descra = matrix_descr(type = sparse_matrix_type_general, mode = sparse_fill_mode_upper, diag = sparse_diag_non_unit)
	type(matrix_descr), parameter :: descrb = matrix_descr(type = sparse_matrix_type_general, mode = sparse_fill_mode_upper, diag = sparse_diag_non_unit)
	integer(c_int), parameter :: k0 = 10
	integer(c_int) :: k
	real(c_double) :: e(k0), ee(k0)
	real(c_double) :: x(k0, cols)
	real(c_double) :: res(k0)
	
	type(sparse_matrix_t) :: acsr
	type(sparse_matrix_t) :: bcsr
	
	integer :: i
	
	! read a
	open(unit=1, file="./coo/mr")
	open(unit=2, file="./coo/mc")
	open(unit=3, file="./coo/mv")
	do i = 1, a_nnz
		read(1, *) a_row_indx(i)
		read(2, *) a_col_indx(i)
		read(3, *) a_values(i)
	end do
	close(unit=1)
	close(unit=2)
	close(unit=3)
	
	! read b
	open(unit=1, file="./coo/kr")
	open(unit=2, file="./coo/kc")
	open(unit=3, file="./coo/kv")
	do i = 1, b_nnz
		read(1, *) b_row_indx(i)
		read(2, *) b_col_indx(i)
		read(3, *) b_values(i)
	end do
	close(unit=1)
	close(unit=2)
	close(unit=3)
	
	a_row_indx(:) = a_row_indx(:) + 1
	a_col_indx(:) = a_col_indx(:) + 1

	b_row_indx(:) = b_row_indx(:) + 1
	b_col_indx(:) = b_col_indx(:) + 1
	
	stat = mkl_sparse_d_create_coo(a, sparse_index_base_one, rows, cols, a_nnz, a_row_indx, a_col_indx, a_values)
	stat = mkl_sparse_d_create_coo(b, sparse_index_base_one, rows, cols, b_nnz, b_row_indx, b_col_indx, b_values)
	stat = mkl_sparse_convert_csr(a, sparse_operation_non_transpose, acsr)
	stat = mkl_sparse_convert_csr(b, sparse_operation_non_transpose, bcsr)
	stat = mkl_sparse_ee_init(pm)
	
	!pm(1) = 0
	!pm(2) = 6
	!pm(3) = 2
	!pm(4) = 512
	!pm(5) = 60 ! 10000 
	!pm(6) = 512
	!pm(7) = 0
	!pm(8) = 0
	!pm(9) = 0
	
	!pm(3) = 1 ! -> Exception thrown at 0x00007FFCE77556EC (mkl_core.dll) in Console1.exe: 0xC0000005: Access violation accessing location 0x0000000000000000.
	!pm(3) = 2 ! -> STAT = 2 SPARSE_STATUS_ALLOC_FAILED
	
	stat = mkl_sparse_d_gv(which, pm, acsr, descra, bcsr, descrb, k0, k, e, x, res)
	
	ee(:) = 1.0 / e(:)
	
end program main
	
!	0	SPARSE_STATUS_SUCCESS			The operation was successful.
!	1	SPARSE_STATUS_NOT_INITIALIZED	The routine encountered an empty handle or matrix array.
!	2	SPARSE_STATUS_ALLOC_FAILED		Internal memory allocation failed.
!	3	SPARSE_STATUS_INVALID_VALUE		The input parameters contain an invalid value.
!	4	SPARSE_STATUS_EXECUTION_FAILED	Execution failed.
!	5	SPARSE_STATUS_INTERNAL_ERROR	An error in algorithm implementation occurred.
!	6	SPARSE_STATUS_NOT_SUPPORTED		The requested operation is not supported.

I have wasted all day looking at this code and I can't tell what I am doing wrong. Any help is greatly appreciated.

[Gennady_F_Intel]

Which version of MKL did you try?

[csouto]

I am using Intel® Parallel Studio XE for Visual Studio. In IntelSWTools directory it seems I have two versions, compilers_and_libraries_2020 and compilers_and_libraries_2020.0.166.

[Gennady_F_Intel]

It might be the reals issue with the latest mkl and we will escalate the issue. In the case, if you need to have the official support, then please submit the request to the Intel Online Service Center.