Fortran Coder

查看: 3658|回复: 5
打印 上一主题 下一主题

[数学库] lapack zheev 出现数组越界错误如何更改

[复制链接]

50

帖子

18

主题

0

精华

熟手

F 币
276 元
贡献
167 点
跳转到指定楼层
楼主
发表于 2022-3-7 23:50:51 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
在跑今天下面的这个代码的时候,算本征值和本征向量,但是报错,问了一下应该是数组越界错误,应该如何改呢
[Fortran] 纯文本查看 复制代码
rogram workfile
    Implicit none
    integer , parameter :: q = 1000
    real , parameter :: t1 = 3.12 , t2 = 0.29 , t3 = -0.0103
    complex , parameter :: i = (0, 1) , t0 = (0,0)
    real , parameter :: PI = acos(-1.0) , a0 = 1.42e-10
    real :: b1(2), b2(2) , a1(2) , a2(2) , m1 , m2
!    complex H(6, 6)
    real , allocatable :: k(:,:,:)
    complex :: f1 , f2
    integer :: x , y
! Output of the Hamiltonian matrix    
!    integer::i1

! Eigenvalue and Eigenvector of the Hamiltonian matrix
    CHARACTER JOBZ , UPLO
    integer :: INFO=0
    integer , parameter :: N=6, LDA=N ! , LWORK=3*N
    real*8 RWORK(3*N-2), W(N)
    complex*16 H(LDA,N), WORK(3*N)

    allocate(k(2,q-1,q-1))
    b1 = [ 2*PI/a0*sqrt(3.0) , 2*PI/a0 ]
    b2 = [ 2*PI/a0*sqrt(3.0) , -2*PI/a0 ]
    a1 = [ a0*sqrt(3.0)/2 , a0/2 ]
    a2 = [ a0*sqrt(3.0)/2 , -a0/2 ]
    Forall(x=1:size(k,2), y=1:size(k,3))
            k(:,x,y) = (b1*x+b2*y)/q
    End Forall
    m1 = dot_product(a1 , k(:,1,2))
    m2 = dot_product(a2 , k(:,1,2))
    f1 = 1+exp(i*m1)+exp(i*m2)
    f2 = 1+exp(-i*m1)+exp(-i*m2)
    H = reshape([complex::&
            t0, t1*f1, t0, t0, t0, t0,&
            t1*f2, t0, t0, t2, t0, t3,&
            t0, t0, t0, t1*f2, t0 ,t0,&
            t0, t2, t1*f1,t0, t0, t2, &
            t0, t0, t0, t0, t0, t1*f1,&
            t0, t3, t0, t0, t1*f2, t0],[6 , 6])

! Output of the Hamiltonian matrix
!    do i1=1,6
!        write(*,*)H(i1,:)
!    end do
!print *, H(:,:)

! Eigenvalue and eigenvector of the Hamiltonian matrix
call myzheev(N,W,H)
call ZHEEV('V', 'U', N, H(6,6), N, W, WORK, 3*N, RWORK, INFO)
print *,W

  contains
 subroutine myzheev(N,W,H)
    implicit none
    integer::N, INFO=0
    real*8::W(N),RWORK(3*N-2)
    complex*16::H(N,N),WORK(3*N)
    call ZHEEV('V', 'U', N, H, N, W, WORK, 3*N, RWORK, INFO)
  endsubroutine

End Program workfile

出现了以下报错:
[Fortran] 纯文本查看 复制代码
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7f8e4ef818a0 in ???
#1  0x7f8e4ef80a45 in ???
#2  0x7f8e4ec9451f in ???
	at ./signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0
#3  0x5626f8c2f1e3 in ???
段错误 (核心已转储)

最后还是把结果返还给了我,如果是数组越界的话应该如何修改呢,由于这里用了lapack的包,所以在print的时候也不太好print
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

168

帖子

2

主题

1

精华

大师

Vim

F 币
1021 元
贡献
486 点

规矩勋章

沙发
发表于 2022-3-8 10:18:22 | 只看该作者
你调用了两次zheev,50行那个没有必要,传入了H(6,6)会直接访问到未知的内存区域,造成段错误

50

帖子

18

主题

0

精华

熟手

F 币
276 元
贡献
167 点
板凳
 楼主| 发表于 2022-3-8 14:21:05 | 只看该作者
Transpose 发表于 2022-3-8 10:18
你调用了两次zheev,50行那个没有必要,传入了H(6,6)会直接访问到未知的内存区域,造成段错误 ...

好的,十分感谢!

50

帖子

18

主题

0

精华

熟手

F 币
276 元
贡献
167 点
地板
 楼主| 发表于 2022-3-8 14:25:34 | 只看该作者
Transpose 发表于 2022-3-8 10:18
你调用了两次zheev,50行那个没有必要,传入了H(6,6)会直接访问到未知的内存区域,造成段错误 ...

另外我还有个问题,就是我在输出本征向量的时候,输出的是复数,但是加上一个real后,如real(H(:,:)),输出的便是实型,这两种数据类型是如何转换的呢?谢谢!

168

帖子

2

主题

1

精华

大师

Vim

F 币
1021 元
贡献
486 点

规矩勋章

5#
发表于 2022-3-8 16:55:14 | 只看该作者
real如果参数是复数的话,返回的是复数的实部。最好注明对应的kind值 例如 real(H(:,:),kind=4)

50

帖子

18

主题

0

精华

熟手

F 币
276 元
贡献
167 点
6#
 楼主| 发表于 2022-3-8 19:06:40 | 只看该作者
Transpose 发表于 2022-3-8 16:55
real如果参数是复数的话,返回的是复数的实部。最好注明对应的kind值 例如 real(H(:,:),kind=4) ...

好的谢谢!
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-12-24 03:12

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表