小编Sch*_*der的帖子

我是使用 LAPACK/BLAS 的新手。我想计算一个方程的解：

AU=F

我想知道这部分代码的逻辑错误是什么。我使用求解器输入大小为 ((xdiv-1) (ydiv-1),(xdiv-1) (ydiv-1)) 的矩阵 A。然后，随后求解方程。U=逆（A）* f。

其中 U 和 f 大小相同。(u((xdiv-1) (ydiv-1),1),f((xdiv-1) (ydiv-1),1))。执行矩阵求逆时出现分段错误错误。

这是我的代码：

program main
double precision, allocatable :: A(:,:)
double precision, allocatable :: u(:,:), f(:,:)
double precision mesh(2), dx, dy
integer  xdiv, ydiv
 xdiv=55
 ydiv=55
 mesh(1)=.001
 mesh(2)=.001
 dx=mesh(1)
 dy=mesh(2)

allocate (A((xdiv-1)*(ydiv-1),(xdiv-1)*(ydiv-1))) 
allocate (Ainv((xdiv-1)*(ydiv-1),(xdiv-1)*(ydiv-1)))
allocate (u((xdiv-1)*(ydiv-1),1),f((xdiv-1)*(ydiv-1),1))

         do i =1,(xdiv-1)*(ydiv-1)
         A(i,i)=-2.d0*(1.d0/(dx**2)+1.d0/(dy**2))
         enddo
         do i=1,(xdiv-2)
         do j=1,(ydiv-1)
         A(i+(j-1)*(xdiv-1),i+(j-1)*(xdiv-1)+1)=1.d0/(dx**2)
         A(i+(j-1)*(xdiv-1)+1,i+(j-1)*(xdiv-1))=1.d0/(dx**2)
         enddo
         enddo
         do i=1,(xdiv-1)
         do j=1,(ydiv-2)
         A(i+(j-1)*(xdiv-1),i+(j)*(xdiv-1))=1.d0/(dy**2)
         A(i+(j)*(xdiv-1),i+(j-1)*(xdiv))=1.d0/(dy**2)
         enddo
         enddo

        do i=1,(xdiv-1)
         do j=1,(ydiv-1)
       xcoord = (i-1)*mesh(1) …