Guass消元法求解线性方程组

/*高斯消元法求解线性方程组
    高斯消元的消元计算：
               (k)  (k)
        Mik = Aik  / Akk       (i = k+1,k+2,,.....,n)

         (k+1)  (k)          (k)    
        Aij  = Aij  - MikAkj    (i,j= k+1,k+2,.....,n)
    
         (k+1)  (k)          (k)    
        Bi   = Bi   - MikBk        (i = k+1,k+2,.....,n)

    回代求解：
        
               (n)      (n)
        Xn =  Bn   / Ann

               (i)         (i)
        Xi = (Bi    - )/Aii
        (i = n-1,n-2,....,1)
    

    使用高斯消元法求解线性方程组比用Cramer求解的计算量要小的多。

*/

 

/*函数名称：gauss_elimination_calculate
    函数参数：int (*p)[3]；线性方程组的系数行列式
              int* B；线性方程组右边的常数向量
              int size；线性方程组的阶数
    函数返回值：函数没有返回值，输出方程组的求解结果
*/
void gauss_elimination_calculate(double (*a)[3], double* ipvt, int size)
{
    double *X = new double[size];

    //消元，是系数矩阵成为上三角矩阵
    double mi = 1.0;
    int k = 0;
    for(int i = 1; i < size; i++)
    {
        for(int j = 0; j < i; j++)
        {
            mi = a[i][j]/a[j][j];
            k = j;
            while(k < size)
            {
                a[i][k] = a[i][k] - mi * a[j][k];
                k ++;
            }
            ipvt[i] = ipvt[i] - ipvt[j] * mi;
        }
    }


    //回代，求解线性方程组的结果
    X[size - 1] = ipvt[size - 1]/a[size - 1][size - 1];
    for(int i = size - 2; i >= 0; i--)
    {
        double temp = ipvt[i] - X[i + 1] * a[i][i + 1];
        for(int j = i + 2; j < size; j++)
            temp = temp - X[j] * a[i][j];
        X[i] = temp/a[i][i];
    }

    for(int i = 0; i < size; i++)
        std::cout << X[i] << " ";
    std::cout << '\n';

    delete[] X;
}