Gauss列主元消去法（封裝函式）

//高斯列主元消去法（封裝函式一）
#include<bits/stdc++.h>
#include "windows.h"
using namespace std;
int SolverEqGuass(double **A,double *b,int n,double eps,double *x)
{
	int m,i,j,k;
	double **a =new double *[n];	
	double *btemp= new double[n];
	double d,temp;
	for(i=0;i<n;i++)
	{
		a[i]=new double[n+1];	//系數矩陣+常數項矩陣
	}
	m=n;
	//賦值
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			a[i][j]=A[i][j];
		}
		a[i][n]=b[i];
	} 
	if(m<=0||n<=0)
	{
		cout<<"Input Mistake"<<endl;
	}
	for(k=0;k<n-1;k++)	//找列主元最大值
	{
		double max=0;
		int row=k,num=0;
		for(i=k;i<n;i++)
		{
			if(fabs(a[i][k])>max)
			{
				max=fabs(a[i][k]);
				row=i;	//記錄換行行標
			}
		}
		if(a[row][i]==0)
		{
			cout<<"無法計算"<<endl; 
		}
		if(k!=row)//換行
		{
		 	for(i=0;i<m+1;i++)//該行全部元素交換
		 	{
		 		temp=a[k][i];
		 		a[k][i]=a[row][i];
		 		a[row][i]=temp;
			 }
		}
		for(i=k+1;i<m;i++)
		{
			d=a[i][k]/a[k][k];	//定義解法
			for(j=0;j<n+1;j++)
			{
				a[i][j]=a[i][j]-d*a[k][j];
			}
		}
	} 
	for(i=0;i<n;i++)
	{
		btemp[i]=0.0;
	}
	for(i=m-1;i>=0;i--)	//求解向量
	{
		d=0;
		for(k=0;k<n;k++)
		{
			d=d+btemp[k]*a[i][k];
		}
		btemp[i]=(a[i][n]-d)/a[i][i]; 
	} 
	for(i=0;i<m;i++)
	{
		x[i]=btemp[i];
	}
	for(i=0;i<n;i++)
	{
		delete[]a[i];
	}
	delete[]a;
	delete[]btemp;
}
int main()
{
	LARGE_INTEGER lFrequency, lEndCount, lBeginCount;//LARGE_INTEGER為資料型別
	QueryPerformanceFrequency(&lFrequency);//獲取CPU的時鐘頻率
	QueryPerformanceCounter(&lBeginCount);//獲取CPU計數器數字，放在計時開頭
	double CompuTime;
	double **A,*b,*x;
	int kk=4;	
	A=new double *[kk];
	b=new double[kk];
	x=new double[kk];
	for(int i=0;i<kk;i++)
	{
		A[i]=new double[kk];
	}
	int Loop_num=1e5;
	int n=4;	//根據題目選擇賦值或輸入
	A[0][0]=1;A[0][1]=-1;A[0][2]=2;A[0][3]=-1;
	A[1][0]=2;A[1][1]=-2;A[1][2]=3;A[1][3]=-3;
	A[2][0]=1;A[2][1]=1;A[2][2]=1;A[2][3]=0;
	A[3][0]=1;A[3][1]=-1;A[3][2]=4;A[3][3]=3;
	b[0]=-8;b[1]=-20;b[2]=-2;b[3]=4;
	SolverEqGuass(A,b,n,1.0e-10,x);
	for(int i=0;i<n;i++)
	{
		cout<<"x"<<i+1<<"=";
		printf("%.6f\n",x[i]);	//輸出六位小數
	}
	QueryPerformanceCounter(&lEndCount);//獲得CPU計數器數字，放在計時結尾
	CompuTime = (double)(lEndCount.QuadPart - lBeginCount.QuadPart) / (double)lFrequency.QuadPart;//得到運行時間，單位微秒
	cout<<"Gauss法計算時間為"<<CompuTime; 
}

//封裝函式二
void SolverEqGauss(double** A, double* b, int n, double* x)
{
    double aef = 0, bt = 0, gm = 0, delta = 0, omiga = 0;
    int lamda = 0;
    bool bo = false;
    for (int I = 1;I < n;I++)
    {
        delta = abs(A[I][I]);
        lamda = I;
        for (int k = I+1;k < n + 1;k++)
        {
            if (abs(A[k][I]) > delta)
            {
                delta = abs(A[k][I]);
                lamda = k;
            }
        }
        if (delta < pow(0.1, 8))
            bo = true;
        if (bo)
        {
            cout << "此方程無解";
            break;
        }
        else
        {
            for (int l = I;l < n + 1;l++)
            {
                omiga = A[I][l];
                A[I][l] = A[lamda][l];
                A[lamda][l] = omiga;
            }
            omiga = b[I];
            b[I] = b[lamda];
            b[lamda] = omiga;
        }
        for (int i = I+1;i < n + 1;i++)
        {
            aef = A[i][I] / A[I][I];

            for (int j = I;j < n + 1;j++)
            {
                A[i][j] = A[i][j] - A[I][j]*aef;                
            }
            b[i] = b[i] - b[I]*aef;
        }
    }
    
    for (int u = n;u > 0;u--)
    {
        for (int v = n;v > u ;v--)
        {
            b[u] -= A[u][v] * x[v];
        }
        x[u] = b[u] / A[u][u];
    }
}