Я пытаюсь написать драйвер на 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 в качестве параметра. Однако я хотел бы знать, как это делается правильно в случае каких-либо потенциальных проблем с совместимостью.
Спасибо!
Вам не нужно поставлять матрицу в ARPACK.
Все, что вам нужно сделать, это умножить матрицу на возвращенные векторы (таким образом, обратную связь), пока не будет достигнута желаемая сходимость.
Для получения информации об алгоритмах, вы должны взглянуть на руководство пользователя и особенно на глава о базовых алгоритмах.
Ответ на комментарий: Основной алгоритм является формой Арнольди Итерация. Основной алгоритм показан в википедия и показывает, что матрица А не будет доступна. Ни прямо, ни косвенно.
В частности, алгоритмы начинаются с произвольного нормализованного вектора q_1. Этот вектор возвращается пользователю. Пользователь умножает этот вектор на матрицу A, используя свою любимую подпрограмму (обычно это эффективное рассеянное умножение матрицы на вектор), и возвращает результат в Арнольди Итерация вычислить часть матрицы Гессенберга H (чьи собственные значения обычно сходятся к крайним собственным значениям A) и следующий вектор q_2. Это должно повторяться до тех пор, пока ваши результаты не сойдутся.