LAPACK DGETRF + DGETRI失败

Question

我正在尝试通过LAPACK的DGETRF和DGETRI例程求逆矩阵，但是下面的代码：

program Tester

    !use LapackMatrixOps
    use MathematicalResources

    implicit none

    real :: B(2, 2), A(2, 2), WORK(2)

    integer :: i, j, SIZ, IPIV(2), INFO

    SIZ=2

    do i=1, SIZ
        do j=1, SIZ
            !if(i==j) then
            !    A(i, j)=1
            !else
            !    A(i, j)=0
            !end if
            A(i, j)=rand(0)/2+0.5
            B(i, j)=0
        end do
    end do

    B=A

    call PrintMatrix(A)

    call dgetrf(size(B, 1), size(B, 2), B, size(B, 1), IPIV, INFO)
    print *, "========="
    call PrintMatrix(B)
    print *, IPIV
    call dgetri(size(B, 1), B, size(B, 1), IPIV, WORK, size(WORK), INFO);

    print *, "========="

    call PrintMatrix(B)

end program Tester

编译：

CC=gfortran
CFLAGS=-O2 -W -g
LFLAGS=-L/usr/lib/lapack -llapack -lblas

TARGET=Tester

all: install clean

.PHONY: install
install:
    $(CC) *.f90 $(CFLAGS) $(LFLAGS) -o $(TARGET)

.PHONY: test
test:
    $(CC) *.f90 $(CFLAGS) $(LFLAGS) -o $(TARGET)
    #==============TEST=============#
    ./$(TARGET)

.PHONY: clean
clean:
    rm -rf ./*.o

似乎失败了。当我运行make test时，我得到：

gfortran *.f90 -O2 -W -g -L/usr/lib/lapack -llapack -lblas -o Tester
#==============TEST=============#
./Tester
  0.500003815      0.565768838    
  0.877802610      0.729325056    
 =========
  0.500003815       1.89445039E+28
  0.877802610       1.57469535    
           1           2
 =========
              NaN   6.86847806E-20
              NaN  -1.28451300    
pedro@pedro-300E5M-300E5L:~/Área de Trabalho/AED/openpypm2$ 

我做错了什么？我已经尝试了自定义编译的LAPACK库和ATLAS库二进制文件，它们都得到了相同的结果。我还尝试更改某些参数，例如工作数组的长度，而我所做的一切一直导致这一点。

Answer 1

您定义了单精度数组，但是您调用了双精度例程。分别呼叫sgetrf和sgetri。另外，如果您要使用双精度，请使用以下方法声明变量：

integer, parameter :: dp = kind(1.d0)
real(kind=dp) :: B(2,2), A(2,2), work(2)

无论如何，使用单精度，以下代码：

Program Tester

    implicit none
    real :: B(2, 2), A(2, 2), WORK(2)
    integer :: i, j, SIZ, IPIV(2), INFO

    SIZ=2

    do i=1, SIZ
        do j=1, SIZ
            A(i, j)=rand(0)/2+0.5
            B(i, j)=0
        end do
    end do

    B=A

    print *, A
    call sgetrf(size(B, 1), size(B, 2), B, size(B, 1), IPIV, INFO)
    print *, "========="
    print *, B
    print *, IPIV
    call sgetri(size(B, 1), B, size(B, 1), IPIV, WORK, size(WORK), INFO);
    print *, "========="
    print *, B

end program Tester

使用以下makefile编译：

FC =gfortran

MKL_LIB    = -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -m64 -I${MKLROOT}/include

all: test.f90
        $(FC) -o test.exe $^ $(MKL_LIB) -I$(MKLROOT)/include

产生此输出：

  0.500003815      0.877802610      0.565768838      0.729325056    
 =========
  0.877802610      0.569608510      0.729325056      0.150339082    
           2           2
 =========
  -5.52652836       6.65163040       4.28716612      -3.78882527   

N.B .：如评论中所述，出于方便，我使用了MKL的LAPACK实现。不管实现方式如何，例程的精度和变量都应匹配。每当他们不执行操作时，会发生什么情况都可能取决于实现，但在所有情况下都很可能会成为垃圾。