Wdrożenie metody Jacobi-Davidsona dla problemu wartości własnej sześciennej

Mam duży problem wartości własnej sześciennej:

(A_{0} + λ A_{1} + λ^{2} A_{2} + λ^{3} A_{3}) x = 0.

$\left(\mathbf{A}_0 + \lambda\mathbf{A}_1 + \lambda^2\mathbf{A}_2 + \lambda^3\mathbf{A}_3\right)\mathbf{x} = 0.$

Mógłbym to rozwiązać, przechodząc na liniowy problem wartości własnych, ale spowodowałoby to, że układ byłby tak duży: $3^2$

[\begin{matrix} - A_{0} & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & I \end{matrix}] [\begin{matrix} x \\ y \\ z \end{matrix}] = λ [\begin{matrix} A_{1} & A_{2} & A_{3} \\ I & 0 & 0 \\ 0 & I & 0 \end{matrix}] [\begin{matrix} x \\ y \\ z \end{matrix}],

$\begin{bmatrix} -\mathbf{A}_0 & 0 & 0 \\ 0 & \mathbf{I} & 0 \\ 0 & 0 & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{x} \\ \mathbf{y} \\ \mathbf{z} \end{bmatrix} = \lambda \begin{bmatrix} \mathbf{A}_1 & \mathbf{A}_2 & \mathbf{A}_3 \\ \mathbf{I} & 0 & 0 \\ 0 & \mathbf{I} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{x} \\ \mathbf{y} \\ \mathbf{z} \end{bmatrix},$

gdzie i . Jakie inne techniki są dostępne, aby rozwiązać problem wartości własnej sześciennej? Słyszałem, że istnieje wersja Jacobi-Davidson, która ją rozwiąże, ale nie znalazła implementacji. $\mathbf{y} = \lambda\mathbf{x}$ $\mathbf{z} = \lambda\mathbf{y}$

Ponadto muszę być w stanie celować w określone wartości własne podobnie do metody przesunięcia i odwracania ARPACK i znaleźć powiązane wektory własne.

c++ eigensystem eigenvalues

— OSE
źródło

Jakie są wymiary zaangażowanych macierzy?

— Bill Barth

A_{i}

$\mathbf{A}_i$ jest zamówieniem

10000 \times 10000

$10000\times 10000$ . Mam dwa różne sformułowania tego problemu, jeden w którym

A_{i}

$\mathbf{A}_i$ jest gęsty, a drugi jest rzadki.

— OSE

SLEPc ma procedury dla problemów z kwadratową wartością własną i nieliniowymi problemami z wartością własną, więc możesz być w stanie znaleźć tam, czego potrzebujesz. Ma również funkcje zmiany i odwracania oraz interfejs do ARPACK.

— Geoff Oxberry

Dzięki protokołowi odwrotnej komunikacji ARPACK nie trzeba jawnie przechowywać macierzy : wystarczy podać dwie funkcje, które obliczają: $3n \times 3n$

$\left[ \begin{array}{c} x \\ y \\ z \end{array}\right] \rightarrow \left[ \begin{array}{c} -A_0 x \\ y \\ z \end{array}\right]$ i $\left[ \begin{array}{c} x \\ y \\ z \end{array}\right] \rightarrow \left[ \begin{array}{c} A_1 x + A_2 y + A_3 z \\ y \\ z \end{array}\right]$

(nadal płacisz cenę za przechowywanie wektorów wymiarowych ale nic nie płacisz za macierze). $3\times n$

Jeśli chodzi o inwertowany przekształcać, można zrobić to samo, czyli wdrożenia go samodzielnie za pomocą wywołania zwrotnego, który oblicza zamiast i zastąpić komputerowej z . Aby obliczyć , możesz wstępnie wyliczyć swoją macierz , co oznacza tylko wstępne faktoring (przy użyciu LU, Cholesky'ego lub ich rzadkich wersji, w zależności od struktury macierzy). W przypadku pełnej transformacji z odwróconym przesunięciem myślę, że można zrobić coś podobnego. $x \mapsto M^{-1} x$ $x \mapsto M x$ $\lambda's$ $\lambda^{-1}$ $M^{-1}x$ $M$ $A_0$

— BrunoLevy
źródło