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 может быть передана как матрица, дескриптор функции или анонимная функция 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 могут быть переданы как матрица или как дескриптор функции или анонимная функция 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 — квадратная матрица, с помощью метода сопряжённых градиентов (Conjugate Gradients Squared).

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

• - 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, Итеративные методы для разреженных линейных систем, 2-е издание, 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, Итеративные методы для разреженных линейных систем, 2-е издание, 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 с использованием итерационного метода квазиминимального остатка (без предвычисления).

• - 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, Итерационные методы для разреженных линейных систем, 2-е издание, 2003, SIAM

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