18.5 Специализированные решатели

: x = bicg (A, b) ¶

: x = bicg (A, b, tol) ¶

: x = bicg (A, b, tol, maxit) ¶

: x = bicg (A, b, tol, maxit, M) ¶

: x = bicg (A, b, tol, maxit, M1, M2) ¶

: x = bicg (A, b, tol, maxit, M, [], x0) ¶

: x = bicg (A, b, tol, maxit, M1, M2, x0) ¶

: x = bicg (A, b, tol, maxit, M, [], x0, …) ¶

: x = bicg (A, b, tol, maxit, M1, M2, x0, …) ¶

: [x, flag, relres, iter, resvec] = bicg (A, b, …) ¶

Решает систему линейных уравнений A * x = b с помощью итерационного метода бисопряжённых градиентов.

Входные аргументы:

• A — матрица системы линейных уравнений, которая должна быть квадратной. A может быть передана как матрица, дескриптор функции или функция inline Afun таким образом, что Afun (x, "notransp") = A * x и Afun (x, "transp") = A' * x. Дополнительные параметры для Afun могут быть переданы после x0.
• b — вектор правой части. Он должен быть столбцовым вектором с тем же количеством строк, что и у A.
• tol — требуемая относительная точность для ошибки остатка, b - A * x. Итерации прекращаются, если norm (b - A * x) ≤ tol * norm (b). Если tol опущено или пусто, используется значение 1e-6.
• maxit — максимальное количество итераций; если maxit опущено или пусто, используется значение 20.
• M1, M2 — прекондиционеры. Прекондиционер M задаётся как M = M1 * M2. M1 и M2 могут быть переданы как матрица или как дескриптор функции или функция inline g таким образом, что g (x, "notransp") = M1 \ x или g (x, "notransp") = M2 \ x и g (x, "transp") = M1' \ x или g (x, "transp") = M2' \ x. Если M1 опущено или пусто, прекондиционирование не применяется. Прекондиционированная система теоретически эквивалентна применению метода bicg к системе линейных уравнений inv (M1) * A * inv (M2) * y = inv (M1) * b и inv (M2') * A' * inv (M1') * z = inv (M2') * b и затем устанавливанию x = inv (M2) * y.
• x0 — начальное приближение. Если x0 опущено или пусто, функция по умолчанию устанавливает x0 в нулевой вектор.

Все аргументы, которые следуют за x0, рассматриваются как параметры и передаются надлежащим образом любым функциям (Afun или Mfun) или тем, что были заданы для bicg.

Выходные параметры:

• x — вычисленное приближение решения A * x = b. Если алгоритм не сошёлся, то x — итерация с минимальным остатком.
• flag указывает на результат выхода:
  • 0: Алгоритм сошёлся в пределах заданной точности.
  • 1: Алгоритм не сошёлся и достиг максимального количества итераций.
  • 2: Матрица прекондиционера вырождена.
  • 3: Алгоритм зациклился, т.е. абсолютное значение разницы между текущей итерацией x и предыдущей меньше, чем eps * norm (x,2).
  • 4: Алгоритм не смог продолжить, потому что промежуточные значения стали слишком маленькими или слишком большими для надёжного вычисления.
• relres — отношение конечного остатка к начальному значению, измеренное в евклидовой норме.
• iter — итерация, на которой вычисляется x.
• resvec — вектор, содержащий остаток на каждой итерации. Общее количество выполненных итераций равно length (resvec) - 1.

Рассмотрим тривиальную задачу с трёхдиагональной матрицей

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
              sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A);  # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afun = @(x, string) strcmp (string, "notransp") * (A * x) + ...
                     strcmp (string, "transp") * (A' * x);
Mfun = @(x, string) strcmp (string, "notransp") * (M \ x) + ...
                     strcmp (string, "transp") * (M' \ x);
M1fun = @(x, string) strcmp (string, "notransp") * (M1 \ x) + ...
                     strcmp (string, "transp") * (M1' \ x);
M2fun = @(x, string) strcmp (string, "notransp") * (M2 \ x) + ...
                     strcmp (string, "transp") * (M2' \ x);

ПРИМЕР 1: простое использование bicg

x = bicg (A, b)

ПРИМЕР 2: bicg с функцией, которая вычисляет A*x и A'*x

x = bicg (Afun, b, [], n)

ПРИМЕР 3: bicg с матрицей прекондиционера M

x = bicg (A, b, 1e-6, n, M)

ПРИМЕР 4: bicg с функцией в качестве прекондиционера

x = bicg (Afun, b, 1e-6, n, Mfun)

ПРИМЕР 5: bicg с матрицами прекондиционеров M1 и M2

x = bicg (A, b, 1e-6, n, M1, M2)

ПРИМЕР 6: bicg с функциями в качестве прекондиционеров

x = bicg (Afun, b, 1e-6, n, M1fun, M2fun)

ПРИМЕР 7: bicg с функцией, требующей аргумента

function y = Ap (A, x, string, z)
  ## compute A^z * x or (A^z)' * x
  y = x;
  if (strcmp (string, "notransp"))
    for i = 1:z
      y = A * y;
    endfor
  elseif (strcmp (string, "transp"))
    for i = 1:z
      y = A' * y;
    endfor
  endif
endfunction

Apfun = @(x, string, p) Ap (A, x, string, p);
x = bicg (Apfun, b, [], [], [], [], [], 2);

Ссылка:

Y. Saad, Итеративные методы для разреженных систем линейных уравнений, второе издание, 2003, SIAM.

См. также: bicgstab, cgs, gmres, pcg, qmr, tfqmr.

: x = bicgstab (A, b, tol, maxit, M1, M2, x0, …) ¶

: x = bicgstab (A, b, tol, maxit, M, [], x0, …) ¶

: [x, flag, relres, iter, resvec] = bicgstab (A, b, …) ¶

Решить A x = b с помощью итерационного метода стабилизированного бисопряжённого градиента.

Входные параметры:

• - A — матрица линейной системы, она должна быть квадратной. A может быть передана как матрица, дескриптор функции или встроенная функция Afun таким образом, что Afun(x) = A * x. Дополнительные параметры для Afun передаются после x0.
• - b — вектор правой части. Он должен быть столбцовым вектором с таким же числом строк, как у A.
• - tol — требуемая относительная точность для остаточной ошибки, b - A * x. Итерация останавливается, если norm (b - A * x) ≤ tol * norm (b). Если tol опущено или пусто, используется толерантность 1e-6.
• - maxit — максимальное число внешних итераций, если не задано или установлено в [], используется значение по умолчанию min (20, numel (b)).
• - M1, M2 — прекондиционеры. Прекондиционер M задаётся как M = M1 * M2. Оба M1 и M2 могут быть переданы как матрица или как дескриптор функции или встроенная функция g таким образом, что g(x) = M1 \ x или g(x) = M2 \ x. Используется правое прекондиционирование, т.е., сначала решается A * inv (M) * y = b и затем x = inv (M) * y.
• - x0 — начальное приближение, если не задано или установлено в [], используется значение по умолчанию zeros (size (b)).

Аргументы, следующие за x0, обрабатываются как параметры и передаются должным образом любым функциям (A или M), которые передаются в bicstab.

Выходные параметры:

• - x — вычисленное приближение. Если метод не сходится, то это итерация с минимальным остатком.
• - flag указывает статус завершения:
  • - 0: итерация сошлась в пределах выбранной толерантности
  • - 1: максимальное число итераций достигнуто до сходимости
  • - 2: матрица прекондиционера является вырожденной
  • - 3: алгоритм достиг застоя
  • - 4: алгоритм не может продолжить из-за деления на ноль
• - relres — относительный остаток, полученный как (A*x-b) / norm(b).
• - iter — (возможно половина) итерация, на которой вычисляется x. Если это половина итерации, то это iter + 0.5
• - resvec — вектор, содержащий остаток каждой полу- и полной итерации (Также есть полу-итерации, так как x вычисляется в два этапа на каждой итерации). Используя (length(resvec) - 1) / 2 можно увидеть общее количество (полных) итераций, выполненных.

Рассмотрим тривиальную задачу с треугольной матрицей

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

ПРИМЕР 1: простейшее использование bicgstab

x = bicgstab (A, b, [], n)

ПРИМЕР 2: bicgstab с функцией, которая вычисляет A * x

x = bicgstab (Afun, b, [], n)

ПРИМЕР 3: bicgstab с матрицей-прекондиционером M

x = bicgstab (A, b, [], 1e-06, n, M)

ПРИМЕР 4: bicgstab с функцией в качестве прекондиционера

x = bicgstab (Afun, b, 1e-6, n, Mfun)

ПРИМЕР 5: bicgstab с матрицами-прекондиционерами M1 и M2

x = bicgstab (A, b, [], 1e-6, n, M1, M2)

ПРИМЕР 6: bicgstab с функциями в качестве прекондиционеров

x = bicgstab (Afun, b, 1e-6, n, M1fun, M2fun)

ПРИМЕР 7: bicgstab со входом, требующим аргумент

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfun = @(x, string, p) Ap (A, x, string, p);
x = bicgstab (Apfun, b, [], [], [], [], [], 2);

ПРИМЕР 8: явный пример, чтобы показать, что bicgstab использует правое прекондиционирование

[M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by bicgstab after one iteration
[x_ref, fl] = bicgstab (A, b, [], 1, M)

## right preconditioning
[y, fl] = bicgstab (A / M, b, [], 1)
x = M \ y # compare x and x_ref

Ссылки:

Y. Saad, Итерационные методы для разреженных линейных систем, второе издание, 2003, SIAM

См. также: bicg, cgs, gmres, pcg, qmr, tfqmr.

: x = cgs (A, b, tol, maxit, M1, M2, x0, …) ¶

: x = cgs (A, b, tol, maxit, M, [], x0, …) ¶

: [x, flag, relres, iter, resvec] = cgs (A, b, …) ¶

Решить A x = b, где A — квадратная матрица, используя метод сопряжённых градиентов в квадрате.

Входные аргументы:

• - A — матрица линейной системы, она должна быть квадратной. A может быть передана как матрица, дескриптор функции или встроенная функция Afun таким образом, что Afun(x) = A * x. Дополнительные параметры для Afun передаются после x0.
• - b — вектор правой части. Он должен быть столбцовым вектором с тем же числом строк, что и у A.
• - tol — относительная точность, если не задано или установлено в [], используется значение по умолчанию 1e-6.
• - maxit — максимальное число внешних итераций, если не задано или установлено в [], используется значение по умолчанию min (20, numel (b)).
• - M1, M2 — прекондиционеры. Матрица прекондиционера задаётся как M = M1 * M2. Оба M1 и M2 могут быть переданы как матрица или как дескриптор функции или встроенная функция g таким образом, что g(x) = M1 \ x или g(x) = M2 \ x. Если M1 пусто или не передано, прекондиционеры не применяются. Используется правое прекондиционирование, т.е., сначала решается A*inv(M)*y = b и затем x = inv(M)*y.
• - x0 — начальное приближение, если не задано или установлено в [], используется значение по умолчанию zeros (size (b)).

Аргументы, следующие за x0, обрабатываются как параметры и передаются должным образом любым функциям (A или P), которые передаются в cgs.

Выходные параметры:

• - x — вычисленное приближение. Если метод не сходится, то это итерация с минимальным остатком.
• - flag указывает статус завершения:
  • - 0: итерация сошлась в пределах выбранной толерантности
  • - 1: максимальное число итераций достигнуто до сходимости
  • - 2: матрица прекондиционера является вырожденной
  • - 3: алгоритм достиг застоя
  • - 4: алгоритм не может продолжить из-за деления на ноль
• - relres — относительный остаток, полученный как (A*x-b) / norm(b).
• - iter — итерация, на которой вычисляется x.
• - resvec — вектор, содержащий остаток на каждой итерации. Используя length(resvec) - 1 можно увидеть общее количество выполненных итераций.

Рассмотрим тривиальную задачу с треугольной матрицей

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)'
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

ПРИМЕР 1: простейшее использование cgs

x = cgs (A, b, [], n)

ПРИМЕР 2: cgs с функцией, которая вычисляет A * x

x = cgs (Afun, b, [], n)

ПРИМЕР 3: cgs с матрицей-прекондиционером M

x = cgs (A, b, [], 1e-06, n, M)

ПРИМЕР 4: cgs с функцией в качестве прекондиционера

x = cgs (Afun, b, 1e-6, n, Mfun)

ПРИМЕР 5: cgs с матрицами-прекондиционерами M1 и M2

x = cgs (A, b, [], 1e-6, n, M1, M2)

ПРИМЕР 6: cgs с функциями в качестве прекондиционеров

x = cgs (Afun, b, 1e-6, n, M1fun, M2fun)

ПРИМЕР 7: cgs со входом, требующим аргумент

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfun = @(x, string, p) Ap (A, x, string, p);
x = cgs (Apfun, b, [], [], [], [], [], 2);

ПРИМЕР 8: явный пример, чтобы показать, что cgs использует правое прекондиционирование

[M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by cgs after one iteration
[x_ref, fl] = cgs (A, b, [], 1, M)

## right preconditioning
[y, fl] = cgs (A / M, b, [], 1)
x = M \ y # compare x and x_ref

Ссылки:

Y. Saad, Итерационные методы для разреженных линейных систем, второе издание, 2003, SIAM

См. также: pcg, bicgstab, bicg, gmres, qmr, tfqmr.

: x = gmres (A, b, restart, tol, maxit, M1, M2, x0, …) ¶

: x = gmres (A, b, restart, tol, maxit, M, [], x0, …) ¶

: [x, flag, relres, iter, resvec] = gmres (A, b, …) ¶

Решить A x = b с помощью итерационного метода GMRES с предварительной обработкой и перестартом, также известного как PGMRES(restart).

Входные аргументы:

• - A — матрица линейной системы, она должна быть квадратной. A может быть передана как матрица, дескриптор функции или встроенная функция Afun таким образом, что Afun(x) = A * x. Дополнительные параметры для Afun передаются после x0.
• - b — вектор правой части. Он должен быть вектором-столбцом с тем же числом строк, что и у A.
• - restart — число итераций перед перестартом метода. Если это [] или N = numel (b), то перестрат не применяется.
• - tol — требуемая относительная погрешность преобразованного остатка, inv (M) * (b - a * x). Итерации прекращаются, если norm (inv (M) * (b - a * x)) ≤ tol * norm (inv (M) * B). Если tol опущена или пуста, используется значение 1e-6.
• - maxit — максимальное число внешних итераций. Если не указано или задано как [], используется значение по умолчанию min (10, N / restart). Обратите внимание, что если restart пусто, то maxit является максимальным числом итераций. Если restart и maxit не пусты, то максимальное число итераций равно restart * maxit. Если оба restart и maxit пусты, максимальное число итераций устанавливается в min (10, N).
• - M1, M2 — предварительные условия. Предварительное условие M задаётся как M = M1 * M2. Оба M1 и M2 могут быть переданы как матрица, дескриптор функции или встроенная функция g таким образом, что g(x) = M1 \ x или g(x) = M2 \ x. Если M1 — [] или не задано, предварительное условие не применяется. Используется левая предварительная обработка, то есть решается inv(M) * A * x = inv(M) * b вместо A * x = b.
• - x0 — начальное приближение. Если не задано или равно [], используется значение по умолчанию zeros (size (b)).

Аргументы, следующие за x0, рассматриваются как параметры и передаются соответствующим образом любой из функций (A или M или M1 или M2), которые переданы gmres.

Выходные данные:

• - x — вычисленное приближение. Если метод не сходится, это итерация с минимальным остатком.
• - flag — код выхода:

  0 : итерация сошлась в пределах заданной погрешности

  1 : максимальное число итераций превышено

  2 : матрица предварительного условия вырождена

  3 : алгоритм достиг стагнации (относительная разница между двумя

  последовательными итерациями меньше eps)

• - relres — значение относительного преобразованного остатка приближения x.
• - iter — вектор, содержащий количество внешних и внутренних итераций, выполненных для вычисления x. Это:
  • iter(1): число внешних итераций, т.е. сколько раз метод был перезапущен. (Если restart пусто или N, то 1, иначе 1 ≤ iter(1) ≤ maxit).
  • iter(2): число итераций, выполненных перед перестартом, т.е. метод перезапускается, когда iter(2) = restart. Если restart пусто или N, то 1 ≤ iter(2) ≤ maxit.

  Для большей ясности, приближение x вычисляется на итерации (iter(1) - 1) * restart + iter(2). Поскольку выходной результат x соответствует решению с минимальным преобразованным остатком, общее число итераций, выполненных методом, задаётся length (resvec) - 1.

• - resvec — вектор, содержащий преобразованные относительные остатки на каждой итерации, включая 0-ю итерацию norm (A * x0 - b).

Рассмотрим тривиальную задачу с треугольной матрицей

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

ПРИМЕР 1: самое простое использование gmres

x = gmres (A, b, [], [], n)

ПРИМЕР 2: gmres с функцией, вычисляющей A * x

x = gmres (Afun, b, [], [], n)

ПРИМЕР 3: использование gmres с перестартом

x = gmres (A, b, restart);

ПРИМЕР 4: gmres с матрицей предварительного условия M с перестартом и без него

x = gmres (A, b, [], 1e-06, n, M)
x = gmres (A, b, restart, 1e-06, n, M)

ПРИМЕР 5: gmres с функцией в качестве предварительного условия

x = gmres (Afun, b, [], 1e-6, n, Mfun)

ПРИМЕР 6: gmres с матрицами предварительных условий M1 и M2

x = gmres (A, b, [], 1e-6, n, M1, M2)

ПРИМЕР 7: gmres с функциями в качестве предварительных условий

x = gmres (Afun, b, 1e-6, n, M1fun, M2fun)

ПРИМЕР 8: gmres в качестве входной функции, требующей аргумент

function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfun = @(x, p) Ap (A, x, p);
x = gmres (Apfun, b, [], [], [], [], [], [], 2);

ПРИМЕР 9: явный пример, демонстрирующий, что gmres использует левую предварительную обработку

[M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by gmres after two iterations
[x_ref, fl] = gmres (A, b, [], [], 1, M)

## left preconditioning
[x, fl] = gmres (M \ A, M \ b, [], [], 1)
x # compare x and x_ref

Ссылка:

Y. Saad, Итерационные методы для разреженных линейных систем, второе издание, 2003, SIAM

См. также: bicg, bicgstab, cgs, pcg, pcr, qmr, tfqmr.

: x = qmr (A, b, rtol, maxit, M1, M2, x0) ¶

: x = qmr (A, b, rtol, maxit, P) ¶

: [x, flag, relres, iter, resvec] = qmr (A, b, …) ¶

Решить A x = b с помощью итерационного метода Quasi-Minimal Residual (без прогнозирования).

• - rtol — относительная погрешность, если не указано или равно [], используется значение по умолчанию 1e-6.
• - maxit — максимальное число внешних итераций, если не указано или равно [], используется значение по умолчанию min (20, numel (b)).
• - x0 — начальное приближение, если не указано или равно [], используется значение по умолчанию zeros (size (b)).

A может быть передана как матрица или как дескриптор функции или встроенная функция f таким образом, что f(x, "notransp") = A*x и f(x, "transp") = A'*x.

Предварительное условие P задаётся как P = M1 * M2. Оба M1 и M2 могут быть переданы как матрица, дескриптор функции или встроенная функция g таким образом, что g(x, "notransp") = M1 \ x или g(x, "notransp") = M2 \ x и g(x, "transp") = M1' \ x или g(x, "transp") = M2' \ x.

Если вызов осуществляется с более чем одним выходным параметром

• - flag — код выхода:
  • - 0: итерация сошлась в пределах выбранной точности
  • - 1: максимальное число итераций было достигнуто до сходимости
  • - 3: алгоритм достиг стагнации

  (значение 2 не используется, но пропущено для совместимости).

• - relres — конечное значение относительного остатка.
• - iter — число выполненных итераций.
• - resvec — вектор, содержащий нормы остатка на каждой итерации.

Ссылки:

1. R. Freund и N. Nachtigal, QMR: метод квазиминимального остатка для неэрмитовых линейных систем, Numerische Mathematik, 1991, 60, стр. 315–339.
2. R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhour, R. Pozo, C. Romine, и H. van der Vorst, Шаблоны для решения линейных систем: строительные блоки для итерационных методов, SIAM, 2-е изд., 1994.

См. также: bicg, bicgstab, cgs, gmres, pcg.

: x = tfqmr (A, b, tol, maxit, M1, M2, x0, …) ¶

: x = tfqmr (A, b, tol, maxit, M, [], x0, …) ¶

: [x, flag, relres, iter, resvec] = tfqmr (A, b, …) ¶

Решить A x = b с помощью метода qmr с транспонированным деревом, основанного на методе cgs.

Параметры входных данных:

• - A — матрица линейной системы, она должна быть квадратной. A может быть передана в виде матрицы, функции-обработчика или встроенной функции Afun таким образом, чтобы Afun(x) = A * x. Дополнительные параметры для Afun передаются после x0.
• - b — вектор правой части. Он должен быть вектором-столбцом с тем же числом строк, что и у A.
• - tol — относительная погрешность; если не задана или установлена в [], используется значение по умолчанию 1e-6.
• - maxit — максимальное число внешних итераций; если не задана или установлена в [], используется значение по умолчанию min (20, numel (b)). Для совместимости, поскольку метод имеет разное поведение при нечетном или четном числе итераций, рассматривается как итерация в tfqmr весь цикл нечетных-четных итераций. То есть, для выполнения одной полной итерации алгоритм выполняет две под-итерации: нечетную и четную.
• - M1, M2 — прекондиционеры. Прекондиционер M задается как M = M1 * M2. M1 и M2 могут быть переданы в виде матрицы, функции-обработчика или встроенной функции g таким образом, чтобы g(x) = M1 \ x или g(x) = M2 \ x. Используется метод правого прекондиционирования, т.е. сначала решается A*inv(M)*y = b и затем x = inv(M)*y, а не A x = b.
• - x0 — начальное приближение; если не задано или установлено в [], используется значение по умолчанию zeros (size (b)).

Аргументы, следующие за x0, обрабатываются как параметры и передаются надлежащим образом в любые функции (A или M), которые передаются в tfqmr.

Параметры вывода:

• - x — вычисленное приближение. Если метод не сходится, возвращается значение с минимальным остатком.
• - flag — показывает статус завершения:
  • - 0: итерации сошлись в пределах выбранной точности
  • - 1: максимальное число итераций достигнуто до сходимости
  • - 2: матрица прекондиционера вырождена
  • - 3: алгоритм достиг стагнации
  • - 4: алгоритм не может продолжить из-за деления на ноль
• - relres — относительный остаток, вычисленный как (A*x-b) / norm (b).
• - iter — итерация, на которой вычислено x.
• - resvec — вектор, содержащий остаток на каждой итерации (включая norm (b - A x0)). Выполнив length (resvec) - 1, можно увидеть общее количество выполненных итераций.

Рассмотрим тривиальную задачу с треугольной матрицей

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)'
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

ПРИМЕР 1: самое простое использование tfqmr

x = tfqmr (A, b, [], n)

ПРИМЕР 2: tfqmr с функцией, которая вычисляет A * x

x = tfqmr (Afun, b, [], n)

ПРИМЕР 3: tfqmr с матрицей прекондиционера M

x = tfqmr (A, b, [], 1e-06, n, M)

ПРИМЕР 4: tfqmr с функцией в качестве прекондиционера

x = tfqmr (Afun, b, 1e-6, n, Mfun)

ПРИМЕР 5: tfqmr с матрицами прекондиционеров M1 и M2

x = tfqmr (A, b, [], 1e-6, n, M1, M2)

ПРИМЕР 6: tfmqr с функциями в качестве прекондиционеров

x = tfqmr (Afun, b, 1e-6, n, M1fun, M2fun)

ПРИМЕР 7: tfqmr со входной функцией, требующей аргумента

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfun = @(x, string, p) Ap (A, x, string, p);
x = tfqmr (Apfun, b, [], [], [], [], [], 2);

ПРИМЕР 8: явный пример, показывающий, что tfqmr использует правый прекондиционер

[M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by tfqmr after one iteration
[x_ref, fl] = tfqmr (A, b, [], 1, M)

## right preconditioning
[y, fl] = tfqmr (A / M, b, [], 1)
x = M \ y # compare x and x_ref

Ссылка:

Y. Saad, Итерационные методы для разреженных линейных систем, второе издание, 2003, SIAM

См. также: bicg, bicgstab, cgs, gmres, pcg, qmr, pcr.