1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253 |
- // Reduced row echelon form via Gauss-Jordan elimination
- // with partial pivoting. This can be used for computing
- // the rank of a matrix.
- // Running time: O(n^3)
- // INPUT: a[][] = an nxn matrix
- // OUTPUT: rref[][] = an nxm matrix (stored in a[][])
- // returns rank of a[][]
- typedef double T;
- typedef vector<T> VT;
- typedef vector<VT> VVT;
- int rref(VVT &a) {
- int n = a.size();
- int m = a[0].size();
- int r = 0;
- for (int c = 0; c < m; c++) {
- int j = r;
- for (int i = r+1; i < n; i++)
- if (fabs(a[i][c]) > fabs(a[j][c])) j = i;
- if (fabs(a[j][c]) < EPS) continue;
- swap(a[j], a[r]);
- T s = 1.0 / a[r][c];
- for (int j = 0; j < m; j++) a[r][j] *= s;
- for (int i = 0; i < n; i++) if (i != r) {
- T t = a[i][c];
- for (int j = 0; j < m; j++) a[i][j] -= t * a[r][j];
- }
- r++;
- }
- return r;
- }
- int main() {
- const int n = 5;
- const int m = 4;
- double A[n][m] = { {16,2,3,13},{5,11,10,8},{9,7,6,12},{4,14,15,1},{13,21,21,13} };
- VVT a(n);
- for (int i = 0; i < n; i++)
- a[i] = VT(A[i], A[i] + n);
- int rank = rref (a);
- // expected: 4
- cout << "Rank: " << rank << endl;
- // expected: 1 0 0 1
- // 0 1 0 3
- // 0 0 1 -3
- // 0 0 0 2.78206e-15
- // 0 0 0 3.22398e-15
- cout << "rref: " << endl;
- for (int i = 0; i < 5; i++){
- for (int j = 0; j < 4; j++)
- cout << a[i][j] << ' ';
- cout << endl;
- }
- }
|