inverse.cpp

1.大型稀疏线性方程组的求解 A*X=b 。 2. 一维数组冒泡法排序算法 4.矩阵求逆 5. 改进的牛顿算法——弦割法

//===================================================
//                                                  =
//                    矩阵求逆                      =
//                                                  =
//                                                  =
//                                                  =
//                                                  =
//===================================================

//	参数说明：
//		N	方阵的阶次
//		a	输入时存放矩阵数据，按行存放，返回时就是该矩阵的逆阵
//		eps	误差控制项，一般取exp(-10)
//
//	返回值说明：
//		返回值就是矩阵的秩。
//		n			矩阵满秩
//		0--(n-1)	矩阵非满秩
//		-1			该矩阵无逆矩阵
		

#include <stdio.h>
#include <math.h>

#define N 4

void dswap(double *x,double *y);
int MatrixInverse(int n,double a[N][N],double eps);

void main()
{
/*	double a[N][N]={{0.68,0.47,0.68,1.21},
					{2.16,0.12,1.97,2.49},
					{0.58,5.16,0.18,0.71},
					{1.96,1.07,1.68,0.71}};*/

	double a[N][N]={{5,7,6,5},
					{7,10,8,7},
					{6,8,10,9},
					{5,7,9,10}};

	MatrixInverse(N,a,1.0e-10);

	for(int i=0;i<N;i++)
	{
		for(int j=0;j<N;j++)
		{
			printf("%10.4f   ",a[i][j]);
		}
		printf("\n");
	}
}

int MatrixInverse(int n,double a[N][N],double eps)
{
	int i,j,k;
	double t,d;
	int *is=new int[N];
	int *js=new int[N];

	if(is==NULL||js==NULL)
	{
		delete []is;
		delete []js;
		return(-1);
	}

	for( k=0;k<n;k++)
	{
		d=0.0;
		for(int i=k;i<n;i++)
		{
			for( j=k;j<n;j++)
			{
				t=fabs(a[i][j]);
				if(t>d)
				{
					d=t;is[k]=i;js[k]=j;
				}
			}
		}

		if(d<=eps)
		{
			delete []is;
			delete []js;
			return(k);
		}

		for( j=0;j<n;j++)
		{
			i=is[k];
			dswap(&a[k][j],&a[i][j]);
		}

		for(i=0;i<n;i++)
		{
			j=js[k];
			dswap(&a[i][k],&a[i][j]);
		}

		a[k][k]=1.0/a[k][k];

		for(j=0;j<n;j++)
		{
			if(j!=k)	a[k][j]*=a[k][k];
		}

		for(i=0;i<n;i++)
		{
			if(i!=k)
			{
				for(int j=0;j<n;j++)
				{
					if(j!=k)	a[i][j]-=a[i][k]*a[k][j];
				}
			}
		}

		for(i=0;i<n;i++)
		{
			if(i!=k)	a[i][k]*=-1.0*a[k][k];
		}
	}

	for( k=n-1;k>=0;k--)
	{
		for(int j=0;j<n;j++)
		{
			i=js[k];
			dswap(&a[k][j],&a[i][j]);
		}

		for(int i=0;i<n;i++)
		{
			j=is[k];
			dswap(&a[i][k],&a[i][j]);
		}
	}

	delete []is;delete []js;
	return(n);
}

void dswap(double *x,double *y)
{
	double temp;
	temp = *x;
	*x=*y;
	*y=temp;
}