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 * xAfcn (x, "transp") = A' * xAfcn 的附加参数可在 x0 之后传递。
  • b 是右侧向量。它必须是列向量,且行数与 A 相同。
  • tol 是残差误差所需的相对容差,即 b - A * x。当 norm (b - A * x)tol * norm (b) 时迭代停止。如果 tol 被省略或为空,则使用 1e-6 的容差。
  • maxit 是允许的最大迭代次数;如果 maxit 被省略或为空,则使用值 20。
  • M1M2 是预条件子。预条件子 M 定义为 M = M1 * M2M1M2 都可以作为矩阵、函数句柄或内联函数 g 传递,使得 g (x, "notransp") = M1 \ xg (x, "notransp") = M2 \ x,且 g (x, "transp") = M1' \ xg (x, "transp") = M2' \ x。如果 M1 被省略或为空,则不应用预条件处理。预条件系统理论上等价于将 bicg 方法应用于线性系统 inv (M1) * A * inv (M2) * y = inv (M1) * binv (M2') * A' * inv (M1') * z = inv (M2') * b,然后设置 x = inv (M2) * y
  • x0 是初始猜测。如果 x0 被省略或为空,则函数默认将 x0 设为零向量。

x0 之后的任何参数都被视为参数,并以适当的方式传递给已提供给 bicg 的任何函数(AfcnMfcn)。

输出参数为:

  • xA * 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*xA'*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 使用预条件子矩阵 M1M2

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, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM, 2003.

另请参阅: 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 * xAfcn 的附加参数在 x0 之后传递。
  • b 是右侧向量。它必须是列向量,且行数与 A 相同。
  • tol 是残差误差所需的相对容差,即 b - A * x。当 norm (b - A * x)tol * norm (b) 时迭代停止。如果 tol 被省略或为空,则使用 1e-6 的容差。
  • maxit 是外部迭代的最大次数,如果未给定或设置为 [],则使用默认值 min (20, numel (b))
  • M1M2 是预条件子。预条件子 M 定义为 M = M1 * M2M1M2 都可以作为矩阵、函数句柄或内联函数 g 传递,使得 g(x) = M1 \ xg(x) = M2 \ x。使用的是右预条件技术,即求解 A * inv (M) * y = b 然后 x = inv (M) * y
  • x0 是初始猜测,如果未给定或设置为 [],则使用默认值 zeros (size (b))

x0 之后的参数被视为参数,并以适当的方式传递给传递给 bicstab 的任何函数(AM)。

输出参数为:

  • 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 使用预条件子矩阵 M1M2

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, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM, 2003.

另请参阅: 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 * xAfcn 的附加参数在 x0 之后传递。
  • b 是右侧向量。它必须是列向量,且行数与 A 相同。
  • tol 是相对容差,如果未给定或设置为 [],则使用默认值 1e-6。
  • maxit 是外部迭代的最大次数,如果未给定或设置为 [],则使用默认值 min (20, numel (b))
  • M1M2 是预条件子。预条件子矩阵定义为 M = M1 * M2M1M2 都可以作为矩阵、函数句柄或内联函数 g 传递,使得 g(x) = M1 \ xg(x) = M2 \ x。如果 M1 为空或未传递,则不应用预条件。使用的是右预条件技术,即求解 A*inv(M)*y = b 然后 x = inv(M)*y
  • x0 是初始猜测,如果未给定或设置为 [],则使用默认值 zeros (size (b))

x0 之后的参数被视为参数,并以适当的方式传递给传递给 cgs 的任何函数(AP)。

输出参数为:

  • 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 使用预条件子矩阵 M1M2

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, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM, 2003.

另请参阅: 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, …)

使用带重启的预条件 GMRES 迭代方法求解 A x = b,也称为 PGMRES(restart)。

输入参数为:

  • A 是线性系统的矩阵,必须是方阵。A 可以作为矩阵、函数句柄或内联函数 Afcn 传递,使得 Afcn(x) = A * xAfcn 的附加参数在 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 就是最大迭代次数。如果 restartmaxit 都不为空,则最大迭代次数为 restart * maxit。如果 restartmaxit 都为空,则最大迭代次数设置为 min (10, N)
  • M1M2 是预条件子。预条件子 M 定义为 M = M1 * M2M1M2 都可以作为矩阵、函数句柄或内联函数 g 传递,使得 g(x) = M1 \ xg(x) = M2 \ x。如果 M1 为 [] 或未给定,则不应用预条件。使用的是左预条件技术,即求解 inv(M) * A * x = inv(M) * b 而不是 A * x = b
  • x0 是初始猜测,如果未给定或设置为 [],则使用默认值 zeros (size (b))

x0 之后的参数被视为参数,并以适当的方式传递给传递给 gmres 的任何函数(AMM1M2)。

输出为:

  • 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 使用预条件子矩阵 M1M2

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, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM, 2003.

另请参阅: 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*xf(x, "transp") = A'*x

预条件子 P 定义为 P = M1 * M2M1M2 都可以作为矩阵、函数句柄或内联函数 g 传递,使得 g(x, "notransp") = M1 \ xg(x, "notransp") = M2 \ x,以及 g(x, "transp") = M1' \ xg(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, 60, pp. 315–339, 1991.
  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, …)

使用基于 cgs 的转置树 qmr(Transpose-Tree qmr)方法求解 A x = b

输入参数为:

  • A 是线性系统的矩阵,必须是方阵。A 可以作为矩阵、函数句柄或内联函数 Afcn 传递,使得 Afcn(x) = A * xAfcn 的附加参数在 x0 之后传递。
  • b 是右侧向量。它必须是列向量,且行数与 A 相同。
  • tol 是相对容差,如果未给定或设置为 [],则使用默认值 1e-6。
  • maxit 是外部迭代的最大次数,如果未给定或设置为 [],则使用默认值 min (20, numel (b))。为了兼容性,由于该方法在迭代次数为奇数或偶数时行为不同,tfqmr 将整个奇偶循环视为一次迭代。也就是说,为了完成一次完整的迭代,算法执行两个子迭代:奇数迭代和偶数迭代。
  • M1M2 是预条件子。预条件子 M 定义为 M = M1 * M2M1M2 都可以作为矩阵、函数句柄或内联函数 g 传递,使得 g(x) = M1 \ xg(x) = M2 \ x。使用的是右预条件技术,即求解 A*inv(M)*y = b 然后 x = inv(M)*y
  • x0 是初始猜测,如果未给定或设置为 [],则使用默认值 zeros (size (b))

x0 之后的参数被视为参数,并以适当的方式传递给传递给 tfqmr 的任何函数(AM)。

输出参数为:

  • 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 使用预条件子矩阵 M1M2

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

示例 6: tfqmr 使用函数作为预条件子

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, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM, 2003.

另请参阅: bicg, bicgstab, cgs, gmres, pcg, qmr.