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

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