#include "Matrix.h" #include int Matrix::gauss(double eps) { int i = 0; int j = 0; while (i < m && j < n) { // 1. Search for maximal element in j-th column // starting from i-th row double amax = fabs((*this)[i][j]); int k = i; for (int l = i+1; l < m; ++l) { if (fabs((*this)[l][j]) > amax) { amax = fabs((*this)[l][j]); k = l; } } if (amax <= eps) { ++j; continue; } if (k != i) { // Swap rows i and k for (int s = j; s < n; ++s) { double tmp = (*this)[i][s]; (*this)[i][s] = (*this)[k][s]; (*this)[k][s] = (-tmp); } } // Annihilate j-th column starting from i+1-th row for (int k = i+1; k < m; ++k) { double coeff = (*this)[k][j]/(*this)[i][j]; for (int s = j; s < n; ++s) { (*this)[k][s] -= coeff*(*this)[i][s]; } } ++i; ++j; } return i; } int Matrix::rref(double eps) { // Reduced row echelon form int rank = gauss(); int i = rank - 1; while (i >= 0) { // Find the first nonzero element in i-th row int j = 0; while (j < cols() && fabs((*this)[i][j]) <= eps) { ++j; } assert(j < cols()); // Divide the i-th row by the first nonzero element a[i][j] double coeff = 1./(*this)[i][j]; (*this)[i][j] = 1.; for (int k = j+1; k < cols(); ++k) { (*this)[i][k] *= coeff; } // Annihilate the j-th column above i-th row for (int l = 0; l < i; ++l) { double coeff = (*this)[l][j]; (*this)[l][j] = 0.; for (int k = j+1; k < cols(); ++k) { (*this)[l][k] -= (*this)[i][k]*coeff; } } --i; // To the next row in reverse order } return rank; } double Matrix::det() const { if (rows() != cols()) { throw std::invalid_argument( "determinant of non-square matrix" ); } Matrix a = *this; int rank = a.gauss(); if (rank < a.rows()) return 0.; double d = 1.; for (int i = 0; i < a.rows(); ++i) { d *= a[i][i]; } return d; } std::ostream& operator<<(std::ostream& s, const Matrix& a) { for (int i = 0; i < a.rows(); ++i) { for (int j = 0; j < a.cols(); ++j) { if (j > 0) s << ' '; // Delimiter s << a[i][j]; } s << std::endl; } return s; } std::istream& operator>>(std::istream& s, Matrix& a) { for (int i = 0; s.good() && i < a.rows(); ++i) { for (int j = 0; s.good() && j < a.cols(); ++j) { double x; if (s >> x) { a[i][j] = x; } } } return s; }