LAPACK DGETRF + DGETRI не работает - PullRequest
Новая
       75

LAPACK DGETRF + DGETRI не работает

3 голосов
/ 06 февраля 2020

Я пытаюсь инвертировать матрицу с помощью подпрограмм 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, и они привели к одним и тем же результатам. Я также пытался изменить некоторые параметры, такие как длина рабочего массива, и все, что я делал, приводило к этому.

1 Ответ

5 голосов
/ 06 февраля 2020

Вы определяете массивы одинарной точности, но вызываете процедуры двойной точности. Позвоните 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

, скомпилированный со следующим make-файлом: 

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

Примечание: как отмечено в комментариях, я использовал реализацию LAPACK от MKL для удобства. Независимо от реализации точность подпрограмм и переменных должны совпадать. То, что происходит, когда они этого не делают, может зависеть от реализации, но, скорее всего, будет бесполезным во всех случаях.

...