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