実対称行列の固有値・固有ベクトル（ヤコビ法）

/************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/*      coded by Y.Suganuma                     */
/************************************************/
#include <stdio.h>

int Jacobi(int, int, double, double **, double **, double **, double **, double **);

int main()
{
	double **A, **A1, **A2, **X1, **X2, eps;
	int i1, i2, ind, ct, n;
					// データの設定
	ct   = 1000;
	eps  = 1.0e-10;
	n    = 3;
	A    = new double * [n];
	A1   = new double * [n];
	A2   = new double * [n];
	X1   = new double * [n];
	X2   = new double * [n];
	for (i1 = 0; i1 < n; i1++) {
		A[i1]  = new double [n];
		A1[i1] = new double [n];
		A2[i1] = new double [n];
		X1[i1] = new double [n];
		X2[i1] = new double [n];
	}

	A[0][0] = 1.0;
	A[0][1] = 0.0;
	A[0][2] = -1.0;

	A[1][0] = 0.0;
	A[1][1] = 1.0;
	A[1][2] = -1.0;

	A[2][0] = -1.0;
	A[2][1] = -1.0;
	A[2][2] = 2.0;
					// 計算
	ind = Jacobi(n, ct, eps, A, A1, A2, X1, X2);
					// 出力
	if (ind > 0)
		printf("収束しませんでした！\n");
	else {
		printf("固有値\n");
		for (i1 = 0; i1 < n; i1++)
			printf(" %f", A1[i1][i1]);
		printf("\n固有ベクトル（各列が固有ベクトル）\n");
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++)
				printf(" %f", X1[i1][i2]);
			printf("\n");
		}
	}

	for (i1 = 0; i1 < n; i1++) {
		delete [] A[i1];
		delete [] A1[i1];
		delete [] A2[i1];
		delete [] X1[i1];
		delete [] X2[i1];
	}
	delete [] A;
	delete [] A1;
	delete [] A2;
	delete [] X1;
	delete [] X2;

	return 0;
}

/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
/*      n : 次数                                             */
/*      ct : 最大繰り返し回数                                */
/*      eps : 収束判定条件                                   */
/*      A : 対象とする行列                                   */
/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/*      return : =0 : 正常                                   */
/*               =1 : 収束せず                               */
/*      coded by Y.Suganuma                                  */
/*************************************************************/
#include <math.h>

int Jacobi(int n, int ct, double eps, double **A, double **A1, double **A2,
           double **X1, double **X2)
{
	double max, s, t, v, sn, cs;
	int i1, i2, k = 0, ind = 1, p = 0, q = 0;
					// 初期設定
	for (i1 = 0; i1 < n; i1++) {
		for (i2 = 0; i2 < n; i2++) {
			A1[i1][i2] = A[i1][i2];
			X1[i1][i2] = 0.0;
		}
		X1[i1][i1] = 1.0;
	}
					// 計算
	while (ind > 0 && k < ct) {
						// 最大要素の探索
		max = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++) {
				if (i2 != i1) {
					if (fabs(A1[i1][i2]) > max) {
						max = fabs(A1[i1][i2]);
						p   = i1;
						q   = i2;
					}
				}
			}
		}
						// 収束判定
							// 収束した
		if (max < eps)
			ind = 0;
							// 収束しない
		else {
								// 準備
			s  = -A1[p][q];
			t  = 0.5 * (A1[p][p] - A1[q][q]);
			v  = fabs(t) / sqrt(s * s + t * t);
			sn = sqrt(0.5 * (1.0 - v));
			if (s*t < 0.0)
				sn = -sn;
			cs = sqrt(1.0 - sn * sn);
								// Akの計算
			for (i1 = 0; i1 < n; i1++) {
				if (i1 == p) {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == p)
							A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                       2.0 * A1[p][q] * sn * cs;
						else if (i2 == q)
							A2[p][q] = 0.0;
						else
							A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
					}
				}
				else if (i1 == q) {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == q)
							A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                       2.0 * A1[p][q] * sn * cs;
						else if (i2 == p)
							A2[q][p] = 0.0;
						else
							A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
					}
				}
				else {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == p)
							A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
						else if (i2 == q)
							A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
						else
							A2[i1][i2] = A1[i1][i2];
					}
				}
			}
								// Xkの計算
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++) {
					if (i2 == p)
						X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
					else if (i2 == q)
						X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
					else
						X2[i1][i2] = X1[i1][i2];
				}
			}
								// 次のステップへ
			k++;
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++) {
					A1[i1][i2] = A2[i1][i2];
					X1[i1][i2] = X2[i1][i2];
				}
			}
		}
	}

	return ind;
}