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);

Ссылки:

1. 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

Ссылки:

1. 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

Ссылки:

1. 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

Ссылки:

1. 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 с помощью итерационного метода квази-минимального остатка (без предвычислений).

• - 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 с использованием метода Transpose-Tree 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

Ссылки:

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

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