guasselimination.cpp

1.求方程组的解算法： Guass迭代法、Doolittle分解法 2.求矩阵特征值算法： 幂法、反幂法、LU分解法

#include <stdio.h>

void GuassCommon(double *pA, double *pB, double *pX, int n)
{
	double akk, mik, temp;
	int k, i, j;

	//elimination
	for (k = 0; k < n - 1; k++)
	{
		akk = pA[k * n + k];
		if (akk < 0.00001 && akk > -0.00001)//akk = 0
		{
			printf("\nError: a[k, k] = 0\n");
			break;
		}
		for (i = k + 1 ; i < n; i++)
		{
			mik = pA[i * n + k] / akk;
			for (j = k + 1; j < n; j++)
			{
				temp = pA[i * n + j];
				pA[i * n + j] = temp - mik * pA[k * n + j];
				temp = pB[i];
			}
			pB[i] = temp - mik * pB[k];
			pA[i * n + k] = 0;//
		}
		/*debug
		for (i = 0; i < n; i++)
		{
		for (j = 0; j  < n; j++)
		{
		printf("a%d%d = %f ", i + 1, j + 1, pA[i * n + j]);
		if ((j + 1) % n == 0)
		{
		printf("\n");
		}
		}
		}
		for (i = 0; i < n; i++)
		{
		printf("b%d = %f ", i + 1, pB[i]);
		if ((i + 1) % n == 0)
		{
		printf("\n");
		}
		}*/
	}
	
	//get result
	pX[n - 1] = pB[n - 1] / pA[(n - 1) * n + (n - 1)];
	for (k = n - 2; k >= 0; k--)
	{
		temp = pB[k];
		for (j = k + 1; j  < n; j++)
		{
			temp -= (pA[k * n + j] * pX[j]);
		}
		pX[k] = temp / pA[k * n + k];
	}
}

void DoolittleCommon(double *pA, double *pB, double *pX, int n)
{
	double *pT = new double[n * n];
	double *pY = new double[n];
	double temp;
	int i, j, k, t;
	
	for (j = 0; j < n; j++)
	{
		pT[j] = pA[j];
		if (j != 0)
		{
			pT[j * n] = pA[j * n] / pA[0];
		}
	}
	
	//LU decompose
	for (k = 0; k < n; k++)
	{
		for (j = k; j < n; j++)
		{
			temp = pA[k * n + j];
			for (t = 0; t < k; t++)
			{
				temp -= (pT[k * n + t] * pT[t * n + j]);
			}
			pT[k * n + j] = temp;
		}
		for (i = k + 1; i < n; i++)
		{
			temp = pA[i * n + k];
			for (t = 0; t < k; t++)
			{
				temp -= (pT[i * n + t] * pT[t * n + k]);
			}
			pT[i * n + k] = temp / pT[k * n + k];
		}
	}
	
	//solving equation
	pY[0] = pB[0];
	for (i = 1; i < n; i++)
	{
		temp = pB[i];
		for (t = 0; t < i; t++)
		{
			temp -= (pT[i * n + t] * pY[t]);
		}
		pY[i] = temp;
	}

	pX[n - 1] = pY[n - 1] / pT[(n - 1) * n + n - 1];
	for (i = n -2; i >= 0; i--)
	{
		temp = pY[i];
		for (t = i + 1; t < n; t++)
		{
			temp -= (pT[i * n + t] * pX[t]);
		}
		pX[i] = temp / pT[i * n + i];
	}
	
	delete []pT;
	delete []pY;
}

int max3(int a, int b, int c)
{
	int d;
	d = a > b ? a : b;
	return (d > c ? d : c);
}
int min2(int a, int b)
{
	return (a < b ? a : b);
}
int Doolittle(double *c, double *b, double *xOut, int length, int s, int r)
{
	int i, j, k, t, max, min;
	double temp;
	for (k = 1; k <= length; k++)
	{
		min = min2(k + s, length);
		for (j = k; j <= min; j++)
		{
			temp = 0;
			max = max3(1, k - r, j - s);
			for (t = max; t <= (k - 1); t++)
			{
				temp += ( c[(k - t + s) * length + t - 1] * c[(t - j + s) * length + j - 1] );
			}
			c[(k - j + s) * length + j - 1] = c[(k - j + s) * length + j - 1] - temp;
		}
		min = min2(k + r, length);
		for (i = k + 1; i <= min && k < length; i++)
		{
			temp = 0;
			max = max3(1, i - r, k - s);
			for (t = max; t <= (k - 1); t++)
			{
				temp += ( c[(i - t + s) * length + t - 1] * c[(t - k + s) * length + k - 1] );
			}
			c[(i - k + s) * length + k - 1] = (c[(i - k + s) * length + k - 1] - temp) / c[s * length + k - 1];
		}
	}

	for (i = 0; i < s + r + 1; i ++)
	{
		for (j = 0; j < length; j++)
		{
			printf("%f ", c[i * length + j]);
		}
		printf("\n");
	}
		
	for (i = 2; i <= length; i++)
	{
		temp = 0;
		max = max3(1, 1, i - r);
		for (t = max; t <= (i - 1); t++)
		{
			temp += ( c[(i - t + s) * length + t - 1] * b[t - 1] );
		}
		b[i - 1] = b[i - 1] - temp;
	}

	xOut[length - 1] = b[length - 1] / c[s * length + length - 1];
	for (i = length - 1; i >= 1; i--)
	{
		temp = 0;
		min = min2(i + s, length);
		for (t = i + 1; t <= min; t++)
		{
			temp += c[(i - t + s) * length + t - 1] * xOut[t - 1];
		}
		xOut[i - 1] = (b[i - 1] - temp) / c[s * length + i - 1];
	}

	return 1;
}

int main()
{
	int i, n;
	double X[5];
	n = 5;
/*
	double A1[] = {0.012, 0.01, 0.167,
	1,       0.8334, 5.19,
	3200, 1200, 4.2};
	double B1[] = {0.6781,
	12.1,
	981};*/
	
	double A2[] = {
		3, 2, 1, 0, 0,
		2, 3, 2, 1, 0,
		1, 2, 3, 2, 1,
		0, 1, 2, 3, 2,
		0, 0, 1, 2, 3
	};
	double B2[] = {
		6,
		6,
		6,
		6,
		6
	};
	double C[] = {
		0, 0, 1, 1, 1,
		0, 2, 2, 2, 2,
		3, 3, 3, 3, 3,
		2, 2, 2, 2, 0,
		1, 1, 1, 0, 0
	};
	double B3[] = {
		6,
		6,
		6,
		6,
		6
	};
	double B1[] = {
		6,
		6,
		6,
		6,
		6
	};

	//Guass Common
	GuassCommon(A2, B1, X, n);
	
	for (i = 0; i < n; i++)
	{
		printf("X%d = %f\n", i + 1, X[i]);
	}

	
	//Doolittle
	DoolittleCommon(A2, B3, X, n);
	
	for (i = 0; i < n; i++)
	{
		printf("X%d = %f\n", i + 1, X[i]);
	}

	Doolittle(C, B2, X, 5, 2, 2);
	for (i = 0; i < n; i++)
	{
		printf("X%d = %f\n", i + 1, X[i]);
	}

	return 0;
}