扫除算法（sweep）求逆矩阵

扫除算法（sweep）求逆矩阵

分类： 数值计算 2012-11-21 14:25 225人阅读 评论(0) 收藏 举报

高斯约旦消去法是求矩阵逆的一种常用算法，但使用计算机来求解时，需要开辟另一个内存存放变化时的右矩阵A^-1。

而扫除变化时对高斯约旦消去法的一种改进，它可以节省存储空间，不需要开辟内存单独存放右矩阵。下面使用三阶例子来说明：

使用高斯约旦消去法，当左边的矩阵变为单位矩阵时，右边就得到原来左边矩阵的逆矩阵。如果第1步变第1列为基向量，得

当然，第1步不一定要变第1列为基向量，如果变第2列为基向量得

可以看到，不论变哪列都有如下结论：

1.当左矩阵中的某些列未变成基向量时，右矩阵相应列的基向量形式不变。即左右两矩阵中基向量的总数恒为m，并构成一个单位矩阵。

2.左矩阵的第k列变成基向量后，右矩阵第k列才变成非基向量形式，并且该列除第k个元素为1/a(kk)外，其他元素与左矩阵的第k行相差一个符号。

由消去过程得知，当左矩阵的某些列一旦被变成基向量后，就不再起实际作用。为节省计算机的存储单元，左边单位矩阵的基向量不再存储，而把消去过程中右边刚要出现的非基向量放到做左矩阵刚变成基向量的那列位置上，并且这列元素与左矩阵的其他元素一样参加其后不运算。

其后各步消去也做类似处理，扫除变换公式：（设a(kk)不为零）

程序代码：

[csharp] view plain copy

1. Matrix<Type> sweep(Matrix<Type> &A)  
2. {  
3.     Type d;  
4.     int n = A.rows();  
5.   
6.     for (int k = 1; k <= n; ++k)  
7.     {  
8.         d = 1.0/A(k, k);  
9.         A(k, k) = d;  
10.         for (int i = 1; i <= n; ++i)  
11.         {  
12.             if ( i != k )  
13.                 A(k, i) *= (Type) -d;   
14.         }  
15.         for (int i = 1; i <= n; ++i)  
16.         {  
17.             if ( i != k)  
18.                 A(i, k) *= (Type) d;  
19.         }  
20.         for (int i = 1; i <= n; ++i)  
21.         {  
22.             if ( i != k )  
23.             {  
24.                 for (int j = 1; j <= n; ++j)  
25.                 {  
26.                     if ( j != k )  
27.                         A(i, j) += A(i, k)*A(k, j)/d;   
28.                 }  
29.             }  
30.         }  
31.     }  
32.     return A;  
33. }