Fortran Coder

标题: MKL能求解非对称矩阵部分特征值的问题吗? [打印本页]

作者: 木木林    时间: 2019-6-5 00:42
标题: MKL能求解非对称矩阵部分特征值的问题吗?
    自己写了个有限元程序,最后需要求解一个大型(大概1000*1000)非对称稀疏矩阵的部分特征值及其特征向量,如最大的10个特征值与对应的特征向量。最近研究了下  Intel MKL  Fortran 的参考手册,发现里面有两个地方提到了特征值的求解:


1.LAPCK---->LAPACK Least Squares and Eigenvalue Problems Routines  

   这里面能够求解四种特征值问题:对称EVPs,对称正定的广义EVPs,非对称EVPs和广义非对称EVPs。还提供了计算Computationl Routines 和Driver Routines 两种版本。应该用的是QR迭代法。

2.Extened Eigensolver Routines
  这里面提到了一种Feast 算法,支持CSR和BSR等压缩矩阵格式,可以求解最大,最小特征值,和在计算给定范围内的特征值和特征向量,但是似乎只支持Symmetric和Hermitian两种性质的矩阵。

那么现在的问题是:
     LAPACK的特征值求解器不支持CSR和BSR等压缩存储格式,也不支持求解部分特征对;而Extended Eigensolver 又不支持非对称矩阵。那么
     MKL能够求解  非对称  稀疏矩阵的   最大的若干个特征值及其特征向量吗?


作者: Jackdaw    时间: 2019-6-5 13:05
?ggev
Computes the generalized eigenvalues, and the left
and/or right generalized eigenvectors for a pair of
nonsymmetric matrices.


可以求所有的特征值和特征向量,1000阶很快
作者: 木木林    时间: 2019-6-5 23:52
Jackdaw 发表于 2019-6-5 13:05
?ggev
Computes the generalized eigenvalues, and the left
and/or right generalized eigenvectors for a ...

感谢回复。
我是用1D FEM求一个旋转梁的流固耦合问题,实际中我一般只需要前10-20个特征值和特征向量分别表示固有频率和模态形状够了。

ggev是LAPACK里面的driver routine。这个LAPACK里面的eigensolver 用的是QR迭代法,需要用稠密阵存储矩阵,不像Extended eigensolver 一样支持稀疏矩阵存储格式,如CSR。而且后期有限元网格会加密,可能带来效率问题和特征向量的存储开销问题。如10000*10000,光所有特征向量就需要存储一亿个浮点数,太夸张了。

目前比较清楚的情况是:

   目前的MKL里的Extended eigensolver 用的是FEAST 2.0, 所以不支持非对称和非Hermitian矩阵。
   2015年发布的FEAST 3.0是支持非对称矩阵的,网上看到Intel的人2016年就说注意到了,但直到现在也没加到MKL。

我尝试下看能不能让MKL调FEAST 3.0。

作者: pasuka    时间: 2019-6-7 06:58
木木林 发表于 2019-6-5 23:52
感谢回复。
我是用1D FEM求一个旋转梁的流固耦合问题,实际中我一般只需要前10-20个特征值和特征向量分 ...

直接用ARPACK不是挺好的嘛!
传送门
https://github.com/opencollab/arpack-ng

作者: 木木林    时间: 2019-6-10 20:07
pasuka 发表于 2019-6-7 06:58
直接用ARPACK不是挺好的嘛!
传送门
https://github.com/opencollab/arpack-ng

   感谢回复。
   之前对特征值问题完全不了解,又因为用的是windows,只是对MKL熟悉点,所以一直想在MKL的框架里解决问题,同时也是想利用MKL中高度硬件优化的BLAS、SPBLAS和LAPCK库。  

   最近折腾了一段时间后,也了解到了ARPACK,对于非对称、非Hermitian特征值问题总结如下,以供后人学习参考:

    1. 中小型特征值问题
        这类问题可以采用QR迭代法直接求解所有特征值和特征向量,只支持稠密矩阵,主要是LAPACK中的相应subroutine用的就是这类方法。

    2. 大型特征值问题的部分特征值问题。
        2.1 ARPACK
             ARPACK采用Krylov 子空间方法和implicit restart 技术,大概在1996年前后开发完成。
        2.2 Feast
            Feast借鉴量子力学中的围道积分法(contour integration),对由围道(contour)在复平面中所围的区域内进行特征值搜索,可以同时设置多个围道,并行效果较好。Feast 2.0只支持对称和Hermitian问题,Feast 3.0 支持所有矩阵问题(对称,非对称,Hermitian,非Hermitian  , dense  , banded, sparse)。Feas 3.0的开发人员中有intel的人,所以很多BLAS和LAPACK库默认是从MKL中调用的,稀疏特征值问题还用到了MKL中的PARDISO模块,因此编译的时候最好用Intel Fortran 编译器,大概语法为:
    ifort  /Qmkl   /fpp  /O2  /c  /I.  *.f90

具体理论可参考Feast 2.0【1】和3.0【2】的原始论文:
  【1】 Eric Polizzi,2009,A Density Matrix-based Algorithm for Solving Eigenvalue Problems
  【2】JAMES KESTYNy, ERIC POLIZZIz, AND PING TAK PETER TANG,2016,FEAST EIGENSOLVER FOR NON-HERMITIAN PROBLEMS
  

   

  

作者: Jackdaw    时间: 2019-6-14 15:28
木木林 发表于 2019-6-10 20:07
感谢回复。
   之前对特征值问题完全不了解,又因为用的是windows,只是对MKL熟悉点,所以一直想在MKL ...

厉害,我涉及特征值的问题只是为了求谱条件数,用了MKL里GGSV,发现慢,改用了幂法反幂法
作者: pasuka    时间: 2019-7-10 23:05
木木林 发表于 2019-6-10 20:07
感谢回复。
   之前对特征值问题完全不了解,又因为用的是windows,只是对MKL熟悉点,所以一直想在MKL ...

lz的口气好大哦!
1、大型特征值,多大才算大呢?若是大型稠密矩阵呢?
2、Lanczos法、子空间迭代法、Jacob-Davidson法总不能无视吧?
3、FEAST方法的点子很新颖,将子空间迭代法扩展到复数域。可是判断区间内特征根数量有啥好法子嘛?
4、PARDISO是业界标杆,可是SuperLU、UMFPACK、SPQR也蛮好玩的
5、学习有限元开发,为啥不用C++的Deal II
作者: chdlkl    时间: 2020-5-20 12:17
您好,请问压缩存储稀疏对称矩阵后,求矩阵特征值的问题解决了没
作者: asri623    时间: 2020-8-15 00:16
试了一下在用FEAST进行区间特征值求解,发现效率比Lanczos法高,但对在有些问题中经常找不到特征值,搞不懂原因?请指教!




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2