# ilu

## 语法

``[L,U] = ilu(A)``
``[L,U,P] = ilu(A)``
``W = ilu(A)``
``[___] = ilu(A,options)``

## 说明

``[L,U] = ilu(A)` 用零填充执行稀疏矩阵 `A` 的不完全 LU 分解，并返回下三角矩阵 `L` 和上三角矩阵 `U`。`
``[L,U,P] = ilu(A)` 还返回置换矩阵 `P`，并满足 `L` 和 `U` 是 `P*A` 或 `A*P` 的不完全因子。默认情况下，`P` 是用于不使用主元消去的不完全 LU 分解的单位矩阵。`

``W = ilu(A)` 返回 LU 因子的非零值。输出 `W` 等于 `L + U - speye(size(A))`。`

``[___] = ilu(A,options)` 使用结构体 `options` 指定的选项对 `A` 执行不完全 LU 分解。例如，通过将 `options` 的 `type` 字段设置为 `"ilutp"`，您可以使用主元消去执行不完全 LU 分解。然后，通过将 `milu` 字段设置为 `"row"` 或 `"col"`，可以指定保留修正不完全 LU 分解的行值总和或列值总和。这种选项组合返回置换矩阵 `P`，使得 `L` 和 `U` 是 `"row"` 选项的 `A*P` 的不完全因子（其中 `U` 是列置换的），或 `L` 和 `U` 是 `"col"` 的 `P*A` 的不完全因子（其中 `L` 是行置换的）。`

## 示例

`ilu` 函数提供三种类型的不完全 LU 分解：零填充分解 (`ILU(0)`)、Crout 版本 (`ILUC`)，以及具有阈值调降和主元消去的分解 (`ILUTP`)。

```A = gallery("neumann",1600) + speye(1600); n = nnz(A)```
```n = 7840 ```
`n = nnz(lu(A))`
```n = 126478 ```
`n = nnz(ilu(A))`
```n = 7840 ```

```[L,U] = ilu(A); e = norm(A-(L*U).*spones(A),"fro")/norm(A,"fro")```
```e = 4.8874e-17 ```

`e = norm(A-L*U,"fro")/norm(A,"fro")`
```e = 0.0601 ```

```options.droptol = 1e-6; options.type = "crout"; [L,U] = ilu(A,options);```

`n = nnz(ilu(A,options))`
```n = 51482 ```
`e = norm(A-L*U,"fro")./norm(A,"fro")`
```e = 9.3040e-07 ```

```options.type = "ilutp"; [L,U,P] = ilu(A,options); n = nnz(ilu(A,options))```
```n = 51541 ```
`norm(P*A-L*U,"fro")./norm(A,"fro")`
```ans = 9.4960e-07 ```

```load west0479 A = west0479; c1 = condest(A)```
```c1 = 1.4244e+12 ```

```[P,R,C] = equilibrate(A); B = R*P*A*C; c2 = condest(B)```
```c2 = 5.1036e+04 ```

```options.type = "ilutp"; options.milu = "row"; options.droptol = 0; [L,U,P] = ilu(B,options);```

`e = norm(B*P-L*U,"fro")/norm(B,"fro")`
```e = 1.0345e-16 ```
`nLU = nnz(L)+nnz(U)-size(B,1)`
```nLU = 19118 ```
`nB = nnz(B)`
```nB = 1887 ```

```ntols = 20; tau = logspace(-6,-2,ntols); e = zeros(1,ntols); nLU = zeros(1,ntols); for k = 1:ntols options.droptol = tau(k); [L,U,P] = ilu(B,options); nLU(k) = nnz(L)+nnz(U)-size(B,1); e(k) = norm(B*P-L*U,"fro")/norm(B,"fro"); end figure semilogx(tau,nLU./nB,LineWidth=2) title("Ratio of nonzeros of LU factors with respect to B") xlabel("drop tolerance") ylabel("nnz(L)+nnz(U)-size(B,1)/nnz(B)",FontName="FixedWidth")```

```figure loglog(tau,e,LineWidth=2) title("Relative error of the incomplete LU factorization") xlabel("drop tolerance") ylabel("\$||BP-LU||_F\,/\,||B||_F\$",Interpreter="latex")```

```load west0479 A = west0479;```

`b = full(sum(A,2));`

```tol = 1e-12; maxit = 20;```

• `x` 是计算 `A*x = b` 所得的解。

• `fl0` 是指示算法是否收敛的标志。

• `rr0` 是计算的解 `x` 的相对残差。

• `it0` 是计算 `x` 时所用的迭代序号。

• `rv0`$‖\mathit{b}-\mathit{Ax}‖$ 的由每个二分之一迭代的残差历史记录组成的向量。

```[x,fl0,rr0,it0,rv0] = bicgstab(A,b,tol,maxit); fl0```
```fl0 = 1 ```
`rr0`
```rr0 = 1 ```
`it0`
```it0 = 0 ```

`bicgstab` 未在请求的 20 次迭代内收敛至请求的容差 `1e-12`，因此 `fl0``1`。实际上，`bicgstab` 的行为太差，因此初始估计值 `x0 = zeros(size(A,2),1)` 是最佳解，并会返回最佳解（如 `it0 = 0` 所示）。

```options = struct("type","ilutp","droptol",1e-6); [L,U] = ilu(A,options); [x_precond,fl1,rr1,it1,rv1] = bicgstab(A,b,tol,maxit,L,U); fl1```
```fl1 = 0 ```
`rr1`
```rr1 = 5.9892e-14 ```
`it1`
```it1 = 3 ```

```semilogy((0:numel(rv0)-1)/2,rv0/norm(b),"-o") hold on semilogy((0:numel(rv1)-1)/2,rv1/norm(b),"-o") yline(tol,"r--"); legend("No preconditioner","ILU preconditioner","Tolerance",Location="East") xlabel("Iteration number") ylabel("Relative residual")```

## 输入参数

LU 分解选项，指定为包含最多五个字段的结构体：

`type`

`type` 的有效值为：

• `"nofill"`（默认值）- 执行具有 `0` 填充级别的 ILU 分解（称为 ILU(0)）。如果将 `type` 设置为 `"nofill"`，则仅使用 `milu` 选项；所有其他字段都将被忽略。

• `"crout"` - 执行 ILU 分解的 Crout 版本，称为 ILUC。如果将 `type` 设置为 `"crout"`，则仅使用 `droptol``milu` 选项；所有其他字段都将被忽略。

• `"ilutp"` - 执行带阈值和主元消去的 ILU 分解，称为 ILUTP。仅可以对此类型执行主元消去。

`droptol`

`droptol` 的有效值是任何非负标量。`U` 的非零项满足 `abs(U(i,j)) >= droptol*norm(A(:,j))`，但对角线元素除外，无论它们是否满足此标准，对角线元素都将保留。在使用主元调整 `L` 的元之前，将根据局部调降容差检验这些元，以使 `L` 的非零项满足 `abs(L(i,j)) >= droptol*norm(A(:,j))/U(j,j)`

`milu`

`milu` 的有效值为：

• `"off"` - 不生成修正不完全 LU 分解。

• `"row"` - 生成行值总和修正不完全 LU 分解。补偿上三角因子 `U` 的对角线元素，以便保留原始矩阵的行值总和。即 `A*e = L*U*e`，其中 `e` 是由 1 组成的列向量。

• `"col"` - 生成列总和修正不完全 LU 分解。补偿上三角因子 `U` 的对角线元素，以便保留原始矩阵的列总和。即 `e'*A = e'*L*U`

`udiag`

`udiag` 的有效值为 `0``1`。如果 `udiag``1`，则上三角因子 `U` 的任何零对角线元素都将替换为局部调降容差，以尝试避免奇异因子。

`thresh`

## 输出参量

• `options.type` 指定为 `"nofill"`（默认值）或 `"crout"` 时，则 `L` 返回为一个单位下三角矩阵（即，主对角线上元素为 1 的下三角矩阵）。

• `options.type` 指定为 `"ilutp"` 并且 `options.milu` 指定为 `"col"` 时，`L` 返回为一个行置换的单位下三角矩阵。

• `options.type` 指定为 `"nofill"`（默认值）或 `"crout"` 时，`U` 返回为一个上三角矩阵。

• `options.type` 指定为 `"ilutp"` 并且 `options.milu` 指定为 `"row"` 时，`U` 返回为一个列置换的上三角矩阵。

`options.type` 指定为 `"nofill"`（默认值）或 `"crout"` 时，`P` 返回为一个与 `A` 大小相同的单位矩阵，因为这两种类型都不执行主元消去。

LU 因子的非零值，以稀疏方阵形式返回，其中 `W = L + U - speye(size(A))`

## 参考

[1] Saad, Y. "Preconditioning Techniques", chap. 10 in Iterative Methods for Sparse Linear Systems. The PWS Series in Computer Science. Boston: PWS Pub. Co, 1996.