Fortran Coder

标题: zgetrf函数的使用 [打印本页]

作者: ucas 时间: 2023-12-5 17:34
标题: zgetrf函数的使用
本帖最后由 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

作者: li913 时间: 2023-12-5 21:32
mkl手册写的是 double complex for zgetrf， mkl中lapack不支持四精度的实数或复数。complex(16) 和 complex*16 是不一样的，前者实部和虚部都是4精度，占用两个16字节；后者实部和虚部都是双精度，为两个8字节。

作者: fcode 时间: 2023-12-6 08:36
不建议现代代码使用 complex*8 或 complex*16 这种形式，对应的改写为 complex(4) 或 complex(8)

作者: ucas 时间: 2023-12-6 10:48

li913 发表于 2023-12-5 21:32
mkl手册写的是 double complex for zgetrf， mkl中lapack不支持四精度的实数或复数。complex(16) 和 comple ...

感谢您的回复，已解决

作者: ucas 时间: 2023-12-6 10:50

fcode 发表于 2023-12-6 08:36
不建议现代代码使用 complex*8 或 complex*16 这种形式，对应的改写为 complex(4) 或 complex(8) ...

是的，我一般都采用后者，感谢回复

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)