Fortran Coder

标题: Fortran求特征根问题 [打印本页]

作者: zzgjtx    时间: 2016-9-28 22:28
标题: Fortran求特征根问题
系统: win10
编译:NAG  Fortran Builder & NAG  Fortran
问题描述:再求特征根问题时,首先调用 dpotrf将质量矩阵进行Cholesky factorization,其后将特征根问题利用DSYGST转化为标准特征根问题。
然后利用DSYEVR求解特征根。
[Fortran] 纯文本查看 复制代码

  call dpotrf(uplo,nsum,M_sum,nsum,info1)                    !The Cholesky factorization of the mass matrix
  call dsygst(itype,uplo,nsum,K_sum,nsum,M_sum,nsum,info2)   !Reduces a generalized eigenvalue problem to the standard form

  call dsyevr(jobz,range,uplo,nsum,K_sum,nsum,vl,vu,il,iu,abstol,f_m,Eig_getw,Eig_getz,  &
              nsum,isuppz,Eig_work,lwork,iwork,liwork,info3)
  lwork=min(Lwmax,int(Eig_work(1)))
  liwork=min(Lwmax,iwork(1))  
  call dsyevr(jobz,range,uplo,nsum,K_sum,nsum,vl,vu,il,iu,abstol,f_m,Eig_getw,Eig_getz,  &
              nsum,isuppz,Eig_work,lwork,iwork,liwork,info4)    !Eigenvalue solution

但在求解时发现,总会提示 Execution of the program ended because it received signal [SIGFPE].
经过网络查找,发现是除零异常,signal(SIGFPE, sig_fpe). 但没完全看明白怎么处理。
此外,同样的程序利用IVF进行计算时,没有遇到这个问题。
请问各位有没有什么好的方法解决这个问题。谢谢

作者: pasuka    时间: 2016-9-29 08:03
首先,这个不是程序的问题:
1、LAPACK难道不能求解广义特征值问题吗?
2、有限元质量阵就一定是SPD吗?

作者: zzgjtx    时间: 2016-9-29 10:20
pasuka 发表于 2016-9-29 08:03
首先,这个不是程序的问题:
1、LAPACK难道不能求解广义特征值问题吗?
2、有限元质量阵就一定是SPD吗?

1. 在我的问题里,质量矩阵是正定对称矩阵,我验证过了。
2. 不管ivf还是NAG函数库,都是调用的LAPACK。我没用过你之前说的命令行编译。有什么好的材料可以借鉴么?

作者: pasuka    时间: 2016-9-29 10:38
zzgjtx 发表于 2016-9-29 10:20
1. 在我的问题里,质量矩阵是正定对称矩阵,我验证过了。
2. 不管ivf还是NAG函数库,都是调用的LAPACK。 ...

一步到位改成DSYGVD不是挺好嘛,NAG网站都有现成例子
http://www.nag.com/lapack-ex/node104.html

作者: zzgjtx    时间: 2016-9-29 11:39
pasuka 发表于 2016-9-29 10:38
一步到位改成DSYGVD不是挺好嘛,NAG网站都有现成例子
http://www.nag.com/lapack-ex/node104.html
...

谢谢,我又仔细去看了这个例子。你说的是对的,最近被特征根问题搞的头大了。

问题好多,发现我求解的A*z = λ*B*z这样一个问题,B为实对称正定矩阵,A为实非对称矩阵。
问题又变得复杂了,可能要按照非对称矩阵的解法来求解了。
作者: pasuka    时间: 2016-9-29 21:06
zzgjtx 发表于 2016-9-29 11:39
谢谢,我又仔细去看了这个例子。你说的是对的,最近被特征根问题搞的头大了。

问题好多,发现我求解的A* ...

计算方法和矩阵论补齐的话,再配上戈卢布的矩阵计算,就没有搞不明白的LAPACK程序调用问题




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