/* SerialJacobi.c
シリアル・ヤコビ法

必要な行･列だけの更新をすればより高速になる

2007.1.16
*/

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include<limits.h>
#include<string.h>

#define RANDOMIZE()	srand(time(NULL))
#define RANDOM(x)	(rand()%(x))

#include "sqmatrix.h"

main()
{
	int i,j,k,steps;
	MATRIX a,p,r,a0;
	double x, pi=3.14159265359,e=1e-10;
	
	SIZE=5;   // 行列サイズ
	steps=6;  // 反復回数

	RANDOMIZE();
	a=randomintmat(10);      // ランダムな整数行列
	a=matadd(a,transmat(a)); // ランダムな整数対称行列
	a0=a;
	printf("A :\n");
	matprint(a);
	printf("\n");
	p=unitmat();

	for(k=1;k<=steps;k++){
		for(i=0; i<SIZE; i++){
			for(j=i+1; j<SIZE; j++){
				// 回転角度設定
				if(fabs(a.ent[i][i]-a.ent[j][j])<e){
					x=pi/4.0;
				}
				else{
					x=atan(2*a.ent[i][j]/(a.ent[i][i]-a.ent[j][j]))/2.0;
				}
				// 回転の行列
				r=rotmat(i,j,x);
				// 直交行列の更新
				p=matmul(p,r);
				// A の更新
				a=matmul(transmat(r),matmul(a,r));
			}
		}
		if(k==1){
			printf("1st round : \n");
		} 
		else{
			if(k==2){
				printf("2nd round : \n");
			}
			else{
				if(k==3) printf("3rd round : \n");
				else printf("%dth round : \n",k);
			}
		}
		printf("A :\n");
		matprint(a);
		printf("\n");
		printf("P :\n");
		matprint(p);
		printf("\n");
		//getchar();
	}
	// 検算
	printf("Check\n\n");
	printf("P^{-1} A P :\n");
	matprint(matmul(transmat(p),matmul(a0,p)));
	printf("\n");
	printf("Hit any key\n");
	getchar();

	for(i=0;i<SIZE;i++){
		printf("%d-th eigenvalue and eigenvector :\n\n",i+1);
		printf("%12.7f\n\n",a.ent[i][i]);
		for(j=0;j<SIZE;j++) printf("%12.7f\n",p.ent[j][i]);
		printf("\n");
	}
}