123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106 |
- // Gauss-Jordan elimination with full pivoting.
- //
- // Uses:
- // (1) solving systems of linear equations (AX=B)
- // (2) inverting matrices (AX=I)
- // (3) computing determinants of square matrices
- //
- // Running time: O(n^3)
- //
- // INPUT: a[][] = an nxn matrix
- // b[][] = an nxm matrix
- //
- // OUTPUT: X = an nxm matrix (stored in b[][])
- // A^{-1} = an nxn matrix (stored in a[][])
- // returns determinant of a[][]
- #include <iostream>
- #include <vector>
- #include <cmath>
- using namespace std;
- const double EPS = 1e-10;
- typedef vector<int> VI;
- typedef double T;
- typedef vector<T> VT;
- typedef vector<VT> VVT;
- T GaussJordan(VVT &a, VVT &b) {
- const int n = a.size();
- const int m = b[0].size();
- VI irow(n), icol(n), ipiv(n);
- T det = 1;
- for (int i = 0; i < n; i++) {
- int pj = -1, pk = -1;
- for (int j = 0; j < n; j++) if (!ipiv[j])
- for (int k = 0; k < n; k++) if (!ipiv[k])
- if (pj == -1 || fabs(a[j][k]) > fabs(a[pj][pk])) { pj = j; pk = k; }
- if (fabs(a[pj][pk]) < EPS) { cerr << "Matrix is singular." << endl; exit(0); }
- ipiv[pk]++;
- swap(a[pj], a[pk]);
- swap(b[pj], b[pk]);
- if (pj != pk) det *= -1;
- irow[i] = pj;
- icol[i] = pk;
- T c = 1.0 / a[pk][pk];
- det *= a[pk][pk];
- a[pk][pk] = 1.0;
- for (int p = 0; p < n; p++) a[pk][p] *= c;
- for (int p = 0; p < m; p++) b[pk][p] *= c;
- for (int p = 0; p < n; p++) if (p != pk) {
- c = a[p][pk];
- a[p][pk] = 0;
- for (int q = 0; q < n; q++) a[p][q] -= a[pk][q] * c;
- for (int q = 0; q < m; q++) b[p][q] -= b[pk][q] * c;
- }
- }
- for (int p = n-1; p >= 0; p--) if (irow[p] != icol[p]) {
- for (int k = 0; k < n; k++) swap(a[k][irow[p]], a[k][icol[p]]);
- }
- return det;
- }
- int main() {
- const int n = 4;
- const int m = 2;
- double A[n][n] = { {1,2,3,4},{1,0,1,0},{5,3,2,4},{6,1,4,6} };
- double B[n][m] = { {1,2},{4,3},{5,6},{8,7} };
- VVT a(n), b(n);
- for (int i = 0; i < n; i++) {
- a[i] = VT(A[i], A[i] + n);
- b[i] = VT(B[i], B[i] + m);
- }
-
- double det = GaussJordan(a, b);
-
- // expected: 60
- cout << "Determinant: " << det << endl;
- // expected: -0.233333 0.166667 0.133333 0.0666667
- // 0.166667 0.166667 0.333333 -0.333333
- // 0.233333 0.833333 -0.133333 -0.0666667
- // 0.05 -0.75 -0.1 0.2
- cout << "Inverse: " << endl;
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < n; j++)
- cout << a[i][j] << ' ';
- cout << endl;
- }
-
- // expected: 1.63333 1.3
- // -0.166667 0.5
- // 2.36667 1.7
- // -1.85 -1.35
- cout << "Solution: " << endl;
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++)
- cout << b[i][j] << ' ';
- cout << endl;
- }
- }
|