专用解算器(GNU Octave(版本10.2.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,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 。如果规范(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,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是用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,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(重启)。

输入参数为:

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 : iteration converged to within the specified tolerance

1 : maximum number of iterations exceeded

2 : the preconditioner matrix is singular

3 : algorithm reached stagnation (the relative difference between two

连续迭代小于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,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是包含每次迭代时的残差范数的向量。

参考文献:

R.Freund和N.Nachtigal,QMR: a quasi-minimal residual method for non-Hermitian linear systems,Numeriche Mathematik,1991,60,第315–339页。
R.Barrett、M.Berry、T.Chan、J.Demmel、J.Donato、J.栋加拉、V.Eijkhour、R.Pozo、C.Romine和H.van der Vorst,Templates for the solution of linear systems: Building blocks for iterative methods,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使用基于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,2003年第二版,SIAM

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