Merge branch 'Gathros-new_gauss_c'

Author: leios

contents/gaussian_elimination/code/c/gaussian_elimination.c

@@ -2,7 +2,9 @@
 #include <stdlib.h>
 #include <math.h>
 
-void swap_rows(double * const a, size_t i, size_t pivot, size_t cols) {
+void swap_rows(double *a, const size_t i, const size_t pivot,
+               const size_t cols) {
+
     for (size_t j = 0; j < cols; ++j) {
         double tmp = a[i * cols + j];
         a[i * cols + j] = a[pivot * cols + j];
@@ -11,40 +13,42 @@ void swap_rows(double * const a, size_t i, size_t pivot, size_t cols) {
 }
 
 void gaussian_elimination(double *a, const size_t rows, const size_t cols) {
-    size_t min_dim = (rows < cols)? rows: cols;
+    size_t row = 0;
 
-    for (size_t k = 0; k < min_dim; ++k) {
-        size_t pivot = k;
+    for (size_t col = 0; col < cols - 1; ++col) {
+        size_t pivot = row;
 
-        for (size_t i = k + 1; i < rows; ++i) {
-            if (fabs(a[i * cols + k]) > fabs(a[pivot * cols + k])) {
+        for (size_t i = row + 1; i < rows; ++i) {
+            if (fabs(a[i * cols + col]) > fabs(a[pivot * cols + col])) {
                 pivot = i;
             }
         }
 
-        if (a[pivot * cols + k] == 0) {
+        if (a[pivot * cols + col] == 0) {
             printf("The matrix is singular.\n");
-            exit(0);
+            continue;
         }
 
-        if (k != pivot) {
-            swap_rows(a, k, pivot, cols);
+        if (col != pivot) {
+            swap_rows(a, col, pivot, cols);
         }
 
-        for (size_t i = k + 1; i < rows; ++i) {
-            double scale = a[i * cols + k] / a[k * cols + k];
+        for (size_t i = row + 1; i < rows; ++i) {
+            double scale = a[i * cols + col] / a[row * cols + col];
 
-            for (size_t j = k + 1; j < cols; ++j) {
-                a[i * cols + j] -= a[k * cols + j] * scale;
+            for (size_t j = col + 1; j < cols; ++j) {
+                a[i * cols + j] -= a[row * cols + j] * scale;
             }
 
-            a[i * cols + k] = 0;
+            a[i * cols + col] = 0;
         }
+
+        row++;
     }
 }
 
-void back_substitution(const double * const a, double * const x,
-                       const size_t rows, const size_t cols) {
+void back_substitution(const double *a, double *x, const size_t rows,
+                       const size_t cols) {
 
     for (int i = rows - 1; i >= 0; --i) {
         double sum = 0.0;
@@ -57,13 +61,46 @@ void back_substitution(const double * const a, double * const x,
     }
 }
 
+void gauss_jordan(double *a, const size_t rows, const size_t cols) {
+    int row = 0;
+
+    for (int col = 0; col < cols - 1; ++col) {
+        if (a[row * cols + col] != 0) {
+            for (int i = cols - 1; i > col - 1; --i) {
+                a[row * cols + i] /= a[row * cols + col];
+            }
+
+            for (int i = 0; i < row; ++i) {
+                for (int j = cols - 1; j > col - 1; --j) {
+                    a[i * cols + j] -= a[i * cols + col] * a[row * cols + j];
+                }
+            }
+
+            row++;
+        }
+    }
+}
+
 int main() {
     double a[3][4] = {{3.0, 2.0, -4.0, 3.0},
                       {2.0, 3.0, 3.0, 15.0},
                       {5.0, -3.0, 1.0, 14.0}};
 
     gaussian_elimination((double *)a, 3, 4);
 
+    printf("Gaussian elimination:\n");
+    for (size_t i = 0; i < 3; ++i) {
+        printf("[");
+        for (size_t j = 0; j < 4; ++j) {
+            printf("%f ", a[i][j]);
+        }
+        printf("]\n");
+    }
+
+    printf("\nGauss-Jordan:\n");
+
+    gauss_jordan((double *)a, 3, 4);
+
     for (size_t i = 0; i < 3; ++i) {
         printf("[");
         for (size_t j = 0; j < 4; ++j) {
@@ -72,7 +109,7 @@ int main() {
         printf("]\n");
     }
 
-    printf("\n");
+    printf("\nSolutions are:\n");
 
     double x[3] = {0, 0, 0};
     back_substitution((double *)a, x, 3, 4);

contents/gaussian_elimination/gaussian_elimination.md

@@ -355,7 +355,7 @@ In code, this looks like:
 {% sample lang="jl" %}
 [import:1-45, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-[import:13-44, lang:"c_cpp"](code/c/gaussian_elimination.c)
+[import:15-48, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}
@@ -387,8 +387,7 @@ Here it is in code:
 {% sample lang="jl" %}
 [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-This code does not exist yet in C, so here's Julia code (sorry for the inconvenience)
-[import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
+[import:64-82, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 This code does not exist yet in rust, so here's Julia code (sorry for the inconvenience)
 [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
@@ -419,7 +418,7 @@ In code, this involves keeping a rolling sum of all the values we substitute in
 {% sample lang="jl" %}
 [import:47-67, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-[import:46-58, lang:"c_cpp"](code/c/gaussian_elimination.c)
+[import:50-62, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}