LAPACK例程中的WORK参数有什么用？

Question

我正在计算具有 scipy.linalg.cython_lapack.syev 的对称矩阵的特征值分解。从 doc I found，我需要传递一个名为 WORK:

的数组

WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

但是，我看不到它的作用（无法理解执行后的值是什么），也看不到它的用途。这个参数的作用是什么？

Answer 1

syev 在计算过程中需要额外的 space 并且调用者必须提供此内存（work 数组）。

计算需要最少的额外内存，但如果您有能力分配更多内存，计算将从中获益并且速度更快。

工作数组的最小大小应为 (nb +2)n，其中 nb 是块大小（可以通过 ilaenv 计算）。

通常，一个人会提供更多内存：一个人要么知道最佳大小，因为矩阵的结构是已知的，要么通过调用 syev 和 lwork = -1，这将 return 最佳大小（在 work 数组中）：

double query;
int lwork = -1;
dsyev(..., &query, lwork);//query the work-size (error handling missing) 
lwork=(int) query; //cast the returned double to int to get the size
....

work 的大部分内容只是一些垃圾，我想您可以通过查看数据来深入了解算法的工作 - 但这取决于实现。

然而，work的第一个条目将是工作数组的推荐大小，如果您通过设置lwork = -1查询推荐大小也是如此。

Answer 2

使用 cython 接口从 scipy.linalg.cython_lapack 到 dsyev() 是有意义的：numpy eigh wraps dsyevs and scipy eigh wraps dsyevr()。但是，按照 dsyev() 的 Fortran 原型，必须提供一个数组 WORK。

syev 需要数组 WORK 供内部使用（如果 LWORK = -1 除外）。

LAPACK写在Fortran 77 and this langage does not support dynamic allocation on the heap in its standards！动态分配可能依赖于平台或由特定的编译器扩展提供。因此，编写 LAPACK 以便用户可以使用任何 she/he 想要的东西：静态数组、分配在堆栈上的数组或分配在堆上的数组。

的确，在库中对 WORK 数组的大小进行硬编码会触发两种尴尬情况。要么数组太大，白白增加内存占用，要么数组太小，导致性能不佳或越界错误（分段错误...）。因此，内存管理留给了库的用户。如果 LWORK = -1.

提供了数组的最佳大小，则会向用户提供一些帮助

如果动态分配可用，LAPACK 函数最常见的用法是首先使用 LWORK = -1 执行工作区查询，然后使用 return 值分配 WORK 数组大小正确，最后调用LAPACK的例程得到预期的结果。 LAPACK 的高端包装器（如 LAPACKE）的功能就是这样做的：看看 source of LAPACKE for function LAPACKE_dsyev()！它调用函数 LAPACKE_dsyev_work 两次，后者调用 LAPACK_dsyev（包装 dsyev()）。

包装器仍然具有 LAPACKE_dsyev_work() 等函数，其中参数 work 和 lwork 仍然是必需的。因此，如果通过不在调用之间释放 WORK 来多次调用类似大小的例程，则可以减少分配的数量，但用户必须自己执行此操作（参见此 example). In addition, the source of ILAENV，LAPACK 的函数调用计算 WORK 的优化大小，具有以下文本：

This version provides a set of parameters which should give good, but not optimal, performance on many of the currently available computers. Users are encouraged to modify this subroutine to set the tuning parameters for their particular machine using the option and problem size information in the arguments.

因此，测试大于工作区查询 return 大小的 WORK 可以提高性能。

事实上，LAPACK 中的许多函数都具有 WORK 和 LWORK 参数。如果您通过 grep -r "alloc" . 在文件夹 lapack-3.7.1/SRC 中搜索 alloc，则输出仅包含注释行：

    ./zgejsv.f:*>          Length of CWORK to confirm proper allocation of workspace.
    ./zgejsv.f:*>            In both cases, the allocated CWORK can accommodate blocked runs
    ./zgejsv.f:*>          Length of RWORK to confirm proper allocation of workspace.
    ./zgesdd.f:*       minimal amount of workspace allocated at that point in the code,
    ./zhseqr.f:*     ==== NL allocates some local workspace to help small matrices
    ./dhseqr.f:*     ==== NL allocates some local workspace to help small matrices
    ./dgesdd.f:*       minimal amount of workspace allocated at that point in the code,
    ./shseqr.f:*     ==== NL allocates some local workspace to help small matrices
    ./chseqr.f:*     ==== NL allocates some local workspace to help small matrices
    ./sgesdd.f:*       minimal amount of workspace allocated at that point in the code,
    ./sgejsv.f:*>          Length of WORK to confirm proper allocation of work space.
    ./cgejsv.f:*>          Length of CWORK to confirm proper allocation of workspace.
    ./cgejsv.f:*>            In both cases, the allocated CWORK can accommodate blocked runs
    ./cgejsv.f:*>          Length of RWORK to confirm proper allocation of workspace.
    ./dgejsv.f:*>          Length of WORK to confirm proper allocation of work space.
    ./cgesdd.f:*       minimal amount of workspace allocated at that point in the code,

这表明 LAPACK 的核心不会通过像 allocate 这样的命令处理堆上的动态内存分配，对于大型数组来说 useful：用户必须自己关心。

LAPACK例程中的WORK参数有什么用？

What is the use of the WORK parameters in LAPACK routines?

python

cython

lapack

lapacke