# Why is Intel MKL failing to solve this eigenproblem while SciPy has no problems (SPARSE_STATUS_ALLOC_FAILED)

0

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.

``````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

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]

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

open(unit=1, file="./coo/mr")
open(unit=2, file="./coo/mc")
open(unit=3, file="./coo/mv")
do i = 1, a_nnz
end do
close(unit=1)
close(unit=2)
close(unit=3)

open(unit=1, file="./coo/kr")
open(unit=2, file="./coo/kc")
open(unit=3, file="./coo/kv")
do i = 1, b_nnz
end do
close(unit=1)
close(unit=2)
close(unit=3)

! 0 to 1 based indexing
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.

python
scipy
fortran
intel-mkl
asked on Stack Overflow Sep 10, 2020 by Carlos