homework1.cpp

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

/***********************************************************************

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

#define  MAXLOOPS	2000	//最大迭代次数
#define GetSgn(a)	((a) > 0 ? 1 : ((a) < 0 ? -1 : 0))	//符号函数
#define  GetAbs(a)	(((a) > 0) ? (a) : ((a) * (-1)))	//计算a的绝对值
#define max3(a, b, c)	(((a) > (b) ? (a) : (b)) > (c) ? ((a) > (b) ? (a) : (b)) : (c))	//求三个数的最大值
#define min2(a, b)	((a) < (b) ? (a) : (b))	//求两个数的最小值

/***********************************************************************
 * 函数名称：
 * Ay()
 *
 * 参数：
 * double *aii,矩阵A指针
 * double *inY,列向量y指针
 * double *outU,乘积结果向量u指针
 * int length,列向量y的长度
 *
 * 返回值：
 * int 计算成功返回1，否则返回0
 *
 * 说明：
 * 矩阵A乘以向量y，结果保存在u中
 **********************************************************************/
int Ay(double *aii, double *inY, double *outU, int length);	
/***********************************************************************
 * 函数名称：
 * GetMaxInVector()
 *
 * 参数：
 * double *pIn,输入的列向量指针
 * int length,列向量的长度
 *
 * 返回值：
 * int 模最大的分量的下标
 *
 * 说明：
 * 求向量中模最大的分量，返回分量的下标，计算无穷范数时使用
 **********************************************************************/
int GetMaxInVector(double *pIn, int length);
/***********************************************************************
 * 函数名称：
 * VectorMNumber()
 *
 * 参数：
 * double *pVectorInOut,输入的列向量指针
 * double *pVectorOut,输出的结果向量指针
 * int length,列向量的长度
 * double number,乘的系数
 *
 * 返回值：
 * int 计算成功返回1，否则返回0
 *
 * 说明：
 * 向量数乘,用于向量u归一化
 **********************************************************************/
int VectorMNumber(double *pVectorInOut, double *pVectorOut, int length, double number);
/***********************************************************************
 * 函数名称：
 * MethodPower()
 *
 * 参数：
 * double *pc,如果是幂法则为压缩矩阵C的指针，如果是反幂法则为对C进行LU分解后的压缩矩阵的指针
 * int length,矩阵的宽度
 * int endLoops,最大迭代次数
 * double *resultOut,输出的特征值的指针
 * bool type,如果true则执行幂法，如果为false则执行反幂法
 *
 * 返回值：
 * int 返回实际迭代的次数
 *
 * 说明：
 * 使用幂法或反幂法求模最大或模最小特征值
 **********************************************************************/
int MethodPower(double *pc, int length, int endLoops, double *resultOut, bool type);
/***********************************************************************
 * 函数名称：
 * LU()
 *
 * 参数：
 * double *c,输入要分解的矩阵C的指针
 * int length,矩阵的宽度
 * int s,上半带宽
 * int r,下半带宽
 * bool det,如果true则计算行列式的值并返回
 *
 * 返回值：
 * double 如果det=true，返回计算的行列式值
 *
 * 说明：
 * 对输入压缩矩阵C进行LU分解
 **********************************************************************/
double LU(double *c, int length, int s, int r, bool det);
/***********************************************************************
 * 函数名称：
 * Doolittle()
 *
 * 参数：
 * double *c,输入压缩矩阵C的指针
 * double *b,输入方程右侧列向量的指针
 * double *xOut,输出计算出的方程解向量的指针
 * int length,矩阵的宽度
 * int s,上半带宽
 * int r,下半带宽
 *
 * 返回值：
 * double 如果det=true，返回计算的行列式值
 *
 * 说明：
 * 杜利特尔分解方法求带状线性方程组的解，b的值会改变
 **********************************************************************/
int Doolittle(double *c, double *b, double *xOut, int length, int s, int r);

//全局变量
double *c = 0;//用c压缩存储带状矩阵A
int as = 2;//上半带宽
int ar = 2;//下半带宽
int am = 5;//压缩矩阵高度
int an = 501;//压缩矩阵宽度
double er = pow(10, -12);//精度水平
double bb = 0.16;
double cc = -0.064;

//主函数
int main()
{
	int i, j, loops;
	double r1, rs, r501, d, mk, cond, det;
	double m[39];
	double * tempC = new double[an * am];
	//把带状矩阵A压缩存储为C
	c = new double[am * an];
	for (i = 0; i < an; i++)
	{
		c[0 * an + i] = cc;
		c[1 * an + i] = bb;
		c[as * 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 = MethodPower(tempC, an, MAXLOOPS, &r1, true);
	printf("最小特征值r1 = %1.11e 迭代次数: %d\n", r1, loops);
	//幂法求最大特征值，原点平移r501
	for (i = 0; i < an * am; i++)
	{
		tempC[i] = c[i];
		if (i >= as * an && i < (as + 1) * an)
		{
			tempC[i] -= r1;
		}
	}
	loops = MethodPower(tempC, an, MAXLOOPS, &r501, true);
	r501 += r1;
	printf("最大特征值r501 = %1.11e 迭代次数: %d\n", r501, loops);
	//反幂法求模最小特征值rs
	for (i = 0; i < an * am; i++)
	{
		tempC[i] = c[i];
	}
	LU(tempC, an, as, ar, false);
	loops = MethodPower(tempC, an, MAXLOOPS, &rs, false);
	printf("模最小特征值rs = %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 >= as * an && i < (as + 1) * an)
			{
				tempC[i] -= mk;
			}
		}
		LU(tempC, an, as, ar, false);
		loops = MethodPower(tempC, an, MAXLOOPS, &m[j - 1], false);
		m[j - 1] += mk;
		printf("(%d)最接近%1.11e的特征值为%1.11e 迭代次数: %d\n", j, mk, m[j - 1], loops);
	}
	//求谱范数和行列式的值
	for (i = 0; i < an * am; i++)
	{
		tempC[i] = c[i];
	}
	cond = GetAbs(r1/rs);
	det = LU(tempC, an, as, ar, true);
	printf("A的条件数cond(A)2 = %1.11e\n", cond);
	printf("A的行列式detA = %1.11e\n", det);

	delete[] tempC;
	delete[] c;
	return 0;
}

//幂法求最大特征值; type = true幂法; type = false反幂法
int MethodPower(double *pc, int length, int endLoops, double *resultOut, bool type)
{
	double * y = new double[length];
	double * u = new double[length];
	double * y2b = new double[length];
	double h;
	int i, k;
	//
	double * aii = new double[length];
	for (i = 0; i < length; i++)
	{
		aii[i] = pc[as * length + i];
	}	
	//u0设置初值
	for (i = 0; i < length; i++)
	{
		u[i] = 1;
	}
	double b1, b2;
	for (int t = 1; t <= endLoops; t++)
	{
		k = GetMaxInVector(u, length);
		h = u[k];
		VectorMNumber(u, y, length, 1/GetAbs(h));

		if (type)
		{
			Ay(aii, y, u, length);
		}
		else
		{
			//y的值会被改变，所以这步之后不能用y了
			Doolittle(pc, y, u, length, as, ar);
		}

		b2 = GetSgn(h) * u[k];
		*resultOut = b2;
		if (!type)
		{
			*resultOut = 1 / b2;
		}
		
		if (t >1 && (GetAbs(b2 - b1) / GetAbs(b2) <= er))
		{
			break;
		}
		b1 = b2;
	}
	delete[] aii;
	delete[] y;
	delete[] u;

	return t;
}

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

//矩阵A乘以向量y，结果保存在u中
int Ay(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 - ar >=0)
			outU[i] += cc * inY[j - ar];
		if (j - ar + 1 >= 0)
			outU[i] += bb * inY[j - ar + 1];	
		if (j + as < length)
			outU[i] += cc * inY[j + as];
		if (j + as - 1 < length)
			outU[i] += bb * inY[j + as - 1];
	}
	return 1;
}

//向量数乘,用于向量u归一化
int VectorMNumber(double *pVectorIn, double *pVectorOut, int length, double number)
{
	int i;
	for (i = 0; i < length; i++)
	{
		pVectorOut[i] = pVectorIn[i] * number;
	}
	return 1;
}

//杜利特尔分解方法求方程组的解，会改变b的内容
int 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];
	}
	return 1;
}

//LU分解，如果det=true，返回计算的行列式值
double LU(double *c, int length, int s, int r, bool det)
{
	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];
		}
	}

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

	return 0;
}