4

我正在计算对称矩阵的特征值分解scipy.linalg.cython_lapack.syev。从我找到的文档中，我需要传递一个名为 WORK 的数组：

WORK 是 DOUBLE PRECISION 数组，维度 (MAX(1,LWORK)) 退出时，如果 INFO = 0，WORK(1) 返回最优 LWORK。

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

4

2 回答 2

6

从 scipy.linalg.cython_lapack使用 cython 接口dsyev()是有意义的：numpyeigh包装 dsyevs和scipyeigh包装 dsyevr()。但是，按照 Fortran 的原型dsyev()，必须提供一个数组 WORK。

该数组WORK是内部使用所必需的syev（如果是LWORK = -1）。

LAPACK 是用Fortran 77编写的，这种语言在其标准中不支持堆上的动态分配！动态分配可能依赖于平台或由特定的编译器扩展提供。因此，编写 LAPACK 以便用户可以使用她/他想要的任何东西：静态数组、分配在堆栈上的数组或分配在堆上的数组。

实际上，在库中硬编码WORK数组的大小会触发两种不妥的情况。要么数组太大，无谓地增加内存占用，要么数组太小，导致性能不佳或越界错误（分段错误......）。结果，内存管理留给了库的用户。为用户提供了一些帮助，因为如果LWORK = -1.

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

包装器仍然具有功能，例如LAPACKE_dsyev_work(), 其中参数work和lwork仍然是必需的。因此，如果例程在类似大小的情况下通过不在调用之间解除分配而多次WORK调用，则可以减少分配的数量，但用户必须自己这样做（参见此示例）。此外， 的来源ILAENV，调用 LAPACK 的函数来计算 的优化大小WORK，具有以下文本：

此版本提供了一组参数，这些参数应该在许多当前可用的计算机上提供良好但不是最佳的性能。鼓励用户修改此子例程，以使用参数中的选项和问题大小信息为他们的特定机器设置调整参数。

因此，测试 WORK 的大小大于工作空间查询返回的大小可以提高性能。

事实上，LAPACK 中的许多函数都具有WORK和LWORK参数。如果您alloc在文件夹 lapack-3.7.1/SRC by 中搜索grep -r "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，这对于大型数组很有用：用户必须自己关心。

于 2017-10-12T00:25:57.790 回答
2

syev在计算过程中需要额外的空间，调用者必须提供这个内存（work数组）。

计算所需的额外内存量最少，但是如果您有能力分配更多内存，则计算将从中受益并更快。

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

通常，一个人会提供更多内存：一个人要么知道最佳大小，因为矩阵的结构是已知的，要么通过调用syevwith lwork = -1，这将返回最佳大小（在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
....

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

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

于 2017-10-07T11:40:09.597 回答