C语言实现Crout分解

PS：3、4公式要同时进行，第4公式的j要从k开始。

#include<stdio.h>
#define N 4
#define col 4//列数
#define row 4//行数
void print(double a[row][col]) {

	for (int i = 0; i < row; i++)//找到最大主元，记录行号和列号
	{
		for (int j = 0; j < col; j++)
		{
			printf("%f  ", a[i][j]);
		}
		printf("\n");
	}
	printf("\n");
}

void printXY(double XY[N]) {

	for (int i = 0; i < N; i++)
	{
		printf("%f  ", XY[i]);
		
	}
	printf("\n\n");
}

void step1(double a[row][col], double L[row][col], double  U[row][col])
{
	for (int i = 1; i <= N; i++)
	{
		L[i - 1][0] = a[i - 1][0];
	}
}

void step2(double a[row][col], double L[row][col], double  U[row][col])
{
	for (int j = 1; j <= N; j++)
	{
		U[0][j-1]=a[0][j-1]/L[0][0];
	}
}

void step3_4(double a[row][col], double L[row][col], double  U[row][col])
{
	for (int k = 2; k <= N; k++)
	{
		for (int i = k; i <= N; i++)
		{
			double sum = 0;
			for (int r = 1; r <= k - 1; r++)
			{
				sum = sum + L[i - 1][r - 1] * U[r - 1][k - 1];
			}
			L[i - 1][k - 1] = a[i - 1][k - 1] - sum;
		}

		for (int j = k ; j <= N; j++)
		{
			double sum = 0;
			for (int r = 1; r <= k - 1; r++)
			{
				sum = sum + L[k - 1][r - 1] * U[r - 1][j - 1];
			}
			U[k - 1][j - 1] = (a[k - 1][j - 1] - sum) / L[k - 1][k - 1];
		}
	}
}

void computeX( double  U[row][col], double X[N], double Y[N])
{
	X[N - 1] = Y[N-1];
	for (int k = N - 1; k >= 1; k--)
	{
		double sum = 0;
		for (int r = k + 1; r <= N; r++)
		{
			sum += U[k-1][r-1] * X[r-1];
		}
		X[k-1] = Y[k-1] - sum;
	}
}

void computeY( double L[row][col],double b[N],double Y[N])
{
	Y[0] = b[0] / L[0][0];
	for (int k = 2; k <=N; k++)
	{
		double sum = 0;
		for (int r = 1; r <= k - 1; r++)
		{
			sum += L[k-1][r-1] * Y[r-1];
		}
		Y[k-1] =( b[k-1]-sum)/L[k-1][k-1];
	}
}

int main()
{
	double a[N][N] = {	{ 2,4,2,6},
						{ 4,9,6,15 },
						{ 2,6,9,18 },
						{ 6,15,18,40 } }, 
		b[N] = {9,23,22,47},
		L[N][N] = { {0} }, U[N][N] = { { 0 } }, Y[N] = { 0 }, X[N] = { 0 };
	


	step1(a, L, U);
	step2(a, L, U);
	step3_4(a, L, U);
	computeY(L,b,Y);
	computeX(U,X,Y);

	printf("y=\n");
	printXY(Y);
	printf("X=\n");
	printXY(X);
	printf("L=\n");
	print(L);
	printf("U=\n");
	print(U);
	
	getchar();
}