专用解算器(GNU Octave(版本11.1.0))

前一个矩阵的函数上级线性代数[目录][索引]

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 。如果规范(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, 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 * x。的附加参数Afcn在之后通过x0.
b是右侧向量。它必须是行数与相同的列向量A.
tol是残余误差所需的相对误差范围,b - A * x 。如果规范(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是用as获得的相对残差(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, 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 * 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是用as获得的相对残差(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, 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, …) ¶

解决A x = b使用带重启的预处理GMRES迭代方法,也称为PGMRES(重启)。

输入参数为:

A是线性系统的矩阵,它必须是平方。A可以作为矩阵、函数句柄或内联函数传递Afcn使得Afcn(x) = A * x。的附加参数Afcn在之后通过x0.
b是右侧向量。它必须是行数与相同的列向量A.
restart是方法重新启动之前的迭代次数。如果它是[]或N=numel(b),则不应用restartis。
tol是第二个定位残差所需的相对误差范围,inv (M) * (b - a * x)。迭代操作如果标准(inv (M) * (b - a * x))≤tol*标准(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使用aleft预处理器

[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*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是包含每次迭代时的残差范数的向量。

参考文献:

R. Freund and N. Nachtigal, "QMR: a quasi-minimal residual method for non-Hermitian linear systems", Numerische Mathematik, 60, pp. 315–339, 1991.
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使用基于cgs的Transpose Tree qmr方法。

输入参数为:

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

详见: bicg,bicgstab,cgs,gmres,pcg,qmr,pcr.