gauss.cpp

采用高斯消元法实现的n阶线性方程组求解程序；

#include "stdafx.h"
#include "memory.h"

#define _FABS(x) (x>0?x:(-x))

bool guass(double *g, double *b, double *x, int n) 
{ 
	int i,j,k,p; 
	int *flag; 
	double l; 
	flag = new int[n];
	if( !flag )return false;
	memset(flag,0,sizeof(int)*n);
	
	//按列循环
	for (i=0;i<n;i++)
	{
		//求出i列绝对值最大数
		p=-1; 
		for (j=0;j<n;j++) 
		{
			if(!flag[j] && _FABS((g[j*n+i]))>1e-8)
			{
				if(p==-1) p=j; 
				else if( _FABS((g[p*n+i]))<=_FABS((g[j*n+i])) ) p=j; 
			}
		}   
		
		if (p==-1) return false; 
		if( _FABS((g[p*n+i]))<1e-8 )return false;

		//flag[p]为零，表示p行还没有主元，为1+i，表示此行的主元是在i列上
		flag[p]=1+i; 
		
		//逐行消去第i列系数元
		for (j=0;j<n;j++) 
		{
			if (p!=j)
			{ 
				l=g[j*n+i]/g[p*n+i]; 
				for (k=i;k<n;k++)
				{
					g[j*n+k]-=(g[p*n+k]*l); 
				}
				b[j]-=(b[p]*l); 
			} 
			
		}
	}

	//算出未知数值
	for (j=0;j<n;j++)
	{
		//flag[j]记录了j行的主元所在列为flag[j]－1；
		x[flag[j]-1]=b[j]/g[j*n+flag[j]-1];
	}
	
	delete[] flag;
	return true;
}

//矩阵相乘
void matrix_multiply(double *m1, double *m2, int n, double *r)
{
	int i,j,k;
	double v;
	for( i=0; i<n; i++ )
	{
		for( j=0; j<n; j++ )
		{
			v = 0.0;
			for( k=0; k<n; k++)
			{
				v += (m1[i*n+k]*m2[k*n+j]);
			}

			r[i*n+j] = v;
		}
	}

	return;
}

//置矩阵为单位矩阵
void matrix_toIdentity(double *m, int n)
{
	int i,j;
	for( i=0; i<n; i++ )
	{
		for( j=0; j<n; j++ )
		{
			if( i==j )m[i*n+j] = 1.0;
			else m[i*n+j] = 0.0;
		}
	}

	return;
}

//判断矩阵是否为单位矩阵
bool matrix_isIdentity(double *m, int n)
{
	int i,j;
	for( i=0; i<n; i++ )
	{
		for( j=0; j<n; j++ )
		{
			if( i==j && _FABS((m[i*n+j]-1))>1e-8 )return false;
			else if( i!=j && _FABS((m[i*n+j]))>1e-8 )return false;
		}
	}

	return true;
}

//求逆矩阵（采用初等变换，高斯列主消元法）
bool matrix_reverse(double *m, int n, double *r)
{
	int i,j,k,p; 
	int *flag; 
	double l; 
	double *tm;

	flag = new int[n];
	tm = new double[n*n];
	if( !flag || !tm )
	{
		if( tm )delete[] tm;
		if( flag )delete[] flag;
		return false;
	}

	memset(flag,0,sizeof(int)*n);
	memcpy(tm,m,sizeof(double)*n*n);

	//将r置为单位方阵
	matrix_toIdentity(r,4);
	
	//按列循环
	for(i=0;i<n;i++)
	{
		//求出i列绝对值最大数
		p=-1; 
		for(j=0;j<n;j++) 
		{
			if(!flag[j] && _FABS((m[j*n+i]))>1e-8)
			{
				if(p==-1) p=j; 
				else if( _FABS((m[p*n+i]))<_FABS((m[j*n+i])) ) p=j; 
			}
		}   
		
		if( p==-1 ) return false;
		if( _FABS((m[p*n+i]))<1e-8 )return false;

		//flag[p]为零，表示p行还没有主元，为1+i，表示此行的主元是在i列上
		flag[p]=1+i; 
		
		//逐行消去第i列系数元(不含拥有列主元的那一行)
		for(j=0;j<n;j++) 
		{
			if(p!=j)
			{ 
				l=m[j*n+i]/m[p*n+i]; 
				for(k=i;k<n;k++)
				{
					m[j*n+k]-=(m[p*n+k]*l);
				}

				for(k=0;k<n;k++)
				{				
					r[j*n+k]-=(r[p*n+k]*l); 
				}
			}
		}
		
		//将拥有列主元的那一行的系数归一化
		l=1.0/m[p*n+i]; 
		for(k=i;k<n;k++)
		{
			m[p*n+k]=(m[p*n+k]*l);
		}

		for(k=0;k<n;k++)
		{
			r[p*n+k]=(r[p*n+k]*l); 
		}		
	}

	//置换各行，使 m 矩阵对角化(为单位矩阵)
	for(i=0; i<n; i++)
	{	
		memcpy( m+(flag[i]-1)*n, r+i*n, sizeof(double)*n);
	}

	memcpy(r,m,sizeof(double)*n*n);
	memcpy(m,tm,sizeof(double)*n*n);

	delete[] tm;
	delete[] flag;

	return true;
}

//求逆矩阵（采用初等变换，高斯列主消元法）
bool matrix_reverse2(double *m, int n, double *r)
{
	int i,j,k,p; 
	int *flag; 
	double l; 
	double t;

	flag = new int[n];
	if( !flag )return false;

	memset(flag,0,sizeof(int)*n);

	//将r置为单位方阵
	matrix_toIdentity(r,4);
	
	//按列循环
	for(i=0;i<n;i++)
	{
		//求出i列绝对值最大数
		p=-1; 
		for(j=0;j<n;j++) 
		{
			if(!flag[j] && _FABS((m[j*n+i]))>1e-8)
			{
				if(p==-1) p=j; 
				else if( _FABS((m[p*n+i]))<_FABS((m[j*n+i])) ) p=j; 
			}
		}   
		
		if( p==-1 ) return false;
		if( _FABS((m[p*n+i]))<1e-8 )return false;

		//flag[p]为零，表示p行还没有主元，为1+i，表示此行的主元是在i列上
		flag[p]=1+i; 
		
		//逐行消去第i列系数元(不含拥有列主元的那一行)
		for(j=0;j<n;j++) 
		{
			if(p!=j)
			{ 
				l=m[j*n+i]/m[p*n+i]; 
				for(k=i;k<n;k++)
				{
					m[j*n+k]-=(m[p*n+k]*l);
				}

				for(k=0;k<n;k++)
				{				
					r[j*n+k]-=(r[p*n+k]*l); 
				}
			}
		}
		
		//将拥有列主元的那一行的系数归一化
		l=1.0/m[p*n+i]; 
		for(k=i;k<n;k++)
		{
			m[p*n+k]=(m[p*n+k]*l);
		}

		for(k=0;k<n;k++)
		{
			r[p*n+k]=(r[p*n+k]*l); 
		}		
	}

	//置换各行，使 m 矩阵对角化(为单位矩阵)
	for(i=0; i<n; i++)
	{
		// 找到主元为 i 列的行（记为j）
		for(j=i; j<n; j++)
		{
			if( flag[j]-1==i )break;
		}
		if( j>=n )return false;

		// 与当前行置换
		if( j==i )continue;
		
		for(k=0; k<n; k++)
		{
			t = r[j*n+k]; r[j*n+k] = r[i*n+k]; r[i*n+k] = t;
		}

		flag[j] = flag[i];
	}
	
	delete[] flag;
	return true;
}