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 и 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 может быть передана в виде матрицы, функции-обработчика или встроенной функции 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.