Документация Engee

Сопряженные градиенты (CG)

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

Использование

# IterativeSolvers.cg — Function

cg(A, b; kwargs...) -> x, [history]

То же, что и cg!, но выделяет вектор решений x, инициализированный нулями.

# IterativeSolvers.cg! — Function

cg!(x, A, b; kwargs...) -> x, [history]

Аргументы

  • x: начальное предположение, будет обновляться на месте;

  • A: линейный оператор;

  • b: правая часть.

Ключевые слова

  • statevars::CGStateVariables: имеет 3 массива, аналогичных x, для хранения промежуточных результатов;

  • initially_zero::Bool: если true предполагает iszero(x) так, чтобы одно произведение матрицы на вектор могло быть сохранено при вычислении начального вектора невязки;

  • Pl = Identity(): левый предобуславливатель метода. Должен быть симметричным, положительно определенным как A;

  • abstol::Real = zero(real(eltype(b))), reltol::Real = sqrt(eps(real(eltype(b)))): абсолютный и относительный допуски для условия остановки |r_k| ≤ max(reltol * |r_0|, abstol), где r_k ≈ A * x_k - b является приближенной невязкой в k-й итерации.

note Истинная норма невязки никогда не вычисляется в явном виде во время итераций по соображениям производительности — она может накапливать ошибки округления.

  • maxiter::Int = size(A,2): максимальное количество итераций;

  • verbose::Bool = false: вывод сведений о методе;

  • log::Bool = false: отслеживание нормы невязки в каждой итерации.

Вывод

если log имеет значение false

  • x: приближенное решение.

если log имеет значение true

  • x: приближенное решение.

  • ch: история сходимости.

Ключи ConvergenceHistory

  • :tol => ::Real: допуск остановки.

  • :resnom => ::Vector: норма невязки в каждой итерации.

На GPU

Этот метод должен отлично работать на GPU. В качестве минимального рабочего примера рассмотрим следующее.

using LinearAlgebra, CuArrays, IterativeSolvers

n = 100
A = cu(rand(n, n))
A = A + A' + 2*n*I
b = cu(rand(n))
x = cg(A, b)

Убедитесь, что все векторы состояния хранятся в GPU. Например, при вызове cg!(x, A, b) может возникнуть проблема, когда x хранится в GPU, а b хранится в ЦП — IterativeSolvers.jl не копирует векторы на одно и то же устройство.

Сведения о реализации

Текущая реализация следует довольно стандартному подходу. Отметим, что предобусловленный CG (или PCG) несколько отличается от обычного CG, поскольку в первом случае невязка должна вычисляться явно, а во втором она доступна как побочный продукт. Наша реализация CG гарантирует наличие минимального количества векторных операций.

CG можно использовать как итератор.