问题描述
我正在尝试解决广义特征值问题(我想要特征值和特征向量):
[A]{x} = lambda[B]{x}
或等效形式(有限元方法):
[M]{x} = (1/w^2)[K]{x}
[M]=[A]
和[K]=[B]
(分别为质量和刚度矩阵)的位置。
我的纯文本文件中同时包含[M]
和[K]
。这些矩阵以COO格式存储,基于0的索引,未排序。
我的文件(下载: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
使用文件中的数据,我可以使用SciPy(Python)解决问题。由于这是一个FEA问题,因此我也已经在Abaqus CAE中解决了这个问题,因此我知道在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)
我的Python代码:
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)),cols))
eigenvalues,eigenvectors = eigs(A=a,M=b,k=10)
val = real(eigenvalues) # to remove "+0j",eigenvalues are real
v = 1.0/val
我的问题是,为什么英特尔MKL无法解决此问题?我尝试修改代码,但总是收到错误或异常:
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
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,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)
! 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,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,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.
我整天都在看这段代码,我无法分辨自己在做什么错。任何帮助将不胜感激。
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)