decomposition.cpp

选主元的矩阵杜利特尔分解

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

void main()
{
	double		a[4][4]={{0.8,10.2,3.5,-1.8},{8.3,2.1,-1.2,0.5},{1.2,0.2,-4.,-0.5},{-0.2,0.3,0.4,-2.}};
	double		b[4]={4.79,-3.02,-6.72,8.89};
	double		s[4];
	double		y[4];
	double		x[4];
	double		temp[4];
	int			m[4];
	double		modulemax;
	double		exchange;
	int			k,i,j,t;

	for(k=0;k<=3;k++)
	{
		for(i=k;i<=3;i++)
		{
			s[i]=a[i][k];
			for(t=0;t<=k-1;t++)
			{
				s[i]-=a[i][t]*a[t][k];
			}
		}
		modulemax=fabs(s[k]);
		m[k]=k;
		for(i=k+1;i<=3;i++)
		{
			if(modulemax<s[i])
			{
				m[k]=i;
				modulemax=s[i];
			}
		}
		if(m[k]!=k)
		{
			exchange=s[k];
			s[k]=s[m[k]];
			s[m[k]]=exchange;
			for(t=0;t<=3;t++)
			{
				temp[t]=a[k][t];
			}
			for(t=0;t<=3;t++)
			{
				a[k][t]=a[m[k]][t];
			}
			for(t=0;t<=3;t++)
			{
				a[m[k]][t]=temp[t];
			}
		}
		a[k][k]=s[k];
		for(j=k+1;j<=3;j++)
		{
			for(t=0;t<=k-1;t++)
			{
				a[k][j]-=a[k][t]*a[t][j];
			}
		}
		for(i=k+1;i<=3;i++)
		{
			a[i][k]=s[i]/a[k][k];
		}
	}
	for(k=0;k<=2;k++)
	{
		if(m[k]>k)
		{
			exchange=b[k];
			b[k]=b[m[k]];
			b[m[k]]=exchange;
		}
	}
	y[0]=b[0];
	for(i=1;i<=3;i++)
	{
		y[i]=b[i];
		for(t=0;t<=i-1;t++)
		{
			y[i]-=a[i][t]*y[t];
		}
	}
	x[3]=y[3]/a[3][3];
	for(i=2;i>=0;i--)
	{
		x[i]=y[i];
		for(t=i+1;t<=3;t++)
		{
			x[i]-=a[i][t]*x[t];
		}
		x[i]/=a[i][i];
	}
}