矩阵算法 高斯消元 行列式 矩阵的秩
今天学习一下矩阵的基本算法
高斯消元是解线性方程组的有力工具。
基本思想是通过将增广矩阵经过行初等变化变成简化阶梯形矩阵。
下面采用的是列主元高斯消元法,复杂度为O(n^3)。
很容易根据高斯消元法的过程得出行列式和秩的算法。
代码:
/********************************************************* * ------------------ * * author AbyssalFish * **********************************************************/ #include<cstdio> #include<iostream> #include<string> #include<cstring> #include<queue> #include<vector> #include<stack> #include<vector> #include<map> #include<set> #include<algorithm> #include<cmath> //#include<bits/stdc++.h> using namespace std; typedef long long ll; const double EPS = 1e-8; typedef vector<double> vec; typedef vector<vec> mat; //O(n^3) vec gauss_jordan(const mat& A, const vec& b) { int n = A.size(); mat B(n,vec(n+1)); //Augment Matrix for(int i = 0; i < n; i++) for(int j = 0; j < n; j++) B[i][j] = A[i][j]; for(int i = 0; i < n; i++) B[i][n] = b[i]; for(int i = 0; i < n; i++){ int piv = i; //取最大以便判断无解或者无穷多解 for(int j = i; j < n; j++){ if(abs(B[j][i]) > abs(B[piv][i])) piv = j; } if(i != piv) swap(B[i],B[piv]); if(abs(B[i][i]) < EPS) return vec(); //假想把系数变成1,只计算对后面有影响的部分 for(int j = n; j > i; j--) B[i][j] /= B[i][i]; for(int j = 0; j < n; j++) if(i != j) { for(int k = i+1; k <= n; k++) B[j][k] -= B[j][i]*B[i][k]; } } vec x(n); for(int i = 0; i < n; i++) x[i] = B[i][n]; return x; } double determinant(const mat& A) { int n = A.size(); mat B = A; double det = 1; int sign = 0; for(int i = 0; i < n; i++){ int piv = i; for(int j = i; j < n; j++){ if(abs(B[j][i]) > abs(B[piv][i])) piv = j; } if(i != piv) swap(B[i],B[piv]), sign ^= 1; if(abs(B[i][i]) < EPS) return 0; det *= B[i][i]; for(int j = i+1; j < n; j++) B[i][j] /= B[i][i]; for(int j = i+1; j < n; j++) { for(int k = i+1; k < n; k++) B[j][k] -= B[j][i]*B[i][k]; } } return sign?-det:det; } int rank_of_mat(const mat& A) { int n = A.size(); mat B = A; for(int i = 0; i < n; i++){ int piv = i; for(int j = i; j < n; j++){ if(abs(B[j][i]) > abs(B[piv][i])) piv = j; } if(i != piv) swap(B[i],B[piv]); if(abs(B[i][i]) < EPS) return i; for(int j = n; --j > i;) B[i][j] /= B[i][i]; for(int j = i+1; j < n; j++) { for(int k = i+1; k < n; k++) B[j][k] -= B[j][i]*B[i][k]; } } return n; } double read(){ double t; scanf("%lf",&t); return t; } //#define LOCAL int main() { #ifdef LOCAL freopen("in.txt","r",stdin); #endif vec alpha; mat A; int n; puts("input a matrix"); while(~scanf("%d",&n) && n <= 0) puts("input invalid"); A.resize(n); for(int i = 0; i < n; i++){ for(int j = 0; j < n; j++){ A[i].push_back(read()); } } puts("input a vector"); for(int i = 0; i < n; i++) alpha.push_back(read()); puts("\nsolution"); vec sol = gauss_jordan(A,alpha); if(sol.size()){ for(int i = 0; i < n; i++) printf("%lf%c", sol[i], i!=n-1?' ':'\n'); }else { puts("not unique"); } puts("\ndeterminant"); printf("%lf\n", determinant(A)); puts("\nrank"); printf("%d\n", rank_of_mat(A)); return 0; }
测试样例
实际使用中如果同时需要这些信息只计算一遍B矩阵,从而提高效率。