Dev C++ из метода простой итерации в метод Гаусса

[Денис999]

Задана система линейных уравнений: MA * X = B где А - матрица размера NxN (N - количество уравнений в системе - задается
пользователем); элементы матрицы заполняются пользователем;
В - вектор размером Nx1 (элементы вектора вводятся пользователем);
X - вектор размером Nx1 - вектор неизвестных величин, относительно которых решается уравнение.
сделал методом простой итерации:

#include <windows.h>
#include <process.h>
#include <stdio.h>
#include <math.h>
#define N 4
#define NumberOfProcessors 2
#define EPSILON 1e-5
float Matrix[N][N + 1] =
{
  { 1.0f, 0.1f, 0.1f, 0.1f, 1.0f },
  { 0.1f, 1.2f, 0.1f, 0.1f, 1.0f },
  { 0.2f, 0.1f, 1.0f, 0.1f, 1.0f },
  { 0.2f, 0.2f, 0.2f, 1.0f, 1.0f }
};
float NewX[N];
float CurX[N] = { 1.0, 1.0, 1.0, 1.0 };
void InitializeMatrix()
{
  printf("Source Matrix\n");
  for (unsigned i = 0; i < N; i++)
  {
    for (unsigned j = 0; j <= N; j++)
    {
      // calculate (E - A)
      if (i == j)
        Matrix[i][j] = 1.0f - Matrix[i][j];
      else
        if (j < N)
          Matrix[i][j] = -Matrix[i][j];
      printf("%7.4f ", Matrix[i][j]);
    }
    printf("\n");
  }
}
enum Commands { Continue, Stop };
struct IterState
{
  HANDLE Starters[NumberOfProcessors];
  HANDLE evIterationFinished[NumberOfProcessors];
};
IterState State;
unsigned Command = Continue;
unsigned __stdcall ThreadFunction(void *pData)
{
  unsigned Branch = (unsigned)pData;
  unsigned nFrom = Branch * N / NumberOfProcessors;
  unsigned nTo   = (Branch + 1) * N / NumberOfProcessors;
  unsigned i;
  for (;;)
  {
    WaitForSingleObject(State.Starters[Branch], INFINITE);
    if (Command == Stop)
      return 0;
    for (i = nFrom; i < nTo; i++)
    {
      NewX[i] = Matrix[i][0] * CurX[0];
      for (unsigned j = 1; j < N; j++)
        NewX[i] += Matrix[i][j] * CurX[j];
      NewX[i] += Matrix[i][N];
    }
    // synchronization
    SetEvent(State.evIterationFinished[Branch]);
    if (Branch == 0)
    {
      WaitForMultipleObjects(NumberOfProcessors,
          State.evIterationFinished, TRUE, INFINITE);
      float d = -1;
      for (i = 0; i < N; i++)
        if (fabs(NewX[i] - CurX[i]) > d)
          d = (float)fabs(NewX[i] - CurX[i]);
      if (d < EPSILON)
        Command = Stop;
      memcpy(CurX, NewX, sizeof(NewX));
      for (i = 0; i < NumberOfProcessors; i++)
        SetEvent(State.Starters[i]);
    }
  }
}
int main(int argc, char* argv[])
{
  UNREFERENCED_PARAMETER(argc);
  UNREFERENCED_PARAMETER(argv);
  
  InitializeMatrix();
  HANDLE hThreads[NumberOfProcessors];
  unsigned i, tid;
  for (i = 0; i < NumberOfProcessors; i++)
  {
    State.Starters[i] = 
      CreateEvent(NULL, FALSE, TRUE, NULL);
    State.evIterationFinished[i] = 
      CreateEvent(NULL, FALSE, FALSE, NULL);
  }
  
  for (i = 0; i < NumberOfProcessors; i++)
    hThreads[i] = (HANDLE)_beginthreadex(
      NULL, 0, ThreadFunction,
      (void *)i, 0, &tid);
  WaitForMultipleObjects(NumberOfProcessors, hThreads, 
    TRUE, INFINITE);
  for (i = 0; i < NumberOfProcessors; i++)
  {
    CloseHandle(State.Starters[i]);
    CloseHandle(State.evIterationFinished[i]);
    CloseHandle(hThreads[i]);
  }
  for (i = 0; i < N; i++)
    printf("x%d = %7.4f\n", i, NewX[i]);
  return 0;
}

Подскажите как ее решить теперь методом Гаусса?