Метода Гаусса с выбором главного элемента по столбцу C/C++

Метод Гаусса достаточно хорошо документирован, поэтому я в очередной раз избавлю себя и вас от теории. К тому же у данного метода существует множество модификаций. Эта, пожалуй, самая простая из всех.


Алгоритм
1) Находим максимальный элемент в первом столбце, затем меняем местами строки L с первой строкой
2) Находим максимальный элемент во втором столбце, меняем L2 со второй строкой и т.д.

Сразу определимся с выделением памяти. Оно будет динамическим. В этой модификации нужно менять местами строки. Куда проще сделать это с помощью трёх операций присваивания, меняя местами адреса, нежели N операций присваивания, цикле перегоняя кучу цифр.
matrix=new float*[n];
for (int i=0; i<n; i++){
matrix[i]=new float[n+1]; 
}
А ещё логичнее держать правую часть (Вектор B) под рукой. Т.е. N+1ым столбцом нашей исходной матрицы. Не работать же с ним отдельно в другом массиве! Много чести будет!
Прямой ход. В результате получим приведённую к треугольному виду матрицу. Заранее извиняюсь за говнокод.
int L=0; //номер строки
float max_el=0.0;
int max_el_row;
float *tmp_address;
float tmpf=0.0, tmpf2=0.0, t3=0.0;
for (int j=0; j<_n; j++)
{
        for (int h=j; h<_n; h++)
        {
    if (fabs(matrix[h][j])>max_el)
            {
            max_el_row=h;
            max_el=fabs(matrix[h][j]);
            }
        }
    max_el=0.0;
    tmp_address=matrix[max_el_row];
    matrix[max_el_row]=matrix[L];
    matrix[L]=tmp_address;
for (int k=1+L; k<_n; k++)  //! k=1; L=0,1,2...
{
    tmpf=matrix[k][j]/matrix[L][j];
    for (int i=0; i<=_n; i++)
    {
        if (fabs(tmpf)==0) break;
        tmpf2=(matrix[L][i]*tmpf);
        t3=matrix[k][i]-tmpf2;
        matrix[k][i]=t3;
    }
}
L++;
}
Обратный ход. В результате будут найдены корни СЛАУ. Хочу напомнить что вектор B хранится у меня в последнем столбце матрицы.
float x[20];
for(int i=_n-1; i>= 0; i--)
{
   double s=0;
        for(int j=_n-1; j > i; j--)
        {
            s+=x[j]*matrix[i][j];
        }
        x[i]=(matrix[i][_n]-s)/matrix[i][i];
        if (x[i]==0) x[i]=fabs(x[i]);
}