doolittle.cpp

Doolittle的编程实现

# include <stdio.h>
# include <math.h>
# define N 10
# define M 10
void main ()
{
	int n,m,i,j,k,r,I;
	float c,w,p=0;
	float x[N][M],s[N],a[N][N+M];
	printf("输入系数矩阵的阶n:");
	scanf("%d",&n);
	printf("输入右端项的列数m:");
	scanf("%d",&m);
	for(i=1;i<n+1;i++)
	{
		for(j=1;j<n+m+1;j++)
		{
			printf("a[%d][%d]=",i,j,a[i][j]);
			scanf("%f",&a[i][j]);
		}
	}
	//选主元
	for(r=1;r<n;r++)
	{
		for(i=r;i<=n;i++)
		{
			s[i]=a[i][r];
		    for(k=1;k<r;k++)
			s[i]=s[i]-a[i][k]*a[k][r];
		    a[i][r]=s[i];
		}
     for(i=r;i<=n;i++)
	{
	   	if(fabs(s[i])>p)
		{
		   	p=s[i];
	    	I=i;
		}
	}
	 //换行
	if(I!=r)	
	{
		for(j=1;j<=n+m;j++)
		{
		  	w=a[r][j];
		    a[r][j]=a[I][j];
		  	a[I][j]=w;
		}
	}
	for(i=r+1;i<=n;i++)
	    a[i][r]=a[i][r]/a[r][r];
		for(j=r+1;j<=n+m;j++)
		{
		   	c=a[r][j];
		   	for(k=1;k<r;k++)
		     	c=c-a[r][k]*a[k][j];
		   	a[r][j]=c;
		}
	}
    for(j=n;j<=n+m;j++)
	{
		c=a[n][j];
	    for(k=1;k<n;k++)
			c=c-a[n][k]*a[k][j];
    	a[n][j]=c;
	}
	//分解
    for(j=1;j<=m;j++)
	{
        a[n][n+j]=a[n][n+j]/a[n][n];
       	for(i=n-1;i>=1;i--)
		{
		   	for(k=i+1;k<=n;k++)
			{
		       a[i][n+j]=a[i][n+j]-a[i][k]*a[k][n+j];
			}
       	   	a[i][n+j]=a[i][n+j]/a[i][i];
		}
	}
	//回代
	printf("输出解:\n");
	for(j=1;j<=m;j++)
	{
		for(i=1;i<=n;i++)
		{
	   	   printf("x[%d][%d]=%f\n",i,j,a[i][n+j]);
		}
	}
}