// Gaussian Elimination #include #include #include // Gaussian Elimination: Transform a Matrix to Row Echelon Form // using elementary transformations of matrix rows. // Return the rank of matrix (i.e. the number of nonzero rows // in Row Echelon Form of matrix) // Parameters: // double *a -- array of matrix elements in row major order // int m -- number of matrix rows // int n -- number of matrix columns // double eps -- precision of algorithm (default eps = 1e-12) int gauss(double *a, int m, int n, double eps = 1e-12); // Input/output a matrix bool readMatrix(FILE* f, double* &a, int& m, int& n); bool readMatrix(const char* path, double* &a, int& m, int& n); bool writeMatrix(FILE* f, const double* a, int m, int n); bool writeMatrix(const char* path, const double* a, int m, int n); // Determinant of a matrix in Row Echelon Form double refDet(const double* a, int m, int n); // Reverse pass after the gaussian elimination. // Transform echelon form of matrix to the // Reduced Row Echelon Form // using elementary transformations of matrix rows. // Example: // Initial Matrix Row Echelon Form // 1 2 3 4 13.0000 14.0000 15.0000 16.0000 // 5 6 7 8 --> 0.0000 -0.9231 -1.8462 -2.7692 // 9 10 11 12 0.0000 0.0000 0.0000 0.0000 // 13 14 15 16 0.0000 0.0000 0.0000 0.0000 // // Reduced Row Echelon Form // 1.0000 0.0000 -1.0000 -2.0000 // --> 0.0000 1.0000 2.0000 3.0000 // 0.0000 0.0000 0.0000 0.0000 // 0.0000 0.0000 0.0000 0.0000 void reversePass(double *a, int m, int n, int rank, double eps = 1e-12); int main() { double *a; int m, n; if (!readMatrix("input.txt", a, m, n)) { perror("Could not read the matrix"); delete[] a; return -1; } writeMatrix(stdout, a, m, n); // Print the initial maatrix int rank = gauss(a, m, n); printf("rank of matrix: %d\n", rank); printf("Row Echelon Form of matrix:\n"); writeMatrix(stdout, a, m, n); // Print the matrix writeMatrix("output.txt", a, m, n); // and write in into the file if (m == n) { double d = refDet(a, m, n); printf("Determinant: %f\n", d); } reversePass(a, m, n, rank); printf("Reverse Pass (Reduced Row Echelon Form):\n"); writeMatrix(stdout, a, m, n); delete[] a; return 0; } int gauss(double *a, int m, int n, double eps) { int i = 0; int j = 0; while (i < m && j < n) { // Search for the maximal element in j-th column double maxelem = fabs(a[i*n + j]); int k = i; for (int l = i+1; l < m; ++l) { if (fabs(a[l*n + j]) > maxelem) { maxelem = fabs(a[l*n + j]); k = l; } } if (maxelem <= eps) { ++j; continue; } // printf( // "max element %f in %d-th row and %d-th column\n", // maxelem, k, j // ); if (k != i) { // Swap i-th and k-th rows for (int s = j; s < n; ++s) { double tmp = a[i*n + s]; a[i*n + s] = a[k*n + s]; a[k*n + s] = (-tmp); } // printf("Swap rows %d and %d\n", k, i); } // writeMatrix(stdout, a, m, n); // Zero the j-th column of matrix starting from i+1-th row double r = a[i*n + j]; // Resolving element for (int l = i+1; l < m; ++l) { double c = a[l*n + j]/r; // Subtract i-th row multiplied by coeff. c from the l-th row a[l*n + j] = 0.; for (int s = j+1; s < n; ++s) { a[l*n + s] -= a[i*n + s]*c; } } // printf("Annihilated %d-th column\n", j); // writeMatrix(stdout, a, m, n); ++i; ++j; } return i; } bool readMatrix(FILE* f, double* &a, int& m, int& n) { a = 0; if (fscanf(f, "%d%d", &m, &n) < 2) return false; a = new double[m*n]; for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { if (fscanf(f, "%lg", &(a[i*n + j])) < 1) { delete[] a; a = 0; return false; } } } return true; } bool readMatrix(const char *path, double* &a, int& m, int& n) { FILE *f = fopen(path, "rt"); if (f == NULL) return false; bool res = readMatrix(f, a, m, n); fclose(f); return res; } bool writeMatrix(FILE* f, const double* a, int m, int n) { fprintf(f, "%d %d\n", m, n); for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { if (j > 0) fprintf(f, " "); // Delimiter fprintf(f, "%12.4f", a[i*n + j]); } fprintf(f, "\n"); // Delimiter } return true; } bool writeMatrix(const char *path, const double* a, int m, int n) { FILE *f = fopen(path, "wt"); if (f == NULL) return false; bool res = writeMatrix(f, a, m, n); fclose(f); return res; } double refDet(const double* a, int m, int n) { assert(m == n); double d = a[0]; for (int i = 1; i < m; ++i) { d *= a[i*n + i]; } return d; } void reversePass(double *a, int m, int n, int rank, double eps) { assert(rank <= m); int i = rank - 1; while (i >= 0) { // Search for the first nonzero element in i-th rows int j = 0; while (j < n && fabs(a[i*n + j]) <= eps) ++j; if (j >= n) assert(false); int j0 = j; double r = a[i*n + j0]; assert(fabs(r) > eps); while (j < n) { a[i*n + j] /= r; ++j; } // Annihilate j-th column in rows 0..i-1 for (int k = 0; k < i; ++k) { r = a[k*n + j0]; a[k*n + j0] = 0.; for (j = j0 + 1; j < n; ++j) { a[k*n + j] -= a[i*n + j]*r; } } --i; } }