Метод Йордана - Гаусса с выбором главного элемента по всей матрице

void            SwapString(double ** &A,int i,int k)
{
    if(i == k) return;
    double * p = A[i];
        A[i] = A[k];
        A[k] = p;
}
 
void            SwapColumn(double ** &A,int n,int k,int l)
{
    if(k == l) return;
 
    for(int i = 0;i < n;i++)
    {
        double t = A[i][k];
        A[i][k] = A[i][l];
        A[i][l] = t;
    }
}
 
 
// A - исходная матрица
// A_inv - матрица которую нужно обратить 
void JS(double ** A,double ** &A_inv,int n)
{
        E(A_inv,n);               
 
    int * st = new int[n];
    int * ct = new int[n];
 
    for(int i = 0;i < n;i++)
    {
        double max = A[i][i];
        int imax = i,jmax = i;
 
        for(int z = i;z < n;z++)
            for(int j = i;j < n;j++)
                if(A[z][j] > max)
                {
                    imax = z;
                    jmax = j;
                    max = A[z][j];
                }
        st[i] = imax;
        ct[i] = jmax;
 
        SwapString(A,imax,i);
        SwapColumn(A,n,jmax,i);
        
        SwapString(A_inv,imax,i);
        SwapColumn(A_inv,n,jmax,i);
        
 
        if(A[i][i] != 0)
        {
            double t = 1./A[i][i];
            for(int j = i;j < n;j++)
            {
                A[i][j] = A[i][j]*t;
            }
 
            for(int j = 0;j < n;j++)
            {
                A_inv[i][j] = A_inv[i][j]*t;
            }
 
            for(int k = 0;k < i;k++)
            {
                double t = A[k][i];
                for(int l = i;l < n;l++)
                {
                    A[k][l] = A[k][l] - A[i][l]*t;
                }
 
                for(int l = 0;l < n;l++)
                {
                    A_inv[k][l] = A_inv[k][l] - A_inv[i][l]*t;
                }
            }
 
            for(int k = i + 1;k < n;k++)
            {
                double t = A[k][i];
                for(int l = i;l < n;l++)
                {
                    A[k][l] = A[k][l] - A[i][l]*t;
                }
 
                for(int l = 0;l < n;l++)
                {
                    A_inv[k][l] = A_inv[k][l] - A_inv[i][l]*t;
                }
            }
        }
    }
    for(int i = n - 1;i >= 0;i--)
    {
        SwapColumn(A_inv,n,st[i],i);
    }
 
    for(int i = n - 1;i >= 0;i--)
    {
        SwapString(A_inv,ct[i],i);
    }
 
    delete [] st;
    delete [] ct; 
} 
 
// A - исходная матрица
// B - столбец правых частей
// для СЛАУ 
 
void JS_SLE(double ** A,double * B,int n)
{
 
//  int * st = new int[n];
    int * ct = new int[n];
 
    for(int i = 0;i < n;i++)
    {
        double max = A[i][i];
        int imax = i,jmax = i;
 
/*      for(int z = i;z < n;z++)                           // по всей матрице
         for(int j = i;j < n;j++)
           if(A[z][j] > max)
           {
            imax = z;
            jmax = j;
            max = A[z][j];
           }
*/
 
/*      // по текущей строке
                imax = i;
        for(int j = i;j < n;j++)
         if(A[i][j] > max)
         {
           jmax = j;
           max = A[i][j];
         }
*/                
/*                                // по текущему столбцу
                jmax = i;
        for(int k = i;k < n;k++)
         if(A[k][i] > max)
         {
            imax = k;
            max = A[k][i];
         }
*/                
 
        ct[i] = jmax;
 
        SwapString(A,imax,i);
        SwapColumn(A,n,jmax,i);
        
        double t = B[imax];
            B[imax] = B[i];
            B[i] = t;   
 
        if(A[i][i] != 0)
        {
            double t = 1./A[i][i];
            for(int j = i;j < n;j++)
            {
                A[i][j] *= t;
            }
 
            B[i] *= t;
 
            for(int k = 0;k < i;k++)
            {
                double s = A[k][i];
                for(int l = i;l < n;l++)
                {
                    A[k][l] -= A[i][l]*s;
                }
 
                B[k] -= B[i]*s;
            }
 
            for(int k = i + 1;k < n;k++)
            {
                double t = A[k][i];
                for(int l = i;l < n;l++)
                {
                    A[k][l] -= A[i][l]*t;
                }
 
                B[k] -= B[i]*t;
            }
        }
    }
    
    for(int i = n - 1;i >= 0;i--)
    {
        double t = B[ct[i]];
            B[ct[i]] = B[i];
            B[i] = t;   
    }
 
    delete [] ct; 
} 
 
double      Norm(double ** p,int n)
{
    double max = p[0][0];
    for(int i = 0;i < n;i++)
        for(int j = 0;j < n;j++)
            if(p[i][j] > max) max = p[i][j];
    return max;
}
 
int main()
{
 
}