在 Fortran90 中使用 LAPACK、BLAS 计算矩阵的逆

Question

我是使用 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)
       ycoord = (j-1)*mesh(2)
       xd=sin(2.0*3.14*xcoord)
       yd=sin(2.0*3.14*ycoord)
       f((i-1)*(xdiv-1)+j,1)= 8.0*3.1415*3.1415*xd*yd
         enddo
        enddo
    do i=1,(xdiv-1)*(ydiv-1)
     do j=1,(xdiv-1)*(ydiv-1)
      Ainv(i,j)=A(i,j)
     end do
    end do

 call DGETRF((xdiv-1)*(ydiv-1), (xdiv-1)*(ydiv-1), Ainv,  & 
  (xdiv-1)*(ydiv-1), ipiv, info)

 call DGETRI((xdiv-1)*(ydiv-1), Ainv, (xdiv-1)*(ydiv-1), ipiv, &
  work,(xdiv-1)*(ydiv-1), info)

 call DGEMV(N, (xdiv-1)*(ydiv-1), (xdiv-1)*(ydiv-1), 1.d0, Ainv, &
  (xdiv-1)*(ydiv-1), f, 1.d0, 0.d0, u, 1.d0)

  do i=1,(xdiv-1)*(ydiv-1)
   write(*,*) "u", u(i,1)
  enddo
 end program main

基本上，我计算 LU 分解，然后将其取反，然后相乘。请帮我找到错误，并建议一个更好的方法来进行这个计算（如果有的话）。

注意：变量定义/赋值的某些部分可能看起来多余或低效，因为这是更大代码的一部分。我只是提取了它的一部分，专注于矩阵求逆问题。

Answer 1

对于一般平方矩阵的求逆，请查看此函数（使用 LAPACK）

  function inv(A) result(Ainv)
    implicit none
    real,intent(in) :: A(:,:)
    real            :: Ainv(size(A,1),size(A,2))
    real            :: work(size(A,1))            ! work array for LAPACK
    integer         :: n,info,ipiv(size(A,1))     ! pivot indices

    ! Store A in Ainv to prevent it from being overwritten by LAPACK
    Ainv = A
    n = size(A,1)
    ! SGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call SGETRF(n,n,Ainv,n,ipiv,info)
    if (info.ne.0) stop 'Matrix is numerically singular!'
    ! SGETRI computes the inverse of a matrix using the LU factorization
    ! computed by SGETRF.
    call SGETRI(n,Ainv,n,ipiv,work,n,info)
    if (info.ne.0) stop 'Matrix inversion failed!'
end function inv

此外，您可以在此处找到其他有用的 Fortran 子程序和矩阵运算函数。

基本上，这遵循与您相同的过程，除了如某些注释中所述，您传递给 LAPACK 子例程的参数不正确。看看那些并修复它们或只是使用这个功能。

编辑：请注意，这是针对单精度数据类型。您可以轻松地将其调整为双精度。