Reading up on numerical methods and c++

main.cpp

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void triangulate(float* a, float* b, int n);
void backsubst(float* a, float* b, int n);
void print(float* a, int rows, int cols);

int main() {
  printf("\nGauss Seidel Elimination - Ch. 10 - Direct Solution of Linear Equation, page 204, An Introduction To Numerical Methods in C++\n\n");

  float a[] = {
    1.0f, 2.0f, 3.0f, 4.0f,
    2.0f, 3.0f, 4.0f, 6.0f,
    3.0f, 4.0f, 2.0f, 5.0f,
    4.0f, 6.0f, 5.0f, 7.0f
  };

  float b[] = { 8.0f, 6.0f, 4.0f, 2.0f };
  float c[4] = { 0 } ;

  triangulate(a, b, 4);
  backsubst(a, b,4);

  print(a, 4, 4);
  print(b, 4, 1);

  return EXIT_SUCCESS;
}


void triangulate(float* a, float* b, int n) {

  for (int j = 0; j < n-1; ++j) { // rows
    float diag = a[j * n + j];
    for(int i = j + 1; i < n; ++i) { 
      float factor = a[i * n + j] / diag;
      b[i] -= b[j] * factor;
      for(int k = 0; k < n; ++k) {
        a[i * n + k] -= a[j * n + k] * factor;
      }
    }
  }
}

void backsubst(float* a, float* b, int n) {

  for(int i = n - 1; i >=0 ; --i) {

    float diag = a[i * n + i];
    if(fabs(diag) < 0.001) {
      printf("Diag too small.\n");
      return;
    }

    float sum = b[i];
    for(int k = i + 1; k < n; ++k) {
      sum -= a[i * n + k] * b[k];
    }

    b[i] = sum / diag;
  }
}


void print(float* a, int rows, int cols) {
  for(int j = 0; j < rows; ++j) {
    for(int i = 0; i < cols; ++i) {
      printf("%2.2f ", a[j * cols + i]);
    }
    printf("\n");
  }
  printf("-\n");
}

run.sh

#!/bin/sh
g++ -g main.cpp -o gauss
./gauss