gauss.cpp

利用高斯全选主元消去解线性方程组

// Gauss.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <math.h>
#include <iostream>
#include <iomanip>

using namespace std;

void Gauss(double **pCoff, double *pConst,int n);

void main()
{
	int i,j;
    double **a,*b;
	a = new double *[3];
	for (i=0;i<3;i++)
	{
		a[i] = new double [3];
	}
	//double a[3][3]={{0.012,0.01,0.167},{1.0,0.8334,5.91},{3200,1200,4.2}};
	cout << "请输入系数矩阵：" <<endl;
	for (i=0;i<3;i++)
	{
		for (j=0;j<3;j++)
			cin >> a[i][j];
	}
	cout << "您所输入的系数矩阵为：" << endl;
	for (i=0;i<3;i++)
	{
		for (j=0;j<3;j++)
			cout << setw(8) << a[i][j];
			cout << endl;
	}
	//double b[3]={0.6781,12.1,981};

    b = new double [3];
    cout << "请输入常数矩阵：" << endl;
	for (i=0;i<3;i++)
	{
		cin >> b[i];
	}
	cout << "您所输入的常数矩阵为：" << endl;
	for (i=0;i<3;i++)
		cout << b[i] << endl;

	Gauss(a,b,3);
	for (i=0;i<3;i++)
		cout << b[i] <<endl;

}

void Gauss(double **pCoff,double *pConst,int n)
{
	double t,d;
	int i,j,k,row,flag,*pCol;

	// 临时缓冲区，存放列数
	pCol = new int [n];

	//消元
	flag = 1;

	for (k=0;k<n-1;k++)
	{
		d = 0.0;
		for (i=k;i<n;i++)
		{
			for(j=k;j<n;j++)
			{
				t = fabs(pCoff[i][j]);
				if (t>d)
				{
					d = t;
					pCol[k] = j;
					row = i;
				}
			}
		}

		if (d==0.0)
			flag = 0;

		//列交换
		else 
		{
			if (pCol[k]!=k)
			{
				for (i=0;i<n;i++)
				{
					t = pCoff[i][k];
					pCoff[i][k] = pCoff[i][pCol[k]];
					pCoff[i][pCol[k]] = t;
				}
			}
			//行交换
			if (row != k)
			{
				for (j=k;j<n;j++)
				{
					t = pCoff[k][j];
					pCoff[k][j] = pCoff[row][j];
					pCoff[row][j] = t;
				}
				
			   t=pConst[k]; 
			   pConst[k]=pConst[row]; 
			   pConst[row]=t;
			}
		}
	

        //没有解
		if(flag == 0)
		{
			delete [] pCol;
			cout << "没有解！?" <<endl;
		}

		d = pCoff[k][k];
		for (j=k;j<n;j++)
		{
			pCoff[k][j] = pCoff[k][j]/d;
		}
		
		pConst[k] = pConst[k]/d;

		for (i=k+1;i<n;i++)
		{
			for(j=k+1;j<n;j++)
			{
				pCoff[i][j] = pCoff[i][j]-pCoff[i][k]*pCoff[k][j];
			}
			pConst[i] = pConst[i]-pCoff[i][k]*pConst[k];
		}
	}

		//求解失败
		d = pCoff[n-1][n-1];
		if(d==0)
		{
			delete [] pCol;
			cout << "没有解！" <<endl;
		}

		//代入求解
		pConst[n-1] = pConst[n-1]/d;
		for (i=n-2;i>=0;i--)
		{
			t = 0.0;
			for (j=i+1;j<n;j++)
			{
				t = t+pCoff[i][j]*pConst[j];
			}
			pConst[i] =  pConst[i]-t;
		}

		//调整解的位置
    	pCol[n-1]=n-1;
        for (k=n-1;k>=0;k--)
		{
    		if (pCol[k]!=k)
			{ 
    			t=pConst[k]; 
    			pConst[k]=pConst[pCol[k]]; 
    			pConst[pCol[k]]=t;
			}
		}


		//显示结果
		/*cout << "解为：" << endl;
		for (i=0;i<n;i++)
		{
			cout << pConst[i];
			cout << endl;
		} */

		//回收内存
		delete [] pCol;

}