Windows环境下声振有限元非对称稀疏广义特征值问题的Python求解求助

我最近在处理一个声振有限元的广义特征值问题，卡壳好一阵了——我的刚度矩阵（stiff）和质量矩阵（mass）都是非对称的稀疏矩阵，但Scipy里的scipy.sparse.linalg.eigs和eigsh都明确只支持对称/厄米特矩阵的场景，这直接导致我的Python代码跑不起来。

先说说我之前的尝试：我一开始写了这两段代码，想求解最小的20个特征值（nev=20）：

import scipy
import numpy as np
import pickle

if __name__ == '__main__':
    with open('../k_full.pickle', 'rb') as f:
        print('loading matrix K')
        k_global = pickle.load(f)
    with open('../m_full.pickle', 'rb') as f:
        print('loading matrix M')
        m_global = pickle.load(f)

    # 第一种尝试：指定求最小模的特征值
    eigvalues, eigvectors = scipy.sparse.linalg.eigs(A=k_global, k=20, M=m_global, which='SM')
    # 第二种尝试：以0为中心搜索特征值
    eigvalues, eigvectors = scipy.sparse.linalg.eigs(A=k_global, k=20, M=m_global, sigma=0)

结果直接报了收敛错误：

scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence: ARPACK error -1: No convergence (321 iterations, 4/20 eigenvectors converged)

后来查Scipy的文档才搞明白：eigsh本身就只支持实对称或复厄米特矩阵；eigs虽然是通用特征值求解器，但当带M矩阵（广义特征值问题场景）时，也要求M是对称/厄米特矩阵，而我的问题里矩阵显然不满足这个条件。

对比之下，同样的问题我用Matlab是能完美解决的，Matlab代码如下：

k = load('k_global.mat');
m = load('m_global.mat');
[V,D] = eigs(k.array, m.array, 20,0);
D = diag(D)
natural_frequency = sqrt(D)/(2*pi)

返回的结果完全符合声振问题的预期，既有实特征值，也有虚部极小的复特征值：

0.00000000000000 + 0.000356591992067188i
0.000668911454165071 + 0.00000000000000i
0.00000000000000 + 0.000973128785222090i
0.00222975851379527 + 0.00000000000000i
0.00246434216130016 + 0.00000000000000i
0.00000000000000 + 0.00372951940564144i
8.06883871646537 + 0.00000000000000i
64.7482150103242 + 0.00000000000000i
234.453670549319 + 0.00000000000000i
268.154072409059 + 0.00000000000000i
312.537263749716 + 0.00000000000000i
356.103849178590 + 0.00000000000000i
389.038117338274 + 0.00000000000000i
412.048267727649 + 0.00000000000000i
473.729345964820 + 0.00000000000000i
2996.35112385098 + 0.00000000000000i
3240.96766107255 + 0.00000000000000i
4186.42444133727 + 0.00000000000000i
4585.99172192305 + 0.00000000000000i
4794.52737053778 + 0.00000000000000i

补充下我的问题细节：这是典型的声振耦合问题，对应的矩阵是分块形式的耦合结构；我用的是Windows系统，刚度和质量矩阵都是稀疏矩阵，它们的非零元素分布结构也已经整理完成。

现在就想请教下，有没有办法在Windows环境下，用Python处理这种非对称的稀疏广义特征值问题，得到和Matlab一致的结果？

备注：内容来源于stack exchange，提问作者SleepyGary