fortran - LAPACK DGETRF+DGETRI 失败
问题描述
我试图通过 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 库二进制文件,它们得到了相同的结果。我还尝试更改某些参数,例如工作数组的长度,而我所做的一切都导致了这一点。
解决方案
您定义单精度数组,但调用双精度例程。打电话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
注意:如评论中所述,出于方便,我使用了 MKL 的 LAPACK 实现。无论实现如何,例程和变量的精度都应该匹配。当他们不这样做时会发生什么可能取决于实现,但很可能在所有情况下都是垃圾。
推荐阅读
- php - 如何从另一个php文件的php文件中的会话中获取值
- c# - 视图可以选择性地填充出现在 ASP.NET mvc 中视图之外的部分吗?
- python - 如何在 Windows 10 的 Pycharm 上使用 ThunderSVM(GPU 模式)
- assembly - 多条 nop 指令不会始终比单条 nop 指令花费更长的时间
- javascript - 为什么是vue invalid with v-bind:src?
I want to dynamically switch html content, so I used vue-loader src to import, but v-bind:src doesn't take effect at all.
<template src="./app.html"></template>
- typescript - 如果可选参数被指定为特定值,则返回类型“从不”
- android - 您必须提供 layout_height 属性
- c# - 我正在尝试为使用按钮做出的选择编写一个过滤器
- r - ggplot中的斜体希腊字母
- php - 使用 Angular JS 和 PHP 上传文件