Вызов подпрограммы обратной связи матрицы-вектора ARPACK

Question

Я пытаюсь написать драйвер на C++ для вычисления собственных значений для асимметричной действительной разреженной матрицы, используя функции Фортрана, предлагаемые ARPACK, но у меня возникли некоторые проблемы с подходом обратной связи.

Как правило, я пытаюсь решить нормальное уравнение собственных значений:

A*v = lambda*v

и любое взаимодействие с матрицей A осуществляется в ARPACK через функцию 'av':

av(n, workd[ipntr[0]], workd[ipntr[1]])

который умножает вектор, содержащийся в массиве 'workd', начиная с местоположения 'ipntr[0]', и вставляет результат в массив 'workd', начиная с местоположения 'ipntr[1]'. Примеры такого подхода приведены в руководстве по адресу http://www.caam.rice.edu/software/ARPACK/ а также в коде ARPACK/EXAMPLES/SIMPLE/dnsimp.f.

Что я хотел бы знать, так это как на самом деле задействовать матрицу А? Если он не передан в рутину, то как можно найти его действие по предоставленному вектору?

В примере кода dnsimp.f их матрица A вычисляется в функции "av" и "выводится из стандартной дискретизации центральной разности двумерного оператора конвекции-диффузии". Тем не менее, я считаю, что это конкретная проблема? Также не кажется слишком полезным кодировать вывод матрицы A в функцию. Я не могу найти много информации по этому вопросу из руководства.

Кажется, это не слишком большая проблема, поскольку, поскольку это определенная пользователем функция, я могу просто изменить определение "av", чтобы включить матрицу A в качестве параметра. Однако я хотел бы знать, как это делается правильно в случае каких-либо потенциальных проблем с совместимостью.

Спасибо!

Answer 1

Вам не нужно поставлять матрицу в ARPACK.

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

Для получения информации об алгоритмах вы должны взглянуть на руководство пользователя и особенно на главу об основных алгоритмах.

Ответ на комментарий: Базовый алгоритм является формой итерации Арнольди. Основной алгоритм показан в википедии и показывает, что матрица A не будет доступна. Ни прямо, ни косвенно.

В частности, алгоритмы начинаются с произвольного нормализованного вектора q_1. Этот вектор возвращается пользователю. Пользователь умножает этот вектор на матрицу A, используя свою любимую подпрограмму (обычно это эффективное рассеянное умножение матрицы на вектор), и возвращает результат в итерацию Арнольди для вычисления части матрицы Гессенберга H (чьи собственные значения обычно сходятся к крайним собственным значениям). а) и следующий вектор q_2. Это должно быть повторено, пока ваши результаты не сойдутся.