I am computing an eigenvalue decomposition of a symmetric matrix with scipy.linalg.cython_lapack.syev. From the doc I found, I need to pass an array called WORK:
WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
However, I can't see what it does (can't understand what the values after execution are), nor what it's used for. What is the purpose of this parameter?
Using the cython interface to dsyev() from scipy.linalg.cython_lapack makes sense: numpy eigh wraps dsyevs and scipy eigh wraps dsyevr(). But, following the Fortran prototype of dsyev(), an array WORK must be provided.
The array WORK is required by syev for internal use (expect if LWORK = -1).
LAPACK is written in Fortran 77 and this langage does not support dynamic allocation on the heap in its standards! Dynamic allocation might have been plateform-dependent or provided by specific compiler extentions. Consequently, LAPACK is written so that the user could use whatever she/he wants: static arrays, arrays allocated on the stack or array allocated on the heap.
Indeed, hard-coding the size for the WORK array in the library would trigger two ackward situations. Either the array is too big, increasing the memory footprint for nothing, or the array is too small, leading to bad performance or out-of-bound errors (segmentation faults...). As a result, the memory management is left to the user of the library. Some help is provided to the user as the optimal size for the array is provided if LWORK = -1.
If dynamic allocation is available, the most common use of LAPACK functions is to first perform a workspace query using LWORK = -1, then use the return value to allocate a WORK array of the correct size and finally call the routine of LAPACK to get the expected result. High-end wrappers of LAPACK such as LAPACKE features function doing just that: take a look at the source of LAPACKE for function LAPACKE_dsyev()! It calls twice the function LAPACKE_dsyev_work, which calls LAPACK_dsyev (wrapping dsyev()).
Wrappers still feature functions such as LAPACKE_dsyev_work(), where the arguments work and lwork are still required. The number of allocations can therefore be reduced if the routine is called multiple times on similar sizes by not deallocating WORK between calls, but the user must do that himself (see this example). In addition, the source of ILAENV, the function of LAPACK called to compute the optimzed size of WORK, features the following text:
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.
As a result, testing sizes of WORK larger than the size returned by the workspace query could improve performances.
Indeed, lots of functions in LAPACK feature the WORK and LWORK arguments. If you search for alloc in folder lapack-3.7.1/SRC by grep -r "alloc" .
, the output only features comment lines:
./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,
It shows that the core of LAPACK does not handle dynamic memory allocation on the heap by commands like allocate, which is useful for large arrays: the user must care for that himself.