/* ThresholdJacobi.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,eps;
	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();

	// しきい値の設定
	eps=0;
	for(i=0; i<SIZE; i++)
		for(j=i+1; j<SIZE; j++)
			eps+=fabs(a.ent[i][j]);
	eps=((2.0*eps)/SIZE)/(SIZE-1);

	for(k=1;k<=steps;k++){
		for(i=0; i<SIZE; i++){
			for(j=i+1; j<SIZE; j++){
				if(fabs(a.ent[i][j])>=eps){
					// 回転角度設定
					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));
				}
			}
		}
		eps/=10.0;
		printf("%d-th step : \n",k);
		printf("A_%d :\n",k);
		matprint(a);
		printf("\n");
		printf("P_%d :\n",k);
		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");
	}
}