gfortran 错误：lapack 中的 zgesvd答案

【问题标题】：gfortran error: zgesvd in lapackgfortran 错误：lapack 中的 zgesvd
【发布时间】：2014-05-17 04:00:33
【问题描述】：

我正在对方阵 A 进行 svd 分解，其中 A=U S Vdag，在 fortran 代码中，该行显示为

lwork = -1

call zgesvd( 'A', 'A', A%d, A%d, A%m, A%d, S, U%m, U%d, Vdag%m, Vdag%d,work, lwork, rwork, info )

lwork = int(work(1));deallocate(work); allocate(work(lwork))

call zgesvd( 'A', 'A', A%d, A%d, A%m, A%d, S, U%m, U%d, Vdag%m, Vdag%d,work, lwork, rwork, info )

当我用 gfortran 编译时，它没有错误或警告就完成了。但是，当我运行程序时，它会显示错误消息：

" ** 在进入 ZGESVD 参数编号 11 时具有非法值"

我不知道出了什么问题。

供参考，参数的定义：

type cmatrix
   integer(4) d
   complex(8), allocatable :: m(:,:)
end type

type (cmatrix) A,U,Vdag

allocate(A%m(dim,dim),U%m(dim,dim),Vdag%m(dim,dim))

A%d = dim; U%m = dim; Vdag%d = dim

real(8) S(dim)

提前致谢！小宇

附言需要说明的是，这样的程序用ifort编译时运行流畅，但是gfortran给出如上图的运行时错误

--- 问题解决了！

似乎问题在于 ifortran 和 gfortran 如何分配内存。我在代码中定义了 USV 类型：

type USV
   integer is_alloc
   type (cmatrix) U,V
   real(8), allocatable :: S(:)
end USV

初始化时

type(USV) Test_usv(:)

allocate(Test_usv(3)),

使用 intel fortran 编译器时 is_alloc 的值为 0，而对于 gfortran 则为任意数字。我需要使用这个值作为分配 U V 矩阵的标准：

if (is_alloc.eq.0) then

   allocate(U%m(dim,dim))

end if

【问题讨论】：

第 11 个参数是 Vdag%d。我建议在每次调用之前打印它的值。尝试使用 gfortran 调试选项：-fimplicit-none -Wall -Wline-truncation -Wcharacter-truncation -Wsurprising -Waliasing -Wimplicit-interface -Wunused-parameter -fcheck=all -fbacktrace
你能告诉我们S的定义吗？
谢谢大家！ @M.S.B.：第 11 个参数确实无效。我发现这是一个问题，因为使用 intel fortran 和 gfortran 分配内存的初始化方案不同。

标签： fortran lapack blas

【解决方案1】：

根本问题不是 ifort 和 gfortran 之间的区别。您初始化变量的方法是无效的 Fortran。除非您使用声明、赋值语句等初始化变量，否则它的值是未定义的。解决此问题的一种方法是向类型定义添加默认初始化：

type USV
   integer is_alloc = 0
   type (cmatrix) U,V
   real(8), allocatable :: S(:)
end USV

另一种方法是不自己跟踪分配状态，而是依赖 Fortran 为此目的提供的内在函数：

if (.NOT. allocated (U%m) ) allocate(U%m(dim,dim))

附：最佳实践是不依赖于种类的特定数值。 kind 值是任意的，不一定是类型的字节数。有些编译器使用字节数，有些则不使用。指定字节数和具有可移植代码的一种方法是使用 ISO Fortran 环境提供的类型：

use, intrinsic :: ISO_FORTRAN_ENV

integer(int32) :: d
real(real64), allocatable :: S(:)

类型以位数命名。可用类型的列表在 gfortran 手册的“内部模块”一章中。

【讨论】：