zgetrf函数的使用

ucas · 发表于 2023-12-5 17:34:45

本帖最后由 ucas 于 2023-12-6 21:24 编辑

在使用lapack库zgetrf求解行列式的时候，发现对矩阵a的数据类型定义会影响结果，而且16字节反而错误，只有8字节是正确的。
但是官方的介绍是定义为complex(16)，初学Fortran没多久，请求一下
附代码如下：
定义的三维复数矩阵a，其中det1是解析解，即按照行展开求解的，是正确的参考解。det2是调用zgetrf库函数对角元相乘得到的，即测试解。
当a定义为complex(16)时候，输出结果为：
det1= (24.000000000000000,18.000000000000000)
det2= (135.00000000000000,-135.00000000000000)

当a定义为complex(8)时候，输出结果才一致。

请教一下为什么会出现这种情况，以及普遍的定义是如何？例如所有类型都定义为8字节，还是随着变量的不同而不同。

附测试代码：

[Fortran] 纯文本查看 复制代码

program aaa

    implicit none
    integer, parameter :: n = 3   
    integer :: info,i,sig   
    integer :: ipv(n)            
    complex(16) :: a(n,n)     
    complex(8) :: det1 ,det2   

! 矩阵赋值
    a(1,1)=(1.0,1.0)
    a(1,2)=(2.0,2.0)
    a(1,3)=(0.0,3.0)
    a(2,1)=(3.0,4.0)
    a(2,2)=(5.0,5.0)
    a(2,3)=(5.0,6.0)
    a(3,1)=(7.0,7.0)
    a(3,2)=(8.0,8.0)
    a(3,3)=(9.0,9.0)

!解析解
det1=a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))

! 符号及行列式初始值
sig=1
det2=(1,0)

! 使用 LAPACK 函数求解行列式
    call zgetrf(n, n, a, n, ipiv, info)

do i = 1, 3
    if (ipiv(i) /= i) then
        sig=-1*sig
    end if
    det2=det2*a(i,i)
end do


    write(*,*)  "det1=", det1
    write(*,*)  "det2=", det2*sig


end program aaa