用幂法和反幂法分别计算矩阵按模最大和按模最小的特征值及其特征向量

//
//	 用幂法和反幂法分别计算矩阵按模最大和按模最小的特征值及其特征向量	//
//
#include<stdio.h>
#include<math.h>
#include<iostream>
#define MAX_N 20	//方阵最大维数
#define MAXREPT 100	//迭代次数
#define epsilon 0.00001	//求解精度
using namespace std;

int n;
int i, j, k;
double error;
double xmax, oxmax;
static double a[MAX_N][MAX_N];		//AX=C的A矩阵
static double  L[MAX_N][MAX_N], u[MAX_N][MAX_N];
static double  x[MAX_N], nx[MAX_N];

void power_Method() {

	//初始化特征向量的值
	for (i = 0; i < n; i++) 
		x[i] = 1;

	oxmax = 0;                
	for (i = 0; i < MAXREPT; i++) {			//迭代最大周期

		for (j = 0; j < n; j++) {		//幂乘   x^(k+1) =  A  *   x^k 
			nx[j] = 0;

			for (k = 0; k < n; k++) {		//矩阵*向量
				nx[j] += a[j][k] * x[k];
			}
		}

		xmax = 0.0;
		for (j = 0; j < n; j++) 		//规范化  x[i] / max(x)
			if (fabs(nx[j]) > xmax)
				xmax = fabs(nx[j]);

		for (j = 0; j < n; j++)
			 nx[j] /= xmax;	

		for (j = 0; j < n; j++)
			x[j] = nx[j] ;

		if (fabs(xmax - oxmax) < epsilon) {
			cout << "Solve the problem after "<<i<<" repeats,the max lambda = " << xmax << endl;
			cout << "The vector is:" << endl;
			for (j = 0;j < n; j++)
				cout << nx[j] << endl;
			return;
		}
		oxmax = xmax;
	}
	cout << "Ater " << MAX_N << "repeat,max no result...\n" << endl;
}
void inversePowerMethod();
int main() {
	cout<<"输入矩阵的维数n（dim of AX =C)"<<endl;
	cin >> n;
	if (n > MAX_N) {
		cout << "The input n is larger than MAX_N,please redefine the MAX_N" << endl;
		return 1;
	}
	if (n <= 0) {
		cout << "Please input a proper number between 1 and" << MAX_N << endl;
		return 1;
	}

	//输入矩阵A
	printf("Now input the matrix a(i,j),i,j=0...%d\n", n - 1);
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			cin >> a[i][j];
	// power_Method();
	inversePowerMethod();
}

/*
	P(A-pI)=LU

*/
void inversePowerMethod() {
	for (i = 0; i < n; i++)		//	U矩阵对角元素为1		
		u[i][i] = 1;

	for (k = 0; k < n; k++) {
		for (i = k; i < n; i++) {		//计算L矩阵
			L[i][k] = a[i][k];
			for (j = 0; j <= k - 1; j++) 
				L[i][k] -= L[i][j] * u[j][k];			
		}

		for (j = k + 1; j < n; j++) {		//计算U矩阵
			u[k][j] = a[k][j];
			for (i = 0; i <= k - 1; i++)
				u[k][j] -= L[k][i] * u[i][j];
			u[k][j] /= L[k][k];
		}
	}

	for (i = 0; i < n; i++)
		x[i] = 1;

	oxmax = 0;

	for (i = 0; i < MAXREPT; i++) {
		for (j = 0; j < n; j++){		//反幂乘
			nx[j] = x[j];
			for (k = 0; k <= j - 1; k++)
				nx[j] -= L[j][k] * nx[k];
			nx[j] /= L[j][j];
		}

		for (j = n - 1; j >= 0; j--) {
			x[j] = nx[j];
			for (k = j + 1; k < n; k++) 
				x[j] -= u[j][k] * x[k];			
		}

		xmax = 0;
		for (j = 0; j < n; j++) 
			if (fabs(x[j] > xmax))
				xmax = fabs(x[j]);
		
		for (j = 0; j < n; j++) {
			x[j] /= xmax;
		}

		if (fabs(xmax - oxmax) < epsilon) {
			cout << "Solve .. min lamda=" << 1 / xmax << endl;
			for (i = 0; i < n; i++)
				cout << x[i] << endl;
			return;
		}
		oxmax = xmax;
	}
}

/*
	测试案例
	A=[ 1 1 0.5
		1 1 0.25 
		0.5 0.25 2.0]
*/