Итерационные решатели как итераторы

В расширенных случаях вам может потребоваться получить доступ к внутренним структурам данных решателя, внедрять код, который будет выполняться после каждой итерации, иметь полный контроль над выделениями или уменьшить затраты на инициализацию. Все это можно сделать в рамках использования итераторов в IterativeSolvers.jl.

На этом этапе методы BiCGStab(l), CG, Чебышева, GMRES, MINRES, QMR и стационарные методы реализуются как итераторы. Однако пакет пока не экспортирует сами итераторы и вспомогательные методы.

Принципы реализации итераторов

Перечисленные выше решатели, по сути, являются тонкой оболочкой вокруг итератора. Помимо прочего, они инициализируют итерируемый объект, выполняют цикл по итератору и возвращают результат:

function my_solver!(x, A, b)
    iterable = MySolverIterable(x, A, b)
    for item in iterable end
    return iterable.x
end

Вместо вызова my_solver!(x, A, b) можно также самостоятельно инициализировать итерируемый объект и выполнить цикл for.

Пример: отказ от ненужной инициализации

Метод Якоби для SparseMatrixCSC имеет некоторые затраты при инициализации. Нам нужно не только выделить временный вектор, но и найти индексы диагонали (и проверить, что их значения ненулевые). Текущая реализация инициализирует итерируемый объект следующим образом.

jacobi_iterable(x, A::SparseMatrixCSC, b; maxiter::Int = 10) =
    JacobiIterable{eltype(x), typeof(x)}(OffDiagonal(A, DiagonalIndices(A)), x, similar(x), b, maxiter)

Теперь, если вы применяете итерацию Якоби несколько раз с одной и той же матрицей в течение нескольких итераций, имеет смысл инициализировать итерируемую матрицу только один раз и впоследствии использовать ее повторно:

julia> using LinearAlgebra, SparseArrays, IterativeSolvers

julia> A = spdiagm(-1 => -ones(3), 0 => 2*ones(4), 1 => -ones(3));

julia> b1 = [1.0, 2, 3, 4];

julia> b2 = [-1.0, 1, -1, 1];

julia> x = [0.0, -1, 1, 0];

julia> my_iterable = IterativeSolvers.jacobi_iterable(x, A, b1, maxiter = 2);

julia> norm(b1 - A * x) / norm(b1)
1.2909944487358056

julia> for item in my_iterable
           println("Iteration for rhs 1")
       end
Iteration for rhs 1
Iteration for rhs 1

julia> norm(b1 - A * x) / norm(b1)
0.8228507357554791

julia> # Копирование следующей правой части в итерируемый объект
       copyto!(my_iterable.b, b2);

julia> norm(b2 - A * x) / norm(b2)
2.6368778887161235

julia> for item in my_iterable
           println("Iteration for rhs 2")
       end
Iteration for rhs 2
Iteration for rhs 2

julia> norm(b2 - A * x) / norm(b2)
1.610815496107484

Другие варианты использования

В число других вариантов использования входят следующие:

вычисление (гармонических) значений Ритца из матрицы Хессенберга в GMRES;
сравнение приближенной невязки таких методов, как GMRES и BiCGStab(l), с истинной невязкой во время итераций;
обновление предобуславливателя в гибких методах.