//描述：用來用冪法和反冪法來解矩陣的特征值與特征向量
//輸入：一個用一維數組表示的矩陣、矩陣的階數
//輸出：矩陣的某些特征值與特征向量

#include?
#include?
#include?
#include?
#include?“Power_Method.h“

using?namespace?std;

const??double?EPS=?1e-6;

template
void?InitData（int?&N?valuetype?*matrix?char?*filename）
{
	//從文件filename里讀入matrix
	fstream?infile（filename?ios_base::in）;
	for（int?i=0;?i		infile?>>?matrix[i];

	infile.close（）;
}

template
valuetype?VectorMul（valuetype?*a?valuetype?*b?int?&count）
{
	valuetype?res=0;
	for（int?i=0;?i		res?+=?a[i]*b[i];
	return?res;
}

//計算矩陣與向量的積
//參數：outVec計算結果
template
void?MatrixMulVector（valuetype?*outVec?valuetype?*matrix?valuetype?*vec?int&count）
{
	for（int?row=0;?row	{
		valuetype?sum?=?0;
		for（int?col=0;?col			sum?+=?matrix[row*count?+?col]?*?vec[col];
		outVec[row]?=?sum;
	}
}

//冪法主函數
//參數：outEigenvalue--特征值、outEigenvector--特征向量、matrix--輸入的矩陣、order--矩陣的階、center--輸入的原點平移量
template
void?PowerMethod（valuetype?&outEigenvalue?valuetype?*outEigenvector?valuetype?*matrix?int?&order?valuetype?¢er）
{
	//初始化一些臨時需要的變量
	valuetype?*u0?=?new?valuetype[order];
	valuetype?*u1?=?new?valuetype[order];
	valuetype?*y0?=?new?valuetype[order];
	valuetype?beta0=0?beta1=0;
	size_t?size?=?sizeof（u0）;
	memset（u0?0?order*size）;

	//初始化一個非零的向量u0
	u0[0]?=?1;

	//平移矩陣
	for（int?i=0;?i		matrix[i*order+i]?-=center;

	//****************************************************
	//冪法主過程
	do?
	{
		beta0?=?beta1;		//記錄上次循環之后beta的值
		memcpy（u0?u1?order*sizeof（valuetype））;

		valuetype?lengthOfU?=?sqrt（VectorMul（u0?u0?order））;
		for（int?i=0;?i			y0[i]?=?u0[i]/lengthOfU;

		MatrixMulVector（u1?matrix?y0?order）;			//u1?=?matrix?*?y0

		beta1?=?VectorMul（y0?u1?order）;
	}?while?（abs（beta1-beta0）/abs（beta1）?>?EPS）;

	outEigenvalue?=?beta1?+?center;
	memcpy（outEigenvector?y0?order*sizeof（valuetype））;
	//****************************************************

	//銷毀臨時分配的空間
	delete（u0?u1?y0）;

	//恢復平移過的矩陣
	for（int?i=0;?i		matrix[i*order+i]?+=center;
}

//反冪法求解矩陣matrix的距離center最近的特征值
template
void?InversePowerMethod（valuetype?&outEigenvalue?valuetype?*outEigenvector?valuetype?*matrix?int?&order?valuetype?¢er）
{
	//初始化一些臨時需要的變量
	valuetype?*u0?=?new?valuetype[order];
	valuetype?*u1?=?new?valuetype[order];
	valuetype?*y0?=?new?valuetype[order];
	valuetype?beta0=0?beta1=0;
	size_t?size?=?sizeof（u0）;
	memset（u0?0?order*size）;

	//初始化一個非零的向量u0
	u0[0]?=?1;

	//平移矩陣
	for（int?i=0;?i		matrix[i*order+i]?-=center;

	//把matrix進行LU分解，以方便后面的解線性方程組
	valuetype?**LUMatrix?=?new?valuetype*[order];
	for（int?i=0;?i		LUMatrix[i]?=?new?valuetype[order];
	LUMethod（LUMatrix?matrix?order）;

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

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件???????10200??2010-11-16?14:01??Power_Method.cpp
?????文件????????1442??2010-10-10?23:20??Power_Method.h