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/) Powered by Discuz! X3.2