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字节,还是随着变量的不同而不同。
附测试代码:
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
mkl手册写的是 double complex for zgetrf, mkl中lapack不支持四精度的实数或复数。complex(16) 和 complex*16 是不一样的,前者实部和虚部都是4精度,占用两个16字节;后者实部和虚部都是双精度,为两个8字节。 不建议现代代码使用 complex*8 或 complex*16 这种形式,对应的改写为 complex(4) 或 complex(8) li913 发表于 2023-12-5 21:32
mkl手册写的是 double complex for zgetrf, mkl中lapack不支持四精度的实数或复数。complex(16) 和 comple ...
感谢您的回复,已解决 fcode 发表于 2023-12-6 08:36
不建议现代代码使用 complex*8 或 complex*16 这种形式,对应的改写为 complex(4) 或 complex(8) ...
是的,我一般都采用后者,感谢回复
页:
[1]