专用解算器(GNU Octave(版本9.3.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*标准(b) 如果tol则使用1e-6的公差。
maxit是允许的最大迭代次数;如果maxit则使用值20。
M1,M2是预处理器。预处理器M给定为M=M1 * M2二者都M1和M2可以作为矩阵或函数句柄或内联函数传递g使得g(x,“nottransp”)=M1\x 或g(x,“nottransp”)=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*范数(x2..
- 4:算法无法继续,因为中间值太小或太大,无法进行可靠的计算。
relres是最终残差与其初始值的比率,以欧几里得范数测量。
iter是迭代x计算。
resvec是包含每次迭代时的残差的向量。执行的迭代总数从下式给出长resvec1..

考虑一个三对角矩阵的平凡问题

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 。如果规范(b-A * x)≤tol*标准(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初始猜测,如果未给定或设置为[]默认值零(大小(b))使用。

以下参数x0被视为参数,并以适当的方式传递给任何函数(A或M)传递给bicstab.

输出参数为:

x是计算的近似值。如果该方法不能转换,则它是具有最小残差的迭代方法。
flag表示退出状态:
- 0:迭代收敛到所选公差内
- 1:收敛前达到最大迭代次数
- 2:预处理矩阵是奇异的
- 3:算法达到停滞
- 4:从于被零除,算法无法继续
relres是用as获得的相对残差(A*x-b) /标准b).
iter是(可能是一半)迭代x已计算。如果它是半迭代,那么它是iter+ 0.5
resvec是一个向量,包含每一半和总迭代的残差(也有一半迭代,因为x在每次迭代中以两个步骤计算)。正在执行长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(My.
x0初始猜测,如果未给定或设置为[]默认值zeros (size (b))使用。

以下参数x0被视为参数,并以适当的方式传递给任何函数(A或P)传递给cgs.

输出参数为:

x是计算的近似值。如果该方法不能转换,则它是具有最小残差的迭代方法。
flag表示退出状态:
- 0:迭代收敛到所选公差内
- 1:收敛前达到最大迭代次数
- 2:预处理矩阵是奇异的
- 3:算法达到停滞
- 4:从于被零除,算法无法继续
relres是用as获得的相对残差(A*x-b) /标准b).
iter是迭代x计算。
resvec是包含每次迭代时的残差的向量。正在执行长resvec1.可以看到执行的迭代总数。

让我们考虑一个三对角矩阵的平凡问题

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(重启)。

输入参数为:

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是初始猜测,如果未给定或设置为[],则为默认值零(大小(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次迭代)的预处理相对残差的向量标准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,稀疏线性系统的迭代方法,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:非Hermitian线性系统的一个拟极小残差方法,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,线性系统解的模板:迭代方法的构建块,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) /标准b).
iter是迭代x已计算。
resvec是一个向量,包含每次迭代时的残差(包括norm (b - A x0)).正在执行长resvec1.可以看到执行的迭代总数。

让我们考虑一个三对角矩阵的平凡问题

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.