Как использовать параметры WORK в подпрограммах LAPACK?

Question

Я вычисляю разложение по собственным значениям симметричной матрицы с scipy.linalg.cython_lapack.syev, Из документа, который я нашел, мне нужно передать массив с именем WORK:

WORK - массив DOUBLE PRECISION, размерность (MAX(1, ​​LWORK)). При выходе, если INFO = 0, WORK(1) возвращает оптимальное значение LWORK.

Тем не менее, я не могу понять, что он делает (не могу понять, каковы значения после выполнения), и для чего он используется. Какова цель этого параметра?

Answer 1

Использование интерфейса Cython для dsyev() из scipy.linalg.cython_lapack имеет смысл: numpy eigh заворачивает дсяев и скипы eigh обертывания дсевр (). Но, следуя прототипу Фортрана 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, что полезно для больших массивов: пользователь должен заботиться об этом сам.

Answer 2

syev во время расчета требуется дополнительное место, и вызывающая сторона должна предоставить эту память (work массив).

Для расчета необходим минимальный объем дополнительной памяти, однако, если вы можете выделить больше памяти, вычисление выиграет от этого и будет быстрее.

Минимальный размер рабочего массива должен быть (nb +2)n где nb это размер блока (который может быть рассчитан с помощью ilaenv).

Обычно каждый дает больше памяти: либо знает оптимальный размер, потому что структура матрицы известна, либо вызывает syev с 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
....

Большая частьwork-content это просто мусор, я думаю, вы могли бы получить представление о работе алгоритма, посмотрев на данные - но это зависит от реализации.

Тем не менее, первая запись work будет рекомендованным размером рабочего массива, это также имеет место, если вы запрашиваете рекомендуемый размер, установив lwork = -1,