sparse_2.cpp

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

//=================================================================
//                                                                =
//                 大型稀疏线性方程组的求解   A*X=b               =
//                                                                =
//                                                                =
//                                                                =
//=================================================================                       

/*参数解释：
	N	线性方程组的阶数
	M	系数矩阵A中非零元素的个数
	a1	存放系数矩阵A中非零元素的一维数组
	a2	存放系数矩阵A中非零元素所在列的列号的一维数组
	ma	存放系数矩阵A中各行所具有的非零元素的个数的一维数组
	b	存放方程组右端常数向量的数组，输出时用于存放方程组的解
	eps	控制常数，为较小的在正数，当行主元的绝对值小于此值时，
		认为方程组无解，例如，可取1.0e-10
*/

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

#define N 4
#define M 7

double a1[M]={2.0,6.0,3.0,7.0,1.5,10.0,8.0};	//按行的顺序存放系数矩阵A中非零元素的一维数组
double b[N]={13.0,8.5,2.25,26.0};
int a2[M]={0,3,0,1,2,1,3};	//非零元素的列号
int ma[N]={2,2,1,2};      //各行非零元素的个数

int LinEq_Sparse(int n,int m,double *a1,int *a2,int *ma,double *b,double eps);

void main()
{
	if(LinEq_Sparse(N,M,a1,a2,ma,b,1.0e-10)==0)	 //线性方程组有解
	{
		cout<<"The Results : \n"<<endl;
		for(int i=0;i<N;i++)
			cout<<"x"<<"["<<i<<"]"<<"   "<<b[i]<<endl;
	}
	else 
		cout<<"Can't find solutions! "<<endl;
}

int LinEq_Sparse(int n,int m,double *a1,int *a2,int *ma,double *b,double eps)
{
	double c;
	int j,ik,ii,jj,ll,j0;

	double *wk=new double[N];	
	int *nv=new int[N];
	double *v1=new double[N*N];
	int *v2=new int[N*N];

	ik=0;	//非零元素序号
	for(int i=0;i<n;++i)	//n=4
	{
		for( j=0;j<n;++j)
			wk[j]=0.0;	//工作数组
		c=0.0;
		
		cout<<" ma[i] = "<<ma[i]<<endl;
		for(int l=0;l<ma[i];++l)	//ma[i], 系数矩阵第 i 行非零元素个数 
		{
			j=a2[ik];   //非零元素列号赋给j
			wk[j]=a1[ik];	//从a1数组中取出原A中第i行的非零元素, 按原位置存放于工作数组
			
			cout<<" ik= "<<ik<<"  ";			
			cout<<a2[ik]<<"  ";
			cout<<a1[ik]<<endl;

			ik++;
		}
		cout<<"\n";

	
		jj=0;
		if(i!=0)	//第一行除外
		{
			for(l=0;l<i;++l)	//消元过程
			{
				ii=ma[l];	// ii   第 i 行非零元素的个数

				if(fabs(wk[ii])>1.0e-10)
				{
					if(ll==nv[l])
					{
						for(int k=0;k<ll;++k)
						{
							j=v2[jj];
							wk[j]-=wk[ii]*v1[jj];
							++jj;
						}
					}
					b[i]-=wk[ii]*b[l];
					wk[ii]=0.0;
				}
				else
					jj+=nv[l];
			}
		}





		 l=0;
		for( j=0;j<n;++j)	//选行主元
		{
			if(wk[j]!=0.0&&fabs(wk[j])>fabs(c))      
			{
				c=wk[j];
				j0=j;
			}
		}


		if(fabs(c)<eps)
			return(-1);   //方程组无解，返回 -1


		ma[i]=j0;
		for(j=0;j<n;j++)
		{
			if(wk[j]!=0.0&&j!=j0)
			{
				v1[jj]=wk[j]/c;
				v2[jj]=j;
				jj++;
				l++;
			}
		}
		nv[i]=l;
		b[i]/=c;
	}


	for(i=n-1;i>=0;i--)
	{
		//wk[ma[i]]=b[i];
		wk[ma[i]]=b[i];
		if(i==0)
			continue;
		ii=nv[i-1];
		if(ii!=0)
		{
			for(int l=0;l<ii;l++)
			{
				--jj;
				j=v2[jj];
				b[i-1]-=wk[j]*v1[jj];
			}
		}
	}

	for(i=0;i<n;i++)	//保存解
		b[i]=wk[i];

	cout<<"b[i]="<<b[i]<<endl;

	delete []nv;
	delete []v1;
	delete []v2;
	delete []wk;
	return(0);	//方程组有解，返回 0
}