python - 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

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

syev 需要数组 WORK 供内部使用(预计 LWORK = -1)。

LAPACK 写在 Fortran 77 中和这种语言 does not support dynamic allocation on the heap in its standards !动态分配可能依赖于平台或由特定的编译器扩展提供。因此，编写 LAPACK 以便用户可以使用她/他想要的任何东西:静态数组、分配在堆栈上的数组或分配在堆上的数组。

事实上，在库中硬编码 WORK 数组的大小会触发两种尴尬的情况。要么数组太大，白白增加内存占用，要么数组太小，导致性能不佳或越界错误(段错误...)。因此，内存管理留给了库的用户。如果 LWORK = -1，则会为用户提供一些帮助，因为会提供数组的最佳大小。

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

包装器仍然具有 LAPACKE_dsyev_work() 等函数，其中参数 work 和 lwork 仍然是必需的。因此，如果通过不在调用之间释放 WORK 来多次调用类似大小的例程，则可以减少分配的数量，但用户必须自己执行此操作(参见 example )。此外，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.

因此，测试大于工作区查询返回大小的 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对于大型阵列:用户必须自己处理。