doolittle.cpp

这是数值分析中doolittle算法的实现程序

/************************write by bear********************
***          选主元的doolittle分解法求Ax=b             ***
**********************************************************/

#include "stdio.h"
#include "math.h"
#define n 4

////////////打印函数/////////////
void print(double A[n][n])
{
	printf("A is:\n");
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
		{	
			
			printf("%8.3f,",A[i][j]);
			if(j==n-1) printf("\n");
		}
}
void print1(char ar[2],double *arr)
{
	printf("%s is:\n",ar);
	for(int j=0;j<n;j++)				
	{	
		if(j==0) printf("[");
		printf("%8.3f",*(arr+j));
		if(j!=n-1) printf(",");
		else printf("\t]");
	}
	printf("\n");
}
void print2(char ar[2],int *arr)
{
	printf("%s is:\n",ar);
	for(int j=0;j<n;j++)				
			printf("%d,",*(arr+j));
	printf("\n");
}
///////求数组的最大值////////////
int Max(int k1,double *array)
{
	double max=0;
	int maxnum=0;
	for(int kk=k1-1;kk<n;kk++)
		if(max<fabs(*(array+kk))) {max=fabs(*(array+kk));maxnum=kk;}
	printf("maxnum is %d,max is %5.3f\n",maxnum,max);
	return maxnum;
}
//////////////doolittle分解算法//////////////
void Doolittle(double A[n][n],double b[n])
{
	int M[n];//用于记录第k个主元素所在的行号
	int k;
	int i,t;
	double S[n];//中间量
	double y[n];
	double x[n];
	double temp=0;
	int ik;
	for(k=1;k<=n;k++)
	{
		//计算中间量si
		for(i=k;i<=n;i++)
		{
			for(t=1,temp=0;t<=k-1;t++)
			{
				temp=temp+A[i-1][t-1]*A[t-1][k-1];
			}
			S[i-1]=A[i-1][k-1]-temp;
			A[i-1][k-1]=S[i-1];
			//print();
		}
		ik=Max(k,S);//(2)选行号ik
		M[k-1]=ik;
		printf("ik=%d,k-1=%d\n",ik,k-1);
		if(ik!=k-1)//(3)
		{
			for(int j=1;j<=n;j++)
			{temp=A[k-1][j-1];A[k-1][j-1]=A[ik][j-1];A[ik][j-1]=temp;}
			//print();	
		}
		//////(4)计算ukk,ukj,lik//////
		
		for(int j=k+1;j<=n;j++)
		{
			for(t=1,temp=0;t<=k-1;t++)
			{
				temp=temp+A[k-1][t-1]*A[t-1][j-1];
			}
			A[k-1][j-1]=A[k-1][j-1]-temp;//计算ukj
			A[j-1][k-1]=A[j-1][k-1]/A[k-1][k-1];//计算lik
			//print();
		}
		
	}
	print(A);
	print2("M",M);
	//求Qb
	for(k=1;k<=n;k++)
	{
		int tt;
		tt=M[k-1];
		temp=b[k-1];b[k-1]=b[tt];b[tt]=temp;
	}
	print1("b",b);
	//求解Ly=Qb和Ux=y
	y[0]=b[0];//y1
	for(i=2;i<=n;i++)
	{
		for(t=1,temp=0;t<=i-1;t++)
		{
			temp=temp+A[i-1][t-1]*y[t-1];
		}
		y[i-1]=b[i-1]-temp;//yi
	}
	print1("y",y);
	x[n-1]=y[n-1]/A[n-1][n-1];//xn
	for(i=n-1;i>=1;i--)
	{
		for(t=i+1,temp=0;t<=n;t++)
		{
			temp=temp+A[i-1][t-1]*x[t-1];
		}
		x[i-1]=(y[i-1]-temp)/A[i-1][i-1];//xi
	}
	print1("x",x);
}


int main()
{	
	double A[n][n]={1,8,2,3,-6,-3,8,1,2,4,4,2,10,5,-5,6};
	//double A[n][n]={8.3,2.1,-1.2,0.5,0.8,10.2,3.5,-1.8,1.2,0.2,-4,-0.5,-0.2,0.3,0.4,-2};
	//double b[n]={-3.02,4.79,-6.72,8.89};
	double b[n]={12,40,-50,80};
	Doolittle(A,b);
	return 0;
}