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

Question

我不熟悉使用 LAPACK/BLAS。我想计算一个方程的解：

AU=F

我想知道这部分代码的逻辑错误是什么。我使用求解器输入大小为 ((xdiv-1)(ydiv-1),(xdiv-1)(ydiv-1)) 的矩阵 A。然后，随后求解方程。 U=inverse(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 子例程和矩阵运算函数here。

基本上，这与您执行的过程相同，只是如某些 cmets 中所述，您传递给 LAPACK 子例程的参数不正确。看看这些并修复它们或仅使用此功能。

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