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 может быть передана в виде матрицы, функции-обработчика или встроенной функции Afcn таким образом, чтобы Afcn (x, "notransp") = A * x и Afcn (x, "transp") = A' * x. Дополнительные параметры для Afcn могут быть переданы после 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, рассматриваются как параметры и передаются соответствующим образом любым функциям (Afcn или Mfcn) или тем, которые были предоставлены 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;
Afcn = @(x, string) strcmp (string, "notransp") * (A * x) + ...
                     strcmp (string, "transp") * (A' * x);
Mfcn = @(x, string) strcmp (string, "notransp") * (M \ x) + ...
                     strcmp (string, "transp") * (M' \ x);
M1fcn = @(x, string) strcmp (string, "notransp") * (M1 \ x) + ...
                     strcmp (string, "transp") * (M1' \ x);
M2fcn = @(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 (Afcn, b, [], n)

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

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

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

x = bicg (Afcn, b, 1e-6, n, Mfcn)

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

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

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

x = bicg (Afcn, b, 1e-6, n, M1fcn, M2fcn)

ПРИМЕР 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

Apfcn = @(x, string, p) Ap (A, x, string, p);
x = bicg (Apfcn, 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 может быть передана как матрица, функция-обработчик или встроенная функция Afcn таким образом, что Afcn(x) = A * x. Дополнительные параметры для Afcn передаются после 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;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

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

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

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

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

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

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

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

x = bicgstab (Afcn, b, 1e-6, n, Mfcn)

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

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

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

x = bicgstab (Afcn, b, 1e-6, n, M1fcn, M2fcn)

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

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = bicgstab (Apfcn, 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 может быть передана как матрица, функция-обработчик или встроенная функция Afcn таким образом, что Afcn(x) = A * x. Дополнительные параметры для Afcn передаются после 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;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

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

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

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

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

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

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

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

x = cgs (Afcn, b, 1e-6, n, Mfcn)

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

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

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

x = cgs (Afcn, b, 1e-6, n, M1fcn, M2fcn)

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

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = cgs (Apfcn, 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 может быть передана как матрица, дескриптор функции или встроенная функция Afcn таким образом, что Afcn(x) = A * x. Дополнительные параметры для Afcn передаются после 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;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

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

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

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

x = gmres (Afcn, 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 (Afcn, b, [], 1e-6, n, Mfcn)

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

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

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

x = gmres (Afcn, b, 1e-6, n, M1fcn, M2fcn)

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

function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfcn = @(x, p) Ap (A, x, p);
x = gmres (Apfcn, 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 and N. Nachtigal, QMR: a quasi-minimal residual method for non-Hermitian linear systems, Numerische Mathematik, 1991, 60, pp. 315–339.
2. R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhour, R. Pozo, C. Romine, and H. van der Vorst, Templates for the solution of linear systems: Building blocks for iterative methods, SIAM, 2nd ed., 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 может быть передана в виде матрицы, функции-обработчика или встроенной функции Afcn таким образом, что Afcn(x) = A * x. Дополнительные параметры для Afcn передаются после 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;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

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

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

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

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

ПРИМЕР 3: tfqmr с матрицей предобуславливания M

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

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

x = tfqmr (Afcn, b, 1e-6, n, Mfcn)

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

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

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

x = tfqmr (Afcn, b, 1e-6, n, M1fcn, M2fcn)

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

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = tfqmr (Apfcn, 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.