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

csouto — Fri, 11 Sep 2020 15:10:12 GMT

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

In order to impose Dirichlet boundary conditions, $[M]$ is singular in my case.

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

EDIT

Code correction:

character, parameter :: which = 'L'

Since I want the lowest w^2, i.e, the largest eigenvalues (lambda = 1/w^2).

pm(3) = 1 ! -> Exception thrown at 0x00007FFCE77556EC (mkl_core.dll) in Console1.exe: 0xC0000005: Access violation accessing location 0x0000000000000000. pm(3) = 2 ! -> Stack trace terminated abnormally. forrtl: severe (157): Message not found - Unhandled exception at 0x00007FF92E64ED79 (ntdll.dll) in Console1.exe: 0xC0000374: A heap has been corrupted (parameters: 0x00007FF92E6B77F0).

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

Gennady_F_Intel — Mon, 14 Sep 2020 04:52:17 GMT

here is a duplicated thread.

please check the same one: https://community.intel.com/t5/Intel-oneAPI-Math-Kernel-Library/Why-is-Intel-MKL-failing-to-solve-this-eigenproblem-while-SciPy/m-p/1208958#M30068

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

Gennady_F_Intel — Mon, 14 Sep 2020 04:53:14 GMT

This issue has been closed we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only

topic Re:Why is Intel MKL failing to solve this eigenproblem while SciPy has no problems in Intel® oneAPI Math Kernel Library

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

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

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