幂法与反幂法 C++程序

这是一个用C++写的幂法和反幂法算法，基本没有特别的东西

//描述：用来用幂法和反幂法来解矩阵的特征值与特征向量
//输入：一个用一维数组表示的矩阵、矩阵的阶数
//输出：矩阵的某些特征值与特征向量

#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