主页 / coding / 问题 / 79198541
Accepted
TobiR
Asked: 2024-11-18 10:21:23 +0800 CST

将非连续数组的指针传递给 Fortran 中的 Lapack95 例程是否安全？

  • 772

我想对syevd非连续数组调用 Lapack95 子程序，如下所示：

real :: mat(15000, 15000), vec(15000)
mat=1.d0
associate(eig_vects=>mat(:10000,:10000),eig_vals=>vec(:10000))
   call syevd(eig_vects,eig_evals,'V')
end associate

这样做安全吗？之前，我在将不连续的数组切片传递给子例程时遇到了问题（无法分配该子例程创建的临时数组）。使用指针eig_vects作为参数时是否会出现类似的问题？

arrays

1 个回答

  1. Best Answer

    Lapack95 不是 Lapack 的独立实现，而是经典 Lapack 的包装器，具有更现代的接口。特别是，Lapack95 的虚拟数组参数是假定形状（可能不连续），而在经典 Lapack 中，它们是假定大小（因此是连续的）。

    由于调用链的末尾（经典 Lapack）需要连续的数组，因此如果您想处理不连续的数组，无论如何都必须进行临时复制。唯一的问题是“在哪里”：

    • associate( eig_vects=>mat(:10000,:10000) )不进行任何复制，eigs_vects只是此处的一种别名
    • call syevd(eig_vects,eig_evals,'V')也不会进行任何复制，因为 Lapack95 的接口syevd可以处理不连续的数组。
    • 因此，临时复制很可能发生在对 Lapack 例程的内部调用中*syevd；但您无法控制其执行方式。例如，通常，Intel 编译器使用堆栈来处理临时数组，而对于如此大的矩阵，它将失败。

    正如@VladimirFГероямслава 所建议的，您最好明确分配一个连续的数组来一次性使用：

    real, allocatable :: eigs_vect(:,:)
allocate( eigs_vect, source=mat(:10000,:10000) )

    请注意，它vec(:10000)是连续的，因此您不需要为其分配另一个数组。

    请注意，您也可以直接调用经典的 Lapack 例程，因为它们的接口可以处理这种情况而无需任何复制，这要归功于问题大小和“主要维度”的解耦：

    call ssyevd( 'V', 'U', 10000, mat, 15000, vec & 
           , work, lwork, iwork, liwork, stat )
    • 1

相关问题