dooli.cpp

doolittle 方法解决矩阵算法问题

/************************write by bear********************
***          doolittle分解法求Ax=b                      ***
**********************************************************/

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

////////////打印函数/////////////
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");
}
//////////////doolittle分解算法//////////////
void Doolittle(double A[n][n],double b[n])
{
	int k;
	int i,t;
	double y[n];
	double x[n];
	double temp=0;
	for(k=1;k<=n;k++)
	{
		//////(4)计算ukj,lik//////		
		for(int j=k;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
		}
		for(i=k+1;i<=n;i++)
		{
			for(t=1,temp=0;t<=k-1;t++)
			{
				temp=temp+A[i-1][t-1]*A[t-1][k-1];
			}
			A[i-1][k-1]=(A[i-1][k-1]-temp)/A[k-1][k-1];//计算lik
		}
		
	}
	print(A);
	//求解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]={8.1,2.3,-1.5,0.5,-6.23,0.87,2.5,1.5,10.2};
	//double A[n][n]={3.2,-1.4,1,6.4,2.8,3,5.5,5.2,-4};
	//double b[n]={2.5,1.8,7.2};
	double b[n]={6.1,2.3,1.8};
	Doolittle(A,b);
	return 0;
}