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

#define MAXLOOPS	2000	//最大迭代次數
#define AU 			2		//上半帶寬
#define AD 			2		//下半帶寬
#define AM 			5		//壓縮矩陣高度
#define AN 			501		//壓縮矩陣寬度
#define ER 			1e-12	//精度水平
#define bb 			0.16
#define cc 			-0.064
#define GetSgn(a)	((a) > 0 ? 1 : ((a) < 0 ? -1 : 0))									//符號函數
#define max3(a, b, c)	(((a) > (b) ? (a) : (b)) > (c) ? ((a) > (b) ? (a) : (b)) : (c))	//求三個數的最大值
#define min2(a, b)	((a) < (b) ? (a) : (b))												//求兩個數的最小值

void main()
{
	int PowERMethod(double *pc, int length, int endLoops, double *resultOut, int type);
	int MaxNumInVector(double *pIn, int length);
	void AMultiY(double *aii, double *inY, double *outU, int length);
	void UnitaryVector(double *pVectorIn, double *pVectorOut, int length, double numbER);
	void Doolittle(double *c, double *b, double *xOut, int length, int s, int r);
	double LU(double *z, int length, int s, int r, int flag);
	
	int true=1;
	int false=0;
	int i, j, loops;
	double c[AN*AM];	//用c壓縮存儲帶狀矩陣A
	double r1, rs, r501, d, mk, cond, det;
	double m[39];
	double tempC[AN*AM];
	
	//把帶狀矩陣A壓縮存儲為C
	for (i = 0; i < AN; i++)
	{
		c[0 * AN + i] = cc;
		c[1 * AN + i] = bb;
		c[2 * AN + i] = ( 1.64 - 0.024 * ( i + 1 ) ) * sin( 0.2 * ( i + 1 ) ) - 0.64 * exp( 0.1 / ( i + 1 ) );
		c[3 * AN + i] = bb;
		c[4 * AN + i] = cc;
	}
	c[ 0 ] = 0;
	c[ 1 ] = 0;
	c[ 1 * AN + 0 ] = 0;
	c[ 3 * AN + AN - 1 ] = 0;
	c[ 4 * AN + AN - 1 ] = 0;
	c[ 4 * AN + AN - 2 ] = 0;
	
	//冪法求最小特征值r1
	for (i = 0; i < AN * AM; i++)
	{
		tempC[i] = c[i];
	}
	loops = PowERMethod(tempC, AN, MAXLOOPS, &r1, true);
	printf("最小特征值r1\t=%1.11e 迭代次數:%d\n", r1, loops);
	
	//冪法求最大特征值，原點平移r501
	for (i = 0; i < AN * AM; i++)
	{
		tempC[i] = c[i];
		if (i >= AU * AN && i < (AU + 1) * AN)
		{
			tempC[i] -= r1;
		}
	}
	loops = PowERMethod(tempC, AN, MAXLOOPS, &r501, true);
	r501 += r1;
	printf("最大特征值r501\t= %1.11e 迭代次數: %d\n", r501, loops);
	
	//反冪法求模最小特征值rs
	for (i = 0; i < AN * AM; i++)
	{
		tempC[i] = c[i];
	}
	LU(tempC, AN, AU, AD, false);
	loops = PowERMethod(tempC, AN, MAXLOOPS, &rs, false);
	printf("模最小特征值rs\t=%1.11e 迭代次數: %d\n", rs, loops);
		
	//用原點平移的反冪法求rik
	d = (r501 - r1) / 40;
	for (j = 1; j < 40; j++)
	{
		mk = r1 + j * d;
		for (i = 0; i < AN * AM; i++)
		{
			tempC[i] = c[i];
			if (i >= AU * AN && i < (AU + 1) * AN)
			{
				tempC[i] -= mk;
			}
		}
		LU(tempC, AN, AU, AD, false);
		loops = PowERMethod(tempC, AN, MAXLOOPS, &m[j - 1], false);
		m[j - 1] += mk;
		printf("(%d)\t最接近%1.11e的特征值為%1.11e\t迭代次數: %d\n", j, mk, m[j - 1], loops);
	}
	
	//求譜范數和行列式的值
	for (i = 0; i < AN * AM; i++)
	{
		tempC[i] = c[i];
	}
	cond = fabs(r1/rs);
	det = LU(tempC, AN, AU, AD, true);
	printf("A的條件數cond(A)2 = %1.11e\n", cond);
	printf("A的行列式detA = %1.11e\n", det);
}

//冪法求最大特征值; type = true冪法; type = false反冪法
int PowERMethod(double *pc, int length, int endLoops, double *resultOut, int type)
{
	double y[AN];
	double u[AN];
	double h;
	double b1=0; 
	double b2;
	int i, k, t; 
	double aii[AN];
	for (i = 0; i < length; i++)
	{
		aii[i] = pc[AU * length + i];
	}	
	//u0設置初值
	for (i = 0; i < length; i++)
	{
		u[i] = 1;
	}
	
	for ( t = 1; t <= endLoops; t++)
	{
		k = MaxNumInVector(u, length);
		h = u[k];
		UnitaryVector(u, y, length, 1/fabs(h));

		if (type)
		{
			AMultiY(aii, y, u, length);
		}
		else
		{	
			Doolittle(pc, y, u, length, AU, AD);  	//y的值會被改變
		}

		b2 = GetSgn(h) * u[k];
		*resultOut = b2;
		if (!type)
		{
			*resultOut = 1 / b2;
		}
		
		if (t >1 && (fabs(b2 - b1) / fabs(b2) <= ER))
		{
			break;
		}
		b1 = b2;
	}
	return t;
}

//求向量中模最大的分量，返回分量的下標
int MaxNumInVector(double *pIn, int length)
{
	int i, k;
	double temp = pIn[0];
	k = 0;
	for (i = 0; i < length; i++)
	{
		if (fabs(pIn[i]) > fabs(temp))
		{
			temp = pIn[i];
			k = i;
		}
	}
	return k;
}

//矩陣A乘以向量y，結果保存在u中
void AMultiY(double *aii, double *inY, double *outU, int length)
{
	int i, j;
	for (i = 0; i < length; i++)
	{
		outU[i] = aii[i] * inY[i];
		j = i;               
		if (j - AD >= 0)
			outU[i] += cc * inY[j - AD];
		if (j - AD + 1 >= 0)
			outU[i] += bb * inY[j - AD + 1];	
		if (j + AU < length)
			outU[i] += cc * inY[j + AU];
		if (j + AU - 1 < length)
			outU[i] += bb * inY[j + AU - 1];
	}
}

//向量數乘,用于向量u歸一化
void UnitaryVector(double *pVectorIn, double *pVectorOut, int length, double numbER)
{
	int i;
	for (i = 0; i < length; i++)
	{
		pVectorOut[i] = pVectorIn[i] * numbER;
	}
}

//杜利特爾分解方法求方程組的解,結果儲存在b中
void Doolittle(double *c, double *b, double *xOut, int length, int s, int r)
{
	int i, t, max, min;
	double temp;
	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];
	}
}

//LU分解，如果flag=true，返回計算的行列式值
double LU(double *z, int length, int s, int r, int flag)
{
	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 += ( z[(k - t + s) * length + t - 1] * z[(t - j + s) * length + j - 1] );
			}
			z[(k - j + s) * length + j - 1] = z[(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 += ( z[(i - t + s) * length + t - 1] * z[(t - k + s) * length + k - 1] );
			}
			z[(i - k + s) * length + k - 1] = (z[(i - k + s) * length + k - 1] - temp) / z[s * length + k - 1];
		}
	}

	if (flag)
	{
		temp =	1;
		for (i = 0; i < length; i ++)
		{
			temp *= z[s * length + i];
		}
		return temp;
	}
	return 0;
}