Magolor 首页文章幻灯片归档标签
从高斯消元到"QGXZ分解"
|线性代数
字数总计: 7.8k|阅读时长: 41 分钟

学习线性代数心血来潮代码一波。

本文将简单介绍矩阵利用高斯消元完成的各种操作和科技，并用C++不使用任何第三方库手写一个封装有矩阵类型和多种操作的库。

主要内容包括: 矩阵的基本运算，矩阵的转置(transpose)，高斯消元，矩阵的echelon form和reduced echelon form，矩阵的迹(trace)，矩阵的秩(rank)，矩阵的逆(inverse)，矩阵的行列式(determinant)，高斯消元解线性方程组，线性方程组的通解，LU分解，PLU分解，以及由于我孤陋寡闻并没有查到的、大佬@GXZlegend传授的"GXZ分解"和"QGXZ分解"。

这里代码的实现更多采用了类似于"函数式编程"的模式，因此实现过程并没有很在意常数，常数可能比较大，使用时请谨慎。

反正我现在不搞OI了，暂时不用关心常数了，啦啦啦！

UPD 2019.10.25: 结尾处GXZ大佬本尊补充了另一种实现方式的"QGXZ分解"。更加OI，封装性相对较低，适合简便、快速地使用。

UPD 2019.10.27: 这种方法类似于LDU分解，两者形式上似乎不完全一样。


矩阵的基本运算和转置

加法、乘法、数乘、转置。似乎没啥好讲的。

(A+B)ij=Aij+Bij(cA)ij=cAij(AB)ij=∑k=1A.mAikBkj(At)ij=Aji(A+B)_{ij} = A_{ij} + B_{ij} \\ (cA)_{ij} = cA_{ij} \\ (AB)_{ij} = \sum_{k=1}^{A.m} A_{ik}B_{kj} \\ (A^t)_{ij} = A_{ji}

下面定义基本的矩阵类型，封装了上述运算，以及矩阵的读入、输出、切割、右拼接。

特别地，这里给出了矩阵快速幂的重载运算符^，然而，为了通用性，快速幂的负数次方被定义为A−K=A.inv()KA^{-K} = A.\text{inv()}^K，目前这里的A.inv()尚未定义。


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#define db double
#define EPS 1e-6
#define X first
#define Y second
#define rint register int
inline bool Zero(db x){return fabs(x)<=EPS;}
struct Mat{db c[128][128]; int n,m; db _, d; inline void clear(){for(rint i = 0; i < n; fill(c[i],c[i]+m,0), i++);} Mat(){} Mat(int a, int b){n=a,m=b,clear();} inline db& operator ()(int i, int j){return c[i][j];}
    inline void read_fix(){for(rint i = 0, j; i < n; i++) for(j = 0; j < m; cin >> c[i][j], j++);} inline void read(){cin >> n >> m; read_fix();}
    inline void print(){for(rint i = 0, j; i < n; i++) for (j = 0; j < m; printf("% 3.3lf%c",c[i][j],"\n "[j<~-m]), j++);}
    inline Mat cut(int lx, int ly, int rx, int ry){Mat A(rx-lx,ry-ly); for(rint i = lx, j; i < rx; i++) for(j = ly; j < ry; A(i-lx,j-ly) = c[i][j], j++); return A;}
    inline Mat merge(Mat B){assert(n==B.n); Mat A(*this); for(rint i = 0, j; i < n; i++) for(j = 0; j < B.m; A(i,m+j) = B(i,j), j++); A.m += B.m; return A;}
    inline Mat transpose(){Mat A(m,n); for(rint i = 0, j; i < n; i++) for(j = 0; j < m; A(j,i) = c[i][j], j++); return A;}
};
inline Mat operator + (Mat A, Mat B){assert(A.n==B.n&&A.m==B.m); for(rint i = 0, j; i < A.n; i++) for(j = 0; j < A.m; A(i,j) += B(i,j), j++); return A;}
inline Mat operator - (Mat A, Mat B){assert(A.n==B.n&&A.m==B.m); for(rint i = 0, j; i < A.n; i++) for(j = 0; j < A.m; A(i,j) -= B(i,j), j++); return A;}
inline Mat operator * (db c, Mat A){for(rint i = 0, j; i < A.n; i++) for(j = 0; j < A.m; A(i,j) *= c, j++); return A;}
inline Mat operator * (Mat A, Mat B){assert(A.m==B.n); Mat C(A.n,B.m); for(rint i = 0, j, k; i < A.n; i++) for(k = 0; k < A.m; k++) if(!Zero(A(i,k))) for(j = 0; j < B.m; C(i,j) += A(i,k)*B(k,j), j++); return C;}
inline Mat& operator += (Mat &A, Mat B){return A=A+B;} inline Mat& operator -= (Mat &A, Mat B){return A=A-B;} inline Mat& operator *= (Mat &A, Mat B){return A=A*B;} inline Mat& operator *= (Mat &A, db c){return A=c*A;}
inline Mat operator ^ (Mat S, int K){assert(S.n==S.m); if(K<0) S = S.inv(), K = -K; Mat A(S.n,S.n); for(rint i = 0; i < A.n; A(i,i) = 1, i++); for(; K; S *= S, K >>= 1) if(K&1) A *= S; return A;}

OI选手可能会表示不满: 怎么都是double啊？

取模请自己手动魔改…反正我现在不搞OI了，啦啦啦

显然，本文不会涉及辗转相除高斯消元。感兴趣者可以在OI资源分享中找到相关内容。


高斯消元，echelon form和reduced echelon form

初等行变换: 交换两行，一行乘等于一个常数，一行加等于若干倍的另一行。

高斯消元: 通过初等行变换将矩阵变成某个想要的形式。

高斯消元似乎也没啥好讲的，基本操作吧。大体思路就是模仿人手计算方式: 将除了第一个式子以外所有含有第一个变量的式子消掉第一个变量，将除了第二个式子以外所有含有第二个变量的式子消掉第二个变量…以此类推。

矩阵的(row) echelon form: 第ii行最左侧的非零元素的位置在第i+1i+1行最左侧的非零元素的位置的严格左侧。记P(i)P(i)表示第ii行最左侧的非零的下标，即P(i)<P(i+1),∀i∈{0,1,2,...,n−2}P(i) \lt P(i+1), \forall i \in \lbrace 0,1,2,...,n-2 \rbrace(下标从0开始). AA的(row) echelon form以后简记为AEA_{E}.

矩阵的reduced (row) echelon form: 在 (row) echelon form的基础上满足Ai,P(i)=1,Ak,P(i)=0,∀i∈{0,1,2,...,n−2},k≠iA_{i,P(i)} = 1, A_{k,P(i)} = 0, \forall i \in \lbrace 0,1,2,...,n-2 \rbrace, k \ne i. AA的reduced (row) echelon form以后简记为AREA_{RE}.


高斯消元得到echelon form(只正消不反消):

1
2
3
4
5
6
inline pair<Mat,Mat> echelon(Mat B=Mat(0,0)){ assert(n==B.n||!B.n); Mat A(*this); A.d = 1;
    for(rint k = 0, t = 0, i, j; t < n && k < m; k++){
        if(Zero(A(t,k))) for(i = t+1; i < n; i++) if(!Zero(A(i,k))){for(j = k; j < m; swap(A(t,j),A(i,j)), j++); for(j = 0; j < B.m; swap(B(t,j),B(i,j)), j++); A.d = -A.d; break;}
        if(Zero(A(t,k))) continue; for(i = t+1; i < n; i++){for(_ = A(i,k)/A(t,k), j = k; j < m; A(i,j) -= _*A(t,j), j++); for(j = 0; j < B.m; B(i,j) -= _*B(t,j), j++);} t++;
    }   return make_pair(A,B);
}

这里传入矩阵BB表示在将AA消成echelon form的同时对BB进行完全相同的操作。因为初等行变换等价于左乘一个矩阵，这里即分块矩阵的乘法分配: E(A∣B)=(EA∣EB)E(A \vert B) = (EA \vert EB)。调用A.echelon(B)，返回一个pair<Mat,Mat>，即(EA,EB)(EA,EB)。


高斯消元得到reduced echelon form(先正消后反消):

1
2
3
4
5
6
7
inline pair<Mat,Mat> gauss(Mat B=Mat(0,0)){ assert(n==B.n||!B.n); Mat A(*this);
    for(rint k = 0, t = 0, i, j; t < n && k < m; k++){
        if(Zero(A(t,k))) for(i = t+1; i < n; i++) if(!Zero(A(i,k))){for(j = k; j < m; swap(A(t,j),A(i,j)), j++); for(j = 0; j < B.m; swap(B(t,j),B(i,j)), j++); break;}
        if(Zero(A(t,k))) continue; for(j = 0; j < B.m; B(t,j) /= A(t,k), j++); for(j = m-1; j >= k; A(t,j) /= A(t,k), j--);
        for(i = !t; i < n; i++, i += i==t){for(_ = A(i,k)/A(t,k), j = k; j < m; A(i,j) -= _*A(t,j), j++); for(j = 0; j < B.m; B(i,j) -= _*B(t,j), j++);} t++;
    }   return make_pair(A,B);
}

同理，调用A.gauss(B)，返回一个pair<Mat,Mat>，即(EA,EB)(EA,EB)。


这里说明一下: double高斯消元容易产生精度误差，一个优化的方法是每次选列主元的时候选择绝对值最大的。这里因为懒没有采用这种优化。


矩阵的迹、秩和逆

只有方阵才有定义。

矩阵的迹就是对角线上元素的和。这玩意有啥用？

消成echelon form以后非零行的数目就是矩阵的秩(rank)。

nn阶方阵AA可逆当且仅当rank(A)=n\text{rank}(A) = n。

A−1(A∣I)=(I∣A−1)A^{-1}(A \vert I) = (I \vert A^{-1}). 对于一个可逆方阵来说，I=AREI = A_{RE}. 因此只要将AA捆绑着单位阵，消成AREA_{RE}，那么单位阵就会成为AA的逆。即A.inv()就是A.gauss(I).second。


1
2
3
inline db trace(){assert(n==m); _ = 0; for(rint i = 0; i < n; _ += c[i][i], i++); return _;}
inline int rank(){Mat A(this->echelon().X); int r = 0; for(rint i = 0; i < n; r += count_if(A.c[i],A.c[i]+m,Zero)<m, i++); return r;}
inline Mat inv(){assert(n==m&&rank()==n); Mat I(n,n); for(rint i = 0; i < n; I(i,i) = 1, i++); return this->gauss(I).Y;}

矩阵的行列式

只有方阵才有定义。

这里为了方便，暂时使用从11开始的下标:

∀i,det⁡A=∑j=1n(−1)i+jAijdet⁡A~ij\forall i, \det A = \sum_{j=1}^n (-1)^{i+j}A_{ij} \det \tilde A_{ij}

这里A~ij\tilde A_{ij}表示删除第ii行和第jj列后的矩阵。

这样不方便求，但是一行加等于另一行的变换不改变行列式、交换两行的变换使行列式取反，因此可以先高斯消元再求行列式。

而下或上三角矩阵行列式很好求，每次只考虑i=1i = 1:

det⁡A=∑j=1n(−1)1+jA1jdet⁡A~1j=A11det⁡A~11=A11A22det⁡A~22=...=∏i=1nAii\begin{aligned} \det A &= \sum_{j=1}^n (-1)^{1+j}A_{1j} \det \tilde A_{1j} = A_{11} \det \tilde A_{11} \\ &= A_{11}A_{22} \det \tilde A_{22} = ... = \prod_{i=1}^n A_{ii} \end{aligned}

因此先转成echelon form，然后直接对角线乘积就好了。注意一下符号。


1
inline db det(){assert(n==m); Mat A(this->echelon().X); _ = 1; for(rint i = 0; i < n; _ *= A(i,i), i++); return A.d*_;}

线性方程组

一个线性方程组可以表示为矩阵形式Ax=bAx = b，其中xx和bb都是列向量。

有E(A∣b)=(ARE∣b′)E(A \vert b) = (A_{RE} \vert b')。而AREx=b′A_{RE}x = b'则是非常好解的: Ai,P(i)=1A_{i,P(i)} = 1，因此xP(i)=bi′−∑j>P(i)xjx_{P(i)} = b'_i - \sum_{j \gt P(i)} x_j, 其余的xkx_k(∀i,k≠P(i)\forall i, k \ne P(i))可以任取。


这里支持可视化输出线性方程组的通解: 给定已经转化为reduced echelon form的矩阵AREA_{RE}和bb，输出所有可能的xx。

1
2
3
4
5
6
7
inline void solution_set(pair<Mat,Mat> S){ int n = S.X.n, m = S.X.m; Mat A(S.X), b(S.Y); assert(n==b.n&&1==b.m);
    for(rint i = 0; i < n; i++) if(count_if(A.c[i],A.c[i]+m,Zero)==m&&!Zero(b(i,0))){puts("No Solution."); return;}
    for(rint j = 0, i, k; j < m; j++) for(i = 0; i < n; i++) if(!Zero(A(i,j))){
        if(count_if(A.c[i],A.c[i]+m,Zero) < j+count_if(A.c[i]+j,A.c[i]+m,Zero)){printf("x%d is free.\n",j); break;}
        printf("x%d = % 3.3lf",j,b(i,0)); for(k = j+1; k < m; k++) if(!Zero(A(i,k))) printf(" %c %3.3lfx%d","+-"[A(i,k)>0],fabs(A(i,k)),k); puts("."); break;
    }
}

函数传入一个pair<Mat,Mat>即(ARE,b)(A_{RE},b)，输出可能的xx的所有解(自由元独立输出，其余变量用自由元表示)，具体形式可以自行试验，文末有试验。


LU分解和PLU分解

注意到每次解线性方程组复杂度都是高斯消元的复杂度O(n3)O(n^3)。

如果每次bb更换，而AA不变，有没有更好的方法呢？

对于一个可逆方阵来说，Ax=bAx = b即x=A−1bx = A^{-1}b。可以直接求逆，求逆O(n3)O(n^3)，此后每次计算A−1bA^{-1}b只需要O(n2)O(n^2)(矩乘列向量)。

然而对于非方阵怎么做呢？

U.P.D: 对于离线，可以直接构造(A∣b1∣b2...∣bq)(A \vert b_1 \vert b_2 ... \vert b_q)，一次高斯消元解决。

LU分解是这样一种思路，设A=LUA = LU，则:

Ax=b⇔LUx=b⇔{Ly=bUx=yAx = b \Leftrightarrow LUx = b \Leftrightarrow \begin{cases} Ly = b \\ Ux = y \\ \end{cases}

如果可以找到一一组合适的LL和UU，使得后两个方程很好解(O(n2)O(n^2))，就做完了。

LU分解发现，LL为n×nn \times n下三角矩阵，U=AEU = A_E的时候似乎会更好解。

怎么分解A=LU=LAEA = LU = LA_{E}呢？只需要高消就行了: E(A∣I)=(AE∣E)E(A \vert I) = (A_E \vert E)，而A=LU=LAE=LEAA = LU = LA_{E} = LEA，整体即L=E−1L = E^{-1}，因此可以捆绑高消求出EE然后直接求逆即可。


PLU分解几乎没有区别，只是为了方便: 如果分解A=LUA = LU过程中使用了交换操作，那么就会出现顺序问题，可以分解PA=LUPA = LU，其中PP是置换矩阵。

注: 并未代码实现过PLU分解，不做担保。

由于尚未实现PLU，因此在LU分解过程中必须保证AA不会发生交换。这是一个比较苛刻的条件。

留坑待填

考虑在做的时候发生了置换，那么这时候求出来的不再是A=LUA = LU，而是A=(P−1L)UA = (P^{-1}L)U，其中PP是一个置换矩阵。

既然LL是一个下三角，那么P−1LP^{-1}L通过选择排序就可以O(n2)O(n^2)地转换为一个下三角。那么P(P−1L∣I)=(L∣P)P(P^{-1}L \vert I) = (L \vert P)，就求出了PP和LL。

现在已经知道了P,L,UP,L,U，并且有A=P−1LUA = P^{-1}LU，那么Ax=b⇔LUx=PbAx = b \Leftrightarrow LUx = Pb. 只需要每次把bb乘上PP(矩乘列向量当然也是O(n2)O(n^2)的)，然后再代入LU分解的做法即可。

将PbPb看做新的bb，此后只考虑LU分解。


如何计算Ly=bLy = b？

yi=bi−∑k<iAikyky_i = b_i - \sum_{k \lt i} A_{ik}y_k

如何计算Ux=yUx = y？

这里如果要计算通解，就出现了一点问题，因为U=AEU = A_E可能存在很多自由元。例如:

U=[1111111...10111222...20011333...30001444...4]U = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & ... & 1 \\ 0 & 1 & 1 & 1 & 2 & 2 & 2 & ... & 2 \\ 0 & 0 & 1 & 1 & 3 & 3 & 3 & ... & 3 \\ 0 & 0 & 0 & 1 & 4 & 4 & 4 & ... & 4 \\ \end{bmatrix}

这时候首先可以得到x4=y4−4(x5+x6+...+xm)x_4 = y_4 - 4(x_5+x_6+...+x_m). 如果任取自由元都代入00当然没有问题，但是如果要求通解，保留整个公式的话，那么计算x3=y3−x4−3(x5+x6+...+xm)x_3 = y_3 - x_4 - 3(x_5+x_6+...+x_m)，这里的−x4-x_4操作就不再是O(1)O(1)而是O(m)O(m)的了，因为x4x_4本身就是一个无法化简的由多个自由元组成的表达式。

这样复杂度就再次退回到了O(n3)O(n^3)(这里假设n,mn,m同阶)。


“GXZ分解"和"QGXZ分解”

首先不保证以下内容的正确性。其次不保证以下内容的独创性(希望大佬指出这东西叫啥)。

请教了隔壁宿舍大佬@GXZlegend。大佬给出了一种将LU分解的UU矩阵继续拆分来解决这个问题的方法。大致思路就是模仿高斯消元到底——反消。

设A=GXZA = GXZ，其中GG是一个n×nn \times n下三角矩阵(就是原来的LL)，XX是一个n×nn \times n上三角矩阵，Z=AREZ = A_{RE}。而其中XZXZ就是原来的UU，GXGX就高斯消元的变换矩阵EE，事实上就是把高斯消元的正消和反消两步拆成了两个矩阵。

现在就有:

GXZx=b⇔{Gz=bXy=zZx=yGXZx = b \Leftrightarrow \begin{cases} Gz = b \\ Xy = z \\ Zx = y \\ \end{cases}

其中解Gz=bGz = b和Xy=zXy = z的方法和解原来Ly=bLy = b相同，而现在Z=AREZ = A_{RE}，因此要输出通解——

当然是直接调用solution_set函数啊！


1
inline void GXZ(Mat &G, Mat &X, Mat &Z){Mat I(n,n); for(rint i = 0; i < n; I(i,i) = 1, i++); pair<Mat,Mat> S(this->echelon(I)); G = S.Y.inv(); S = S.X.gauss(I); X = S.Y.inv(); Z = S.X;}
1
2
3
4
inline void GXZsolve(Mat G, Mat X, Mat Z, Mat b){ int n = G.n; assert(n==b.n&&1==b.m); Mat z(n,1), y(n,1), x(Z.m,1); db _;
    for(rint i = 0, j; i < n; i++) for(_ = G(i,i), z(i,0) = b(i,0)/_, j = 0; j < i; z(i,0) -= G(i,j)/_*z(j,0), j++);
    for(rint i = n-1, j; ~i; i--) for(_ = X(i,i), y(i,0) = z(i,0)/_, j = n-1; j > i; y(i,0) -= X(i,j)/_*y(j,0), j--); solution_set(make_pair(Z,y));
}

QGXZ分解其实就是PLU分解的思路应用到GXZ分解上。为什么换成了字母Q？，因为GXZ大佬太强了!

首先解出A=(Q−1G)XZA = (Q^{-1}G)XZ，然后选择排序Q(Q−1G∣I)=(G∣Q)Q(Q^{-1}G \vert I) = (G \vert Q)，于是就有:

Q−1GXZx=b⇔{Gz=QbXy=zZx=yQ^{-1}GXZx = b \Leftrightarrow \begin{cases} Gz = Qb \\ Xy = z \\ Zx = y \\ \end{cases} 


1
2
3
4
5
inline void QGXZ(Mat &Q, Mat &G, Mat &X, Mat &Z){
    Q = Mat(n,n); for(rint i = 0; i < n; Q(i,i) = 1, i++); GXZ(G,X,Z);
    for(rint j = n-1, i, k; ~j; j--) for(i = 0; i <= j; i++) if(!Zero(G(i,j)))
        {for(k = 0; k < n; swap(G(i,k),G(j,k)), swap(Q(i,k),Q(j,k)), k++); break;}
}

1
inline void QGXZsolve(Mat Q, Mat G, Mat X, Mat Z, Mat b){GXZsolve(G,X,Z,Q*b);}

U.P.D: "QGXZ分解"的另一种实现

@GXZlegend大佬亲自实现的QGXZ分解。说明见[线性代数xOI/ACM]系数矩阵的QGXZ分解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include <bits/stdc++.h>
#define N 510
#define eps 1e-6
using namespace std;
int pos[N];
double Q[N][N] , G[N][N] , X[N][N] , Z[N][N] , b[N] , c[N] , d[N] , e[N];
int main()
{
    int n , m , q , i , j , k , p = 0 , t;
    double mx;
    scanf("%d%d%d" , &n , &m , &q);
    for(i = 1 ; i <= n ; i ++ )
        for(j = 1 ; j <= m ; j ++ )
            scanf("%lf" , &Z[i][j]);
    for(i = 1 ; i <= n ; i ++ ) Q[i][i] = G[i][i] = X[i][i] = 1;
    for(i = 1 ; i <= m ; i ++ )
    {
        t = 0 , mx = eps;
        for(j = p + 1 ; j <= n ; j ++ )
            if(abs(Z[j][i]) > mx)
                t = j , mx = abs(Z[j][i]);
        if(!t) continue;
        pos[ ++ p] = i;
        for(k = i ; k <= m ; k ++ ) swap(Z[p][k] , Z[t][k]);
        for(k = 1 ; k <= n ; k ++ ) swap(Q[p][k] , Q[t][k]);
        for(k = 1 ; k < p ; k ++ ) swap(G[p][k] , G[t][k]);
        for(j = p + 1 ; j <= n ; j ++ )
        {
            G[j][p] = Z[j][i] / Z[p][i];
            for(k = i ; k <= m ; k ++ )
                Z[j][k] -= Z[p][k] * G[j][p];
        }
    }
    for(i = p ; i ; i -- )
    {
        for(j = i - 1 ; j ; j -- )
        {
            X[j][i] = Z[j][pos[i]] / Z[i][pos[i]];
            for(k = pos[i] ; k <= m ; k ++ )
                Z[j][k] -= Z[i][k] * X[j][i];
        }
    }
    while(q -- )
    {
        for(i = 1 ; i <= n ; i ++ ) scanf("%lf" , &b[i]);
        for(i = 1 ; i <= n ; i ++ )
            for(j = 1 ; j <= n ; j ++ )
                if(Q[i][j] == 1)
                    c[j] = b[i];
        for(i = 1 ; i <= n ; i ++ )
        {
            d[i] = c[i];
            for(j = 1 ; j < i ; j ++ )
                d[i] -= G[i][j] * d[j];
        }
        for(i = n ; i ; i -- )
        {
            e[i] = d[i];
            for(j = n ; j > i ; j -- )
                e[i] -= X[i][j] * e[j];
        }
        for(i = p + 1 ; i <= n ; i ++ )
            if(abs(e[i]) > eps)
                break;
        if(i <= n) puts("No solution!");
        else
        {
            for(i = 1 ; i <= p ; i ++ )
            {
                printf("x[%d]=%lf" , pos[i] , e[i] / Z[i][pos[i]]);
                for(j = pos[i] + 1 ; j <= m ; j ++ )
                    if(abs(Z[i][j]) > eps)
                        printf("%+lfx[%d]" , -Z[i][j] / Z[i][pos[i]] , j);
                puts("");
            }
        }
    }
    return 0;
}

完整库以及测试文件

未经DEBUG慎用！！！

未经DEBUG慎用！！！

未经DEBUG慎用！！！

完整库:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <bits/stdc++.h>
using namespace std;
#define db double
#define EPS 1e-6
#define X first
#define Y second
#define rint register int
inline bool Zero(db x){return fabs(x)<=EPS;}
struct Mat{db c[128][128]; int n,m; db _, d; inline void clear(){for(rint i = 0; i < n; fill(c[i],c[i]+m,0), i++);} Mat(){} Mat(int a, int b){n=a,m=b,clear();} inline db& operator ()(int i, int j){return c[i][j];}
    inline void read_fix(){for(rint i = 0, j; i < n; i++) for(j = 0; j < m; cin >> c[i][j], j++);} inline void read(){cin >> n >> m; read_fix();}
    inline void print(){for(rint i = 0, j; i < n; i++) for (j = 0; j < m; printf("% 3.3lf%c",c[i][j],"\n "[j<~-m]), j++);}
    inline Mat cut(int lx, int ly, int rx, int ry){Mat A(rx-lx,ry-ly); for(rint i = lx, j; i < rx; i++) for(j = ly; j < ry; A(i-lx,j-ly) = c[i][j], j++); return A;}
    inline Mat merge(Mat B){assert(n==B.n); Mat A(*this); for(rint i = 0, j; i < n; i++) for(j = 0; j < B.m; A(i,m+j) = B(i,j), j++); A.m += B.m; return A;}
    inline Mat transpose(){Mat A(m,n); for(rint i = 0, j; i < n; i++) for(j = 0; j < m; A(j,i) = c[i][j], j++); return A;}
    inline pair<Mat,Mat> echelon(Mat B=Mat(0,0)){ assert(n==B.n||!B.n); Mat A(*this); A.d = 1;
        for(rint k = 0, t = 0, i, j; t < n && k < m; k++){
            if(Zero(A(t,k))) for(i = t+1; i < n; i++) if(!Zero(A(i,k))){for(j = k; j < m; swap(A(t,j),A(i,j)), j++); for(j = 0; j < B.m; swap(B(t,j),B(i,j)), j++); A.d = -A.d; break;}
            if(Zero(A(t,k))) continue; for(i = t+1; i < n; i++){for(_ = A(i,k)/A(t,k), j = k; j < m; A(i,j) -= _*A(t,j), j++); for(j = 0; j < B.m; B(i,j) -= _*B(t,j), j++);} t++;
        }   return make_pair(A,B);
    }
    inline pair<Mat,Mat> gauss(Mat B=Mat(0,0)){ assert(n==B.n||!B.n); Mat A(*this);
        for(rint k = 0, t = 0, i, j; t < n && k < m; k++){
            if(Zero(A(t,k))) for(i = t+1; i < n; i++) if(!Zero(A(i,k))){for(j = k; j < m; swap(A(t,j),A(i,j)), j++); for(j = 0; j < B.m; swap(B(t,j),B(i,j)), j++); break;}
            if(Zero(A(t,k))) continue; for(j = 0; j < B.m; B(t,j) /= A(t,k), j++); for(j = m-1; j >= k; A(t,j) /= A(t,k), j--);
            for(i = !t; i < n; i++, i += i==t){for(_ = A(i,k)/A(t,k), j = k; j < m; A(i,j) -= _*A(t,j), j++); for(j = 0; j < B.m; B(i,j) -= _*B(t,j), j++);} t++;
        }   return make_pair(A,B);
    }
    inline int rank(){Mat A(this->echelon().X); int r = 0; for(rint i = 0; i < n; r += count_if(A.c[i],A.c[i]+m,Zero)<m, i++); return r;}
    inline Mat inv(){assert(n==m&&rank()==n); Mat I(n,n); for(rint i = 0; i < n; I(i,i) = 1, i++); return this->gauss(I).Y;}
    inline db det(){assert(n==m); Mat A(this->echelon().X); _ = 1; for(rint i = 0; i < n; _ *= A(i,i), i++); return A.d*_;}
    inline db trace(){assert(n==m); _ = 0; for(rint i = 0; i < n; _ += c[i][i], i++); return _;}
    inline void GXZ(Mat &G, Mat &X, Mat &Z){Mat I(n,n); for(rint i = 0; i < n; I(i,i) = 1, i++); pair<Mat,Mat> S(this->echelon(I)); G = S.Y.inv(); S = S.X.gauss(I); X = S.Y.inv(); Z = S.X;}
    inline void QGXZ(Mat &Q, Mat &G, Mat &X, Mat &Z){
        Q = Mat(n,n); for(rint i = 0; i < n; Q(i,i) = 1, i++); GXZ(G,X,Z);
        for(rint j = n-1, i, k; ~j; j--) for(i = 0; i <= j; i++) if(!Zero(G(i,j)))
            {for(k = 0; k < n; swap(G(i,k),G(j,k)), swap(Q(i,k),Q(j,k)), k++); break;}
    }
};
inline Mat operator + (Mat A, Mat B){assert(A.n==B.n&&A.m==B.m); for(rint i = 0, j; i < A.n; i++) for(j = 0; j < A.m; A(i,j) += B(i,j), j++); return A;}
inline Mat operator - (Mat A, Mat B){assert(A.n==B.n&&A.m==B.m); for(rint i = 0, j; i < A.n; i++) for(j = 0; j < A.m; A(i,j) -= B(i,j), j++); return A;}
inline Mat operator * (db c, Mat A){for(rint i = 0, j; i < A.n; i++) for(j = 0; j < A.m; A(i,j) *= c, j++); return A;}
inline Mat operator * (Mat A, Mat B){assert(A.m==B.n); Mat C(A.n,B.m); for(rint i = 0, j, k; i < A.n; i++) for(k = 0; k < A.m; k++) if(!Zero(A(i,k))) for(j = 0; j < B.m; C(i,j) += A(i,k)*B(k,j), j++); return C;}
inline Mat& operator += (Mat &A, Mat B){return A=A+B;} inline Mat& operator -= (Mat &A, Mat B){return A=A-B;} inline Mat& operator *= (Mat &A, Mat B){return A=A*B;} inline Mat& operator *= (Mat &A, db c){return A=c*A;}
inline Mat operator ^ (Mat S, int K){assert(S.n==S.m); if(K<0) S = S.inv(), K = -K; Mat A(S.n,S.n); for(rint i = 0; i < A.n; A(i,i) = 1, i++); for(; K; S *= S, K >>= 1) if(K&1) A *= S; return A;}
inline void solution_set(pair<Mat,Mat> S){ int n = S.X.n, m = S.X.m; Mat A(S.X), b(S.Y); assert(n==b.n&&1==b.m);
    for(rint i = 0; i < n; i++) if(count_if(A.c[i],A.c[i]+m,Zero)==m&&!Zero(b(i,0))){puts("No Solution."); return;}
    for(rint j = 0, i, k; j < m; j++) for(i = 0; i < n; i++) if(!Zero(A(i,j))){
        if(count_if(A.c[i],A.c[i]+m,Zero) < j+count_if(A.c[i]+j,A.c[i]+m,Zero)){printf("x%d is free.\n",j); break;}
        printf("x%d = % 3.3lf",j,b(i,0)); for(k = j+1; k < m; k++) if(!Zero(A(i,k))) printf(" %c %3.3lfx%d","+-"[A(i,k)>0],fabs(A(i,k)),k); puts("."); break;
    }
}
inline void GXZsolve(Mat G, Mat X, Mat Z, Mat b){ int n = G.n; assert(n==b.n&&1==b.m); Mat z(n,1), y(n,1), x(Z.m,1); db _;
    for(rint i = 0, j; i < n; i++) for(_ = G(i,i), z(i,0) = b(i,0)/_, j = 0; j < i; z(i,0) -= G(i,j)/_*z(j,0), j++);
    for(rint i = n-1, j; ~i; i--) for(_ = X(i,i), y(i,0) = z(i,0)/_, j = n-1; j > i; y(i,0) -= X(i,j)/_*y(j,0), j--); solution_set(make_pair(Z,y));
}
inline void QGXZsolve(Mat Q, Mat G, Mat X, Mat Z, Mat b){GXZsolve(G,X,Z,Q*b);}
#undef db
#undef EPS
#undef X
#undef Y
#undef rint

测试文件:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <bits/stdc++.h>
#include "mat.h" 
using namespace std;
#define db double
#define EPS 1e-6
#define X first
#define Y second
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}
int main(){
    Mat A(3,4);
    A(0,0) = 1; A(0,1) = 1; A(0,2) = 0; A(0,3) = 0;
    A(1,0) = 1; A(1,1) = 0; A(1,2) = 1; A(1,3) = 1;
    A(2,0) = 2; A(2,1) = 1; A(2,2) = 0; A(2,3) = 0;
    A.print(); puts("");
    
    A.echelon().X.print(); puts("");
    
    A.transpose().print(); puts("");
    
    A.transpose().echelon().X.print(); puts("");
    
    A.cut(0,1,3,2).print(); puts("");
    
    A.merge(A.cut(0,1,3,2)).print(); puts("");
    
    Mat B(3,1);
    B(0,0) = 1;
    B(1,0) = 2;
    B(2,0) = 3;
    B.print(); puts("");
    auto S = A.gauss(B);
    S.X.print(); puts("");
    S.Y.print(); puts("");
    solution_set(S); puts("");
    
    Mat C(3,4);
    C(0,0) = 1; C(1,1) = 1; C(2,3) = 1;
    C.echelon().X.print();
    cout << C.rank() << endl; puts("");
    
    Mat D(3,3);
    D(0,0) = 3; D(0,1) = 2; D(0,2) = 1;
    D(1,0) = 5; D(1,1) = 4; D(1,2) = 3;
    D(2,0) = 7; D(2,1) = 3; D(2,2) = 1;
    D.print(); puts("");
    D.echelon().X.print();
    cout << D.trace() << endl;
    cout << D.det() << endl; puts("");
    
    (D^3).print(); puts("");
    (D*D*D).print(); puts("");
    (D^-1).print(); puts("");
    D.inv().print(); puts("");
    (D*D.inv()).print(); puts("");
    (D.inv()*D).print(); puts("");
    
    Mat E(2,2);
    E(1,0) = 1; E(1,1) = 1;
//  (E.inv()).print(); puts("");    //E is not invertible!!!!
    
    Mat Q, G, X, Z;
    A.GXZ(G,X,Z);
    G.print(); puts("");
    X.print(); puts("");
    Z.print(); puts("");
    (G*X*Z).print(); puts("");
    GXZsolve(G,X,Z,B); puts("");
    
    Mat F(3,4); 
    F(0,0) = 0; F(0,1) = 1; F(0,2) = 0; F(0,3) = 0;
    F(1,0) = 0; F(1,1) = 0; F(1,2) = 1; F(1,3) = 3;
    F(2,0) = 2; F(2,1) = 1; F(2,2) = 3; F(2,3) = 2;
    F.print(); puts("");
    
    F.QGXZ(Q,G,X,Z);
    Q.print(); puts("");
    G.print(); puts("");
    X.print(); puts("");
    Z.print(); puts("");
    (Q.inv()*G*X*Z).print(); puts("");
    QGXZsolve(Q,G,X,Z,B); puts("");
    
    E.QGXZ(Q,G,X,Z);
    Q.print(); puts("");
    G.print(); puts("");
    X.print(); puts("");
    Z.print(); puts("");
    (Q.inv()*G*X*Z).print(); puts("");
    Mat H(2,1);
    H(0,0) = 0;
    H(1,0) = 2;
    QGXZsolve(Q,G,X,Z,H); puts("");
    
    return 0;
}

输出结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
 1.000  1.000  0.000  0.000
 1.000  0.000  1.000  1.000
 2.000  1.000  0.000  0.000

 1.000  1.000  0.000  0.000
 0.000 -1.000  1.000  1.000
 0.000  0.000 -1.000 -1.000

 1.000  1.000  2.000
 1.000  0.000  1.000
 0.000  1.000  0.000
 0.000  1.000  0.000

 1.000  1.000  2.000
 0.000 -1.000 -1.000
 0.000  0.000 -1.000
 0.000  0.000  0.000

 1.000
 0.000
 1.000

 1.000  1.000  0.000  0.000  1.000
 1.000  0.000  1.000  1.000  0.000
 2.000  1.000  0.000  0.000  1.000

 1.000
 2.000
 3.000

 1.000  0.000  0.000  0.000
 0.000  1.000  0.000  0.000
 0.000  0.000  1.000  1.000

 2.000
-1.000
-0.000

x0 =  2.000.
x1 = -1.000.
x2 = -0.000 - 1.000x3.
x3 is free.

 1.000  0.000  0.000  0.000
 0.000  1.000  0.000  0.000
 0.000  0.000  0.000  1.000
3

 3.000  2.000  1.000
 5.000  4.000  3.000
 7.000  3.000  1.000

 3.000  2.000  1.000
 0.000  0.667  1.333
 0.000  0.000  2.000
8
4

 233.000  150.000  87.000
 483.000  312.000  181.000
 393.000  253.000  147.000

 233.000  150.000  87.000
 483.000  312.000  181.000
 393.000  253.000  147.000

-1.250  0.250  0.500
 4.000 -1.000 -1.000
-3.250  1.250  0.500

-1.250  0.250  0.500
 4.000 -1.000 -1.000
-3.250  1.250  0.500

 1.000 -0.000  0.000
-0.000  1.000  0.000
-0.000 -0.000  1.000

 1.000 -0.000 -0.000
 0.000  1.000  0.000
 0.000 -0.000  1.000

 1.000  0.000  0.000
 1.000  1.000  0.000
 2.000  1.000  1.000

 1.000  1.000  0.000
 0.000 -1.000  1.000
-0.000 -0.000 -1.000

 1.000  0.000  0.000  0.000
 0.000  1.000  0.000  0.000
 0.000  0.000  1.000  1.000

 1.000  1.000  0.000  0.000
 1.000  0.000  1.000  1.000
 2.000  1.000  0.000  0.000

x0 =  2.000.
x1 = -1.000.
x2 = -0.000 - 1.000x3.
x3 is free.

 0.000  1.000  0.000  0.000
 0.000  0.000  1.000  3.000
 2.000  1.000  3.000  2.000

 0.000  0.000  1.000
 1.000  0.000  0.000
 0.000  1.000  0.000

 1.000  0.000  0.000
 0.000  1.000  0.000
 0.000  0.000  1.000

 2.000  1.000  3.000
 0.000  1.000  0.000
 0.000  0.000  1.000

 1.000  0.000  0.000 -3.500
 0.000  1.000  0.000  0.000
 0.000  0.000  1.000  3.000

 0.000  1.000  0.000  0.000
 0.000  0.000  1.000  3.000
 2.000  1.000  3.000  2.000

x0 = -2.000 + 3.500x3.
x1 =  1.000.
x2 =  2.000 - 3.000x3.
x3 is free.

 0.000  1.000
 1.000  0.000

 1.000  0.000
 0.000  1.000

 1.000  0.000
 0.000  1.000

 1.000  1.000
 0.000  0.000

 0.000  0.000
 1.000  1.000

x0 =  2.000 - 1.000x1.
x1 is free.

总结

GXZx = Qb！！！

%%%%%%


扫描二维码即可在手机上查看这篇文章，或者转发二维码来分享这篇文章:


文章作者: Magolor
文章链接: https://magolor.cn/2019/10/24/2019-10-24-blog-01/
版权声明: 本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Magolor！
OI资源黑科技线性代数高斯消元
扫描二维码在手机上查看或转发二维码以分享Magolor的博客