Dev - C++ - ELIMINACIÖN GAUSSIANA

HOla, estoy realizando mi tarea de programación, debo imlementar la clase sistema lineal
realizar la sobrecarga de operadores <<;>>
aplicar el metodo de gauss para resolver el sistema
he intentado lo sigueinte
lo que hace mi programa es lo siguiente:
primero ingreso los valores de A y de el vector b;
luego
tomo el primer elemento de la matriz
y verifico si es ==0 o !=0;
si lo es intercambio filas por una que no tenga el primer elemento cero
caso contrario
los elementos bajo este numero q es el elemento pivote
los hago cero mediante operaciones como multiplicar filas por valores a propiados y restando
pero eso solo lo hago en la primera caoluman como puedo cambiarlo para hacer esto en cada paso bajo cada pivote para que asi me quede un sistema triangular superior

he pensado realizar una función que en cada paso me extraiga submatrices tambien en cada paso igual mente volver a aplicar lo que tengo verifica r si el pivote es cero o nopero escirbi algo y no pude definir tal función en base a lo que tengo:

saludos


#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
using namespace std;

class sistema_lineal
{
private:
int n;
double** A;
double* b;
double* x;
public:
sistema_lineal(double** B, double* c, int m)
{
n=m;
A=new double* [n];
for (int i=0;i<n;i++)
A[i]=new double [n];
for (int i=0;i<n;i++)
for (int j=0;j<n;j++)
A[i][j]=B[i][j];
b=new double [n];
for (int i=0;i<n;i++)
b[i]=c[i];
x=new double[n];
}
sistema_lineal(int m)
{
n=m;
A=new double* [n];
for (int i=0;i<n;i++)
A[i]=new double [n];
for (int i=0;i<n;i++)
for (int j=0;j<n;j++) // creamos la matriz identidad
if (i==j)
A[i][j]=1;
else
A[i][j]=0;
b=new double [n];
for (int i=0;i<n;i++) // vector b de unos
b[i]=1;
x=new double [n];
}
friend istream& operator >>(istream& is, sistema_lineal S);
friend ostream& operator <<(ostream& os, sistema_lineal S);
void resolver_gauss();
void escribir_solucion();
void extraer_matriz();
};

istream& operator >>(istream& is, sistema_lineal S)
{
cout<<"Ingrese la matriz A :"<<endl;
for (int i=0;i<S.n;i++)
for (int j=0;j<S.n;j++)
is>>S.A[i][j];
cout<<"Ingrese el vector b :"<<endl;
for (int i=0;i<S.n;i++)
is>>S.b[i];
return is;
}

ostream& operator <<(ostream& os, sistema_lineal S)
{
for (int i=0;i<S.n;i++)
{
os<<"|\t";
for (int j=0;j<S.n;j++)
os<<S.A[i][j]<<"\t";
os<<"|\t|\tx"<<i+1<<"\t|\t"<<"|\t"<<S.b[i]<<"\t|"<<endl;
}
return os;
}

void sistema_lineal::resolver_gauss()
{double z,w;
if(A[0][0]==0){ //AQUI verifica si el primer elemento de la matriz es distinto
//de cero si lo es cambia la fila por otra
for(int i=1;i<n;i++)
{w=b[i];
b[i]=b[0];
b[0]=w;
if(A[i][0]!=0)
{
for(int j=0;j<=n;j++)
{
z=A[i][j];

A[i][j]=A[0][j];
A[0][j]=z;
}}}
}
else // sino hace ceros debajo del elemento pivote
//pero esto solo lo hace en la primera columna
//quiero que lo haga en cada columna
for(int i=1;i<n;i++){
double m=A[i][0]/A[0][0];
b[i]=b[i]-(m*b[0]);
for(int j=0;j<n;j++)

{
A[i][j]=A[i][j]-(m*A[0][j]);

}}}
void sistema_lineal::extraer_matriz(){
for(int j=0;j<n-1;j++)
for(int i=0;i<n-1;j++){

{A[i][j]=A[i+1][j+1];
}}
}
void sistema_lineal::escribir_solucion()
{
for (int i=0;i<n;i++)
cout<<"x"<<i+1<<" = "<<b[i]<<endl;
}

int main()
{
int n;
cout << "Ingrese n :" <<endl;
cin>>n;
sistema_lineal S1(n);
cin>>S1;
cout<<endl<<endl;
cout<<S1;
cout<<endl<<endl;
S1.resolver_gauss();
cout<<S1;
S1.escribir_solucion();

return 0;
}