/* SerialJacobi.c
シリアル・ヤコビ法
2013年度 ver.
*/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

// 行列のサイズは大域変数としてここで宣言する。
int size = 4;

////////// マクロ等の定義 //////////

#define trace 0	// 途中経過表示 -> 1 / 非表示 -> 0 に設定せよ
#define N 128	// 行列の最大サイズ
#define RANDOM(n) (rand()%(2*n+1)-n)			// 絶対値 |n| 以下の整数乱数
#define REALRANDOM(x) (x*((2.0*rand()/RAND_MAX)-1))	// 絶対値 |x| 以下の実数乱数

////////// 関数のプロトタイプ宣言 //////////

// 行列 a を出力する関数
void matprint(double a[][N]);

// ベクトル v を出力する関数
void vecprint(double v[]);

// 単位行列を p に格納する関数
void unitmat(double p[][N]);

// (i,j)-座標に関する角度 x の回転の行列を r に格納する関数
void rotmat(int i, int j, double x, double r[][N]);

// 行列 a, b の積 ab を c に格納する関数
void matmul(double a[][N], double b[][N], double c[][N]);

// 行列 a の転置行列を b に格納する関数
void transpose(double a[][N], double b[][N]);

// 序数の出力
void ordinal(int k);

////////// メインルーティン //////////

main()
{
	double a[N][N];		// 実対称行列 A を表す２次元配列
	double a0[N][N];	// A のコピーを格納する２次元配列
	double p[N][N];		// A を対角化する直行行列 P を表す２次元配列
	double tp[N][N];	// P の転置行列を表す２次元配列
	double r[N][N];		// 回転の行列 R を表す２次元配列
	double tr[N][N];	// R の転置行列を表す２次元配列
	int steps = 8;		// 反復回数
	double e = 1e-10;	// 許容誤差
	int i, j, k;
	double x;

	srand(time(NULL));	// 乱数の初期化

	// ランダムな整数対称行列
	for(i = 0; i < size; i++)
		for(j = i; j < size; j++){
			x = RANDOM(10);
			a[i][j] = x;
			a[j][i] = x;
		}
	// A のコピー
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			a0[i][j] = a[i][j];
	printf("A :\n"); matprint(a); printf("\n");

	// P の初期値は単位行列
	unitmat(p);

	for(k = 1; k <= steps; k++){
		for(i = 0; i < size; i++){
			for(j = i + 1; j < size; j++){
				// 回転角度設定
				if(fabs(a[i][i] - a[j][j]) < e)	x = M_PI / 4.0;
				else x = atan(2 * a[i][j] / (a[i][i] - a[j][j])) / 2.0;
				// 回転の行列
				rotmat(i, j, x, r);
				// 直交行列の更新
				matmul(p, r, p);
				// A の更新
				transpose(r, tr);	// tr = (R の転置行列)
				matmul(a, r, a);
				matmul(tr, a, a);	// 新 A = tR A R
				if(trace){
					printf("A :\n"); matprint(a); printf("\n");
					printf("P :\n"); matprint(p); printf("\n");
					getchar();
				}
			}
		}
		ordinal(k); printf(" round : \n",k);
		printf("A :\n"); matprint(a); printf("\n");
		printf("P :\n"); matprint(p); printf("\n");
	}
	// 検算
	transpose(p, tp);	// tp = (P の転置行列)
	matmul(a0, p, a0);
	matmul(tp, a0, a0);	// 新 A0 = tP A0 P
	printf("Check\n\n");
	printf("P^{-1} A P :\n"); matprint(a0); printf("\n");
	printf("Hit any key\n");
	getchar();

	for(i = 0; i < size; i++){
		ordinal(i + 1); printf(" eigenvalue and eigenvector :\n\n");
		printf("%12.7f\n\n", a[i][i]);
		for(j = 0; j < size; j++) printf("%12.7f\n", p[j][i]);
		printf("\n");
	}
}

////////// 関数定義 //////////

// 行列 a を出力する関数
void matprint(double a[][N])
{
	int i, j;
	for(i = 0; i < size; i++){
		for(j = 0; j < size; j++)
			printf("%12.7lf ", a[i][j]);
		printf("\n");
	}
	return;
}

// ベクトル b を出力する関数
void vecprint(double b[])
{
	int i;
	for(i = 0; i < size; i++)
		printf("%12.7lf\n", b[i]);
	return;
}

// 単位行列を p に格納する関数
void unitmat(double p[][N])
{
	int i, j;
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			p[i][j] = 0.;
	for(i = 0; i < size; i++)
		p[i][i] = 1.;
	return;
}

// (i,j)-座標に関する角度 x の回転の行列を r に格納する関数
void rotmat(int i, int j, double x, double r[][N])
{
	double c = cos(x), s = sin(x);
	unitmat(r);
	r[i][i] = c;
	r[i][j] = -s;
	r[j][i] = s;
	r[j][j] = c;
	return;
}

// 行列 a, b の積 ab を c に格納する関数
void matmul(double a[][N], double b[][N], double c0[][N])
{
	int i, j, k;
	double x;
	double c[N][N];	// a や b が c0 と同一でも構わないように別の配列を用意
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++){
			x = 0;
			for(k = 0; k < size; k++) x += a[i][k] * b[k][j];
			c[i][j] = x;
		}
	for(i = 0; i < size; i++)
		for(j = 0; j < size; j++)
			c0[i][j] = c[i][j];
	return;
}

// 行列 a の転置行列を b に格納する関数
void transpose(double a[][N], double b[][N])
{
	int i, j;
	double x, y;
	for(i = 0; i < size; i++)
		for(j = i; j < size; j++){
			x = a[i][j];
			y = a[j][i];
			b[i][j] = y;
			b[j][i] = x;
		}
	return;
}

// 序数の出力
void ordinal(int k)
{	
	int m = k % 100;
	if((m < 11) || (m > 20)){
		switch(k % 10){
			default: printf("%dth", k); break;
			case 1: printf("%dst", k); break;
			case 2: printf("%dnd", k); break;
			case 3: printf("%drd", k); break;
		}
	}
	else printf("%dth", k); 
	return;
}