#include <iostream>
#include <math.h>
#include "main.h"
using namespace std;

void printArray(double a[N]){
	for(int i=0;i<N;i++){
		printf("%-16.8e ",a[i]);
		if((i-3)%4==0)printf("\n");
	}
}

void setA(Matrix *&a){
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			a->setElement(1./(i+j+1.),i,j);
}

double getLamda(double a[N]){
	int i;
	double lmd=a[0];
	for(i=0;i<N;i++)
		if(lmd<a[i])lmd=a[i];
	return lmd;
}

double *unitary(double a[N],double l){
	int i;
	double *r=new double[N];
	for(i=0;i<N;i++)r[i]=a[i]/l;
	return r;
}

double *vectorMinus(double a[N],double b[N]){
	double *r=new double[N];
	for(int i=0;i<N;i++){
		r[i]=a[i]-b[i];
	}
	return r;
}

double mabs(double a[N]){				//向量范數,采用1范數 
	double r=0;
	for(int i=0;i<N;i++){
		r+=fabs(a[i]);
	}
	return r;
}
/* 
double *linearEquation(Matrix *A,double Y[N]){
	//A是系數矩陣,Y為常數項,采用Seidel迭代法 
	double *X=new double[N];
	double XX[N],T[N],abs=1;
	double *XXX=new double[N];
	int i,j,k=1;
	for(i=0;i<N;i++)X[i]=1;		//取初始向量[1,1,1,...,1] 
	while(abs>EPS&&k<100000){ 
		for(i=0;i<N;i++)X[i]=XX[i];
		for(i=0;i<N;i++){
			T[i]=0;
			for(j=0;j<i;j++)T[i]+=A->getElement(i,j)*XX[j];
			for(j=i+1;j<N;j++)T[i]+=A->getElement(i,j)*X[j];
			XX[i]=(Y[i]-T[i])/A->getElement(i,i);			
		}	
		XXX=vectorMinus(XX,X);
		abs=mabs(XXX);
	cout << k++ << endl;
	}
	delete [] XXX;
	return X;
}*/

double *linearEquation(Matrix *A,double b[N]){
	//LU分解
	//先解出Y
	double Y[N],s,*X=new double[N];
	int i,j;
	Matrix *L,*U;
	A->LU(L,U);
	Y[0]=b[0];
	for(i=1;i<N;i++){
		s=0;
		for(j=0;j<i;j++)s+=L->getElement(i,j)*Y[j];
		Y[i]=b[i]-s;
	}//計算Y	
	X[N-1]=Y[N-1]/U->getElement(N-1,N-1); 
	for(i=N-2;i>=0;i--){
		s=0;
		for(j=N-1;j>i;j--)s+=U->getElement(i,j)*X[j];
		X[i]=(Y[i]-s)/U->getElement(i,i);
	}
	delete L;
	delete U;
	return X;
} 

double evalue(Matrix *A,int &k1){
	double *X=new double[N];
	double *Y=new double[N];
	double mu=0,lamda;
	for(int i=0;i<N;i++)X[i]=1;
	lamda=getLamda(X);

	while(fabs(lamda-mu)>EPS){
		mu=lamda;	
		Y=unitary(X,lamda);
		X=(*A)*Y;
		lamda=getLamda(X);
		++k1;
	}
	delete [] X;
	delete [] Y;
	return lamda;
}

/*double re_evalue(Matrix *A,int &k2){
	double *X=new double[N];
	double *Y=new double[N];
	double mu=0,lamda;
	for(int i=0;i<N;i++)X[i]=1;
	lamda=getLamda(X);

	while(fabs((lamda-mu)/mu)>EPS){
		mu=lamda;	
		Y=unitary(X,lamda);
		X=linearEquation(A,Y);
		lamda=getLamda(X);
		++k2;
	}
	delete [] X;
	delete [] Y;
	return lamda;
}//收斂太慢,采用加速迭代*/

double re_evalue(Matrix *A,int &k2){
	double *X=new double[N];
	double *Y=new double[N];
	Matrix *I=new Matrix(N,1.0,1);
	double mu=1,lamda=0,alpha;
	for(int i=0;i<N;i++)X[i]=2;
	alpha=getLamda(X);
	Matrix lmdI=*I*lamda;
	Matrix AA=*A-lmdI;

	while(fabs(1./alpha-1./mu)>EPS){
		mu=alpha;	
		Y=unitary(X,alpha);
		X=linearEquation(&AA,Y);
		alpha=getLamda(X);
		++k2;
	}
	delete [] X;
	delete [] Y;
	delete I;
	return lamda+1./alpha;
}

int main(int argc, char *argv[]){
	int i,k1=0,k2=0;
	double ev1,ev2;
	Matrix *A=new Matrix(N);
	setA(A);//置初始矩陣A					 

	ev1=evalue(A,k1);
	ev2=re_evalue(A,k2);

	printf("矩陣A(i,j)=1/(i+j+1)\n\n");
	printf("A矩陣的按模最大特征值:");
	printf("%15.8e\n",ev1);
	printf("迭代次數:%d\n",k1);
	printf("\n");

	printf("A矩陣的按模最小特征值:");
	printf("%15.8e\n",ev2);
	printf("迭代次數:%d\n",k2);
	printf("\n");

	delete A;	
	return 0;
}