szfx1.cpp

求矩阵的最大特征值

// SZFX1.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <memory.h>
#include <stdio.h>
#include <afx.h>
#include <math.h>
#include <iostream.h>

//对三角矩阵下半宽带
const int g_nrLow = 3;
//对三角矩阵上半宽带
const int g_nsUpp = 3;
//对三角矩阵维数
const int g_nMatrixDimension = 1000;
//对三角矩阵压缩矩阵
double g_cMatrix[7][1000];
//double size
const int g_ndoublesize = sizeof(double);
//记录当前矩阵参数K值
int g_nkVar = 0;
//存储条件数数组
double g_vCondition[81];
//存储模最大特征值数组
double g_vMaxFeature[81];
//存储模最小特征值数组
double g_vMinFeature[81];
//存储最大特征值数组
double g_vMaxFeat[81];

//初始化对三角矩阵的压缩矩阵
void InitMatrix(int tk);
//用幂法求条件数
double ComputeCondition();
//用第一种迭代格式的幂法求模最大的特征值
double ComputeMaxFeature1();
//用第二种迭代格式的幂法求模最大的特征值
double ComputeMaxFeature2();
//用反幂法求模最小的特征值
double ComputeMinFeatrue();
//求向量的2-范数
double ComputeTFanShu(double *p);
//求向量的无穷范数
double ComputeNFanShu(double *p, int *r);
//A--对三角矩阵 r--下半带宽 s--上半带宽
//幂法迭代公式 u = Ay，计算u的第i个分量
double ComputeVector(int i, double *pyVector);
//用doolittle分解法求解对三角矩阵: Au = y中的u向量, 参考Page40
bool ComputeUVectorUseDoolittle(double *pyVector, double *puVector);
//用带平移的反幂法计算最大特征值
double ComputeMax(double dbuOffside);
//初始化带平移反幂法矩阵的压缩矩阵
void InitMatrixEx(int tk, double dbuOffside);

int main(int argc, char* argv[])
{
	try
	{
		//初始化
		memset((void *)g_vCondition, 0, 81*g_ndoublesize);  
		memset((void *)g_vMaxFeature, 0, 81*g_ndoublesize);  
		memset((void *)g_vMinFeature, 0, 81*g_ndoublesize);  
		memset((void *)g_vMaxFeat, 0, 81*g_ndoublesize);  
		for(int i=0; i<=80; i++)
		{
			g_nkVar = i;
			//初始化三角矩阵的压缩矩阵
			InitMatrix(i);
			//计算条件式中 k = i 时的条件数
			g_vCondition[i] = ComputeCondition();
			//求最大特征值
			g_vMaxFeat[i] = ComputeMax(fabs(g_vMaxFeature[i]));
			printf(" k = %d : \n", g_nkVar);
			printf("     |Min|=%.6e  |Max|=%.6e  Cond(A)=%.6e  Max=%.6e\n\n", g_vMinFeature[g_nkVar], g_vMaxFeature[g_nkVar],
				g_vCondition[g_nkVar], g_vMaxFeat[i]);

		}
		double dbMin = 0, dbMax = 0, dbMaxFeat = 0;
		int nMin = 0, nMax = 0, nMaxFeat = 0;
		dbMin = dbMax = g_vCondition[0];
		dbMaxFeat = g_vMaxFeat[0];
		//获取最大和最小条件数以及相应输入值
		for(i=1; i<=80; i++)
		{
			if (g_vCondition[i]>dbMax)
			{
				dbMax = g_vCondition[i];
				nMax = i;
			}
			if (g_vCondition[i]<dbMin)
			{
				dbMin = g_vCondition[i];
				nMin = i;
			}
			if (g_vMaxFeat[i]>dbMaxFeat)
			{
				dbMaxFeat = g_vMaxFeat[i];
				nMaxFeat = i;
			}
		}
		//输出计算结果
		printf("Max(k)=%d  Max[Cond(A)]=%.6e\n", nMax, dbMax);
		printf("Min(k)=%d  Min[Cond(A)]=%.6e\n", nMin, dbMin);
		printf("Max(k)=%d  Max[Feat(A)]=%.6e\n", nMaxFeat, dbMaxFeat);
	}
	catch(...)
	{
	}
	return 0;
}
	
//初始化对三角矩阵的压缩矩阵
void InitMatrix(int tk)
{
	//矩阵初始化
	int r = 3,  s = 3;
	int n = 1000;
	memset(g_cMatrix, 0, (r+s+1)*n*g_ndoublesize);
	for(int ic=1; ic<=7; ic++)
	{
		for(int jc=1; jc<=1000; jc++)
		{
			int ja = jc;
			int ia = ic + jc - g_nsUpp - 1;
			if (((max(ja-3,1)<=ia && ia<ja) || (ja<ia && ia<=min(ja+3,1000))) && (ia>0 && ja>0))
			{
				g_cMatrix[ic-1][jc-1] = -1 * (ia + ja) / (1 + tk*0.05);
			}
			else if(ia == ja)
				g_cMatrix[ic-1][jc-1] = 100 * (ia + 0.1) / (ia + tk*0.05);
		}
	}
	//对三角矩阵A = LU分解的第一步，参看P40
	for(int k=1; k<=n; k++)
	{
		int jMax = min(k+s, n);
		double dbTemp = 0;
		for(int j=k; j<=jMax; j++)
		{
			int tMin = max(1, k-r);
			tMin = max(tMin, j-s);
			dbTemp = 0;
			for(int t=tMin; t<=k-1; t++)
				dbTemp += (g_cMatrix[k-t+s][t-1] * g_cMatrix[t-j+s][j-1]);
			g_cMatrix[k-j+s][j-1] -= dbTemp;
		}
		int iMax = min(k+r, n);
		for(int i=k+1; i<=iMax; i++)
		{
			int tMin = max(1, i-r);
			tMin = max(tMin, k-s);
			dbTemp = 0;
			for(int t=tMin; t<=k-1; t++)
				dbTemp += (g_cMatrix[i-t+s][t-1] * g_cMatrix[t-k+s][k-1]);
			dbTemp = g_cMatrix[i-k+s][k-1] - dbTemp;
			g_cMatrix[i-k+s][k-1] = dbTemp / g_cMatrix[s][k-1];
		}
	}
}

//求条件数
double ComputeCondition()
{
	//to compute max feature
	g_vMaxFeature[g_nkVar] = ComputeMaxFeature1();
	//to compute min feature
	g_vMinFeature[g_nkVar] = ComputeMinFeatrue();
	//to compute cond(A)
	return fabs(g_vMaxFeature[g_nkVar]/g_vMinFeature[g_nkVar]);
}

//用第一种迭代格式的幂法求模最大的特征值
double ComputeMaxFeature1()
{
	//to define two vectors : u, y
	double uVector[1000], yVector[1000];
	//to define one doulbe variant : beta;
	double dbPrevBeta = 0, dbBeta = 0;
	//to initialize uVector and yVector
	memset(uVector, 0, g_ndoublesize * 1000);
	memset(yVector, 0, g_ndoublesize * 1000);
	for(int j=0; j<1000; j++)
		uVector[j] = 1.00;
	bool bEnd = true;        //不懂bool
	double dbTemp = 0;
	do
	{
		double dbFanShu = ComputeTFanShu(uVector);
		for(int i=0; i<1000; i++)
			yVector[i] = uVector[i]/dbFanShu;
		for(i=1; i<=1000; i++)
			uVector[i-1] = ComputeVector(i, yVector);              //不懂i－1
		dbBeta = 0;  //初始化
		for(i=0; i<1000; i++)
			dbBeta += (yVector[i] * uVector[i]);
		bEnd = (fabs((dbBeta - dbPrevBeta) / dbBeta) <= 1e-7);
		dbPrevBeta = dbBeta;
	}while( !bEnd );         
	return dbBeta;
}

//用第二种迭代格式的幂法求模最大的特征值
double ComputeMaxFeature2()
{
	//to define two vectors : u, y
	double uVector[1000], yVector[1000];
	//to define one doulbe variant : beta;
	double dbPrevBeta = 0, dbBeta = 0;
	//to init uVector and yVector
	memset(uVector, 0, g_ndoublesize * 1000);
	memset(yVector, 0, g_ndoublesize * 1000);
	for(int j=0; j<1000; j++)
		uVector[j] = 1;
	bool bEnd = true;
	double dbTemp = 0;
	do
	{
		int r;
		double dbFanShu = ComputeNFanShu(uVector, &r);
		for(int i=0; i<1000; i++)
			yVector[i] = uVector[i]/dbFanShu;
		for(i=1; i<=1000; i++)
			uVector[i-1] = ComputeVector(i, yVector);
		double dbBeta = yVector[r] * uVector[r];
		bEnd = (fabs((dbBeta - dbPrevBeta) / dbBeta) <= 1e-7);
		dbPrevBeta = dbBeta;
	}while( !bEnd );
	return dbBeta;
}

//用反幂法求最小特征值
double ComputeMinFeatrue()
{
	//to define two vectors : u, y
	double uVector[1000], yVector[1000];
	//to define one doulbe variant : beta;
	double dbPrevBeta = 0, dbBeta = 0;
	//to init uVector and yVector
	memset(uVector, 0, g_ndoublesize * 1000);
	memset(yVector, 0, g_ndoublesize * 1000);
	for(int i=0; i<1000; i++)
		uVector[i] = 1.00;
	bool bEnd = true;                    //不懂1
	do
	{
		double dbFanShu = ComputeTFanShu(uVector);
		for(i=0; i<1000; i++)
			yVector[i] = uVector[i]/dbFanShu;
		//to compute Au = y : using doolittle method
		ComputeUVectorUseDoolittle(yVector, uVector);	
		dbBeta = 0;
		for(i=0; i<1000; i++)
			dbBeta += (yVector[i] * uVector[i]);
		double dbTemp = fabs((dbBeta - dbPrevBeta) / dbBeta );
		bEnd = fabs((dbBeta - dbPrevBeta) / dbBeta) <= 1e-7;
		dbPrevBeta = dbBeta;
	}while( !bEnd );
	dbBeta = 1 / dbBeta;
	return dbBeta;
}

//求向量的2-范数
double ComputeTFanShu(double *p)          //不懂2
{
	double dbTemp = 0;
	ASSERT(p!=NULL);                 //不懂3
	for(int i=0; i<1000; i++)
		dbTemp += (p[i]*p[i]);
	dbTemp = sqrt(dbTemp);
	return dbTemp;
}

//求向量的无穷范数
double ComputeNFanShu(double *p, int *r)
{
	ASSERT(p!=NULL);
	double dbRes = 0;
	for(int i=0; i<1000; i++)
	{
		if ((fabs(dbRes) - fabs(p[i])) < 0)
		{
			dbRes = fabs(p[i]);
			*r = i;
		}
	}
	return dbRes;
}

//A--对三角矩阵 r--下半带宽 s--上半带宽
//幂法迭代公式 u = Ay，计算u的第i个分量
double ComputeVector(int i, double *pyVector)
{
	ASSERT(i > 0);             //不懂assert
	double dbRes = 0;
	double dbTemp = 0;
	int jMin = 0, jMax = 0;
	if (i>=1 && i<=g_nrLow+1)
	{
		jMin = 1;
		jMax = g_nsUpp + i; 
	}
	else
	{
		jMin = min(i - 3, 1000);
		jMax = min(jMin + 6, 1000);
	}
	for(int j=jMin; j<=jMax; j++)
	{
		if (i==j)
			dbTemp = 100 * (i + 0.1) / (i + g_nkVar*0.05);
		else
			dbTemp = -1 * (i + j) / (1 + g_nkVar*0.05);
		dbRes += (dbTemp * pyVector[j-1]);
	}
	return dbRes;
}

//用doolittle分解法求解对三角矩阵: Au = y中的u向量, 参考Page40
bool ComputeUVectorUseDoolittle(double *pyVector, double *puVector)
{
	ASSERT(pyVector!=NULL && puVector!=NULL);
	int n = 1000, r = g_nrLow,  s = g_nsUpp;
	bool bRes = true;
	double dbTemp = 0;
	double yNewVector[1000];
	memcpy(yNewVector, pyVector, g_ndoublesize * n);
	for(int i=2; i<=n; i++)
	{
		dbTemp = 0;
		int tMin = max(1, i-r);
		for(int t=tMin; t<=i-1; t++)
			dbTemp += (g_cMatrix[i-t+s][t-1] * yNewVector[t-1]);
		yNewVector[i-1] -= dbTemp;
	}
	puVector[n-1] = yNewVector[n-1] / g_cMatrix[s+1-1][n-1];
	for(i=n-1; i>=1; i--)
	{
		dbTemp = 0;
		int tMax = min(i+s, n);
		for(int t=i+1; t<=tMax; t++)
			dbTemp += (g_cMatrix[i-t+s][t-1] * puVector[t-1]);
		dbTemp = yNewVector[i-1] - dbTemp;
		puVector[i-1] = dbTemp / g_cMatrix[s][i-1];
	}
	return bRes;
}

double ComputeMax(double dbuOffside)
{
	//to define two vectors : u, y
	double uVector[1000], yVector[1000];
	//to define one doulbe variant : beta;
	double dbPrevBeta = 0, dbBeta = 0;
	//初始化用带平移的反幂法求最大特征值的压缩矩阵
	InitMatrixEx(g_nkVar, fabs(g_vMaxFeature[g_nkVar]));
	//to init uVector and yVector
	memset(uVector, 0, g_ndoublesize * 1000);
	memset(yVector, 0, g_ndoublesize * 1000);
	for(int i=0; i<1000; i++)
		uVector[i] = 1.00;
	bool bEnd = true;
	do
	{
		double dbFanShu = ComputeTFanShu(uVector);
		for(i=0; i<1000; i++)
			yVector[i] = uVector[i]/dbFanShu;
		//to compute Au = y : using doolittle method
		ComputeUVectorUseDoolittle(yVector, uVector);	
		dbBeta = 0;
		for(i=0; i<1000; i++)
			dbBeta += (yVector[i] * uVector[i]);
		double dbTemp = fabs((dbBeta - dbPrevBeta) / dbBeta );
		bEnd = fabs((dbBeta - dbPrevBeta) / dbBeta) <= 1e-7;
		dbPrevBeta = dbBeta;
	}while( !bEnd );
	return 1/dbBeta + dbuOffside;
}

//初始化带平移反幂法矩阵的压缩矩阵,迭代公式中(A-uI)u = y,其中u=dbuOffside
void InitMatrixEx(int tk, double dbuOffside)
{
	//矩阵初始化
	int r = 3,  s = 3;
	int n = 1000;
	memset(g_cMatrix, 0, (r+s+1)*n*g_ndoublesize);
	for(int ic=1; ic<=7; ic++)
	{
		for(int jc=1; jc<=1000; jc++)
		{
			int ja = jc;
			int ia = ic + jc - g_nsUpp - 1;
			if (((max(ja-3,1)<=ia && ia<ja) || (ja<ia && ia<=min(ja+3,1000))) && (ia>0 && ja>0))
			{
				g_cMatrix[ic-1][jc-1] = -1 * (ia + ja) / (1 + tk*0.05);
			}
			else if(ia == ja)
				g_cMatrix[ic-1][jc-1] = 100 * (ia + 0.1) / (ia + tk*0.05) - dbuOffside;
		}
	}
	//对三角矩阵A = LU分解的第一步，参看P40
	for(int k=1; k<=n; k++)
	{
		int jMax = min(k+s, n);
		double dbTemp = 0;
		for(int j=k; j<=jMax; j++)
		{
			int tMin = max(1, k-r);
			tMin = max(tMin, j-s);
			dbTemp = 0;
			for(int t=tMin; t<=k-1; t++)
				dbTemp += (g_cMatrix[k-t+s][t-1] * g_cMatrix[t-j+s][j-1]);
			g_cMatrix[k-j+s][j-1] -= dbTemp;
		}
		int iMax = min(k+r, n);
		for(int i=k+1; i<=iMax; i++)
		{
			int tMin = max(1, i-r);
			tMin = max(tMin, k-s);
			dbTemp = 0;
			for(int t=tMin; t<=k-1; t++)
				dbTemp += (g_cMatrix[i-t+s][t-1] * g_cMatrix[t-k+s][k-1]);
			dbTemp = g_cMatrix[i-k+s][k-1] - dbTemp;
			g_cMatrix[i-k+s][k-1] = dbTemp / g_cMatrix[s][k-1];
		}
	}
}