最大固有値と固有ベクトル（べき乗法）

/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      coded by Y.Suganuma             */
/****************************************/
#include <stdio.h>

int power(int, int, int, int, double, double **, double **, double **,
          double *, double **, double *, double *, double *);

int main()
{
	double **A, **B, **C, *r, **v, *v0, *v1, *v2, eps;
	int i1, i2, ind, ct, m, n;
					// データの設定
	ct  = 200;
	eps = 1.0e-10;
	n   = 3;
	m   = 15;

	A   = new double * [n];
	B   = new double * [n];
	C   = new double * [n];
	v   = new double * [n];
	for (i1 = 0; i1 < n; i1++) {
		A[i1] = new double [n];
		B[i1] = new double [n];
		C[i1] = new double [n];
		v[i1] = new double [n];
	}

	r  = new double [n];
	v0 = new double [n];
	v1 = new double [n];
	v2 = new double [n];

	A[0][0] = 11.0;
	A[0][1] = 6.0;
	A[0][2] = -2.0;

	A[1][0] = -2.0;
	A[1][1] = 18.0;
	A[1][2] = 1.0;

	A[2][0] = -12.0;
	A[2][1] = 24.0;
	A[2][2] = 13.0;

	for (i1 = 0; i1 < n; i1++) {
		for (i2 = 0; i2 < n; i2++)
			B[i1][i2] = 0.0;
		B[i1][i1] = 1.0;
		v0[i1]    = 1.0;
	}
					// 計算
	ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 出力
	if (ind == 0)
		printf("固有値が求まりませんでした！\n");
	else {
		for (i1 = 0; i1 < ind; i1++) {
			printf("固有値 %f", r[i1]);
			printf(" 固有ベクトル ");
			for (i2 = 0; i2 < n; i2++)
				printf(" %f", v[i1][i2]);
			printf("\n");
		}
	}

	for (i1 = 0; i1 < n; i1++) {
		delete [] A[i1];
		delete [] B[i1];
		delete [] C[i1];
		delete [] v[i1];
	}
	delete [] A;
	delete [] B;
	delete [] C;
	delete [] v;
	delete [] r;
	delete [] v0;
	delete [] v1;
	delete [] v2;

	return 0;
}

/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      i : 何番目の固有値かを示す      */
/*      n : 次数                        */
/*      m : 丸め誤差の修正回数          */
/*      ct : 最大繰り返し回数           */
/*      eps : 収束判定条件              */
/*      A : 対象とする行列              */
/*      B : nxnの行列，最初は，単位行列 */
/*      C : 作業域，nxnの行列           */
/*      r : 固有値                      */
/*      v : 各行が固有ベクトル（nxn）   */
/*      v0 : 固有ベクトルの初期値       */
/*      v1,v2 : 作業域（n次元ベクトル） */
/*      return : 求まった固有値の数     */
/*      coded by Y.Suganuma             */
/****************************************/
#include <math.h>

int power(int i, int n, int m, int ct, double eps, double **A, double **B,
          double **C, double *r, double **v, double *v0, double *v1, double *v2)
{
	double l1, l2, x;
	int i1, i2, i3, ind, k, k1;
					// 初期設定
	ind = i;
	k   = 0;
	l1  = 0.0;
	for (i1 = 0; i1 < n; i1++) {
		l1     += v0[i1] * v0[i1];
		v1[i1]  = v0[i1];
	}
	l1 = sqrt(l1);
					// 繰り返し計算
	while (k < ct) {
						// 丸め誤差の修正
		if (k%m == 0) {
			l2 = 0.0;
			for (i1 = 0; i1 < n; i1++) {
				v2[i1] = 0.0;
				for (i2 = 0; i2 < n; i2++)
					v2[i1] += B[i1][i2] * v1[i2];
				l2 += v2[i1] * v2[i1];
			}
			l2 = sqrt(l2);
			for (i1 = 0; i1 < n; i1++)
				v1[i1] = v2[i1] / l2;
		}
						// 次の近似
		l2 = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			v2[i1] = 0.0;
			for (i2 = 0; i2 < n; i2++)
				v2[i1] += A[i1][i2] * v1[i2];
			l2 += v2[i1] * v2[i1];
		}
		l2 = sqrt(l2);
		for (i1 = 0; i1 < n; i1++)
			v2[i1] /= l2;
						// 収束判定
							// 収束した場合
		if (fabs((l2-l1)/l1) < eps) {
			k1 = -1;
			for (i1 = 0; i1 < n && k1 < 0; i1++) {
				if (fabs(v2[i1]) > 0.001) {
					k1 = i1;
					if (v2[k1]*v1[k1] < 0.0)
						l2 = -l2;
				}
			}
			k    = ct;
			r[i] = l2;
			for (i1 = 0; i1 < n; i1++)
				v[i][i1] = v2[i1];
			if (i == n-1)
				ind = i + 1;
			else {
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						C[i1][i2] = 0.0;
						for (i3 = 0; i3 < n; i3++) {
							x          = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
							C[i1][i2] += x * B[i3][i2];
						}
					}
				}
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++)
						B[i1][i2] = C[i1][i2];
				}
				for (i1 = 0; i1 < n; i1++) {
					v1[i1] = 0.0;
					for (i2 = 0; i2 < n; i2++)
						v1[i1] += B[i1][i2] * v0[i2];
				}
				for (i1 = 0; i1 < n; i1++)
					v0[i1] = v1[i1];
				ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
			}
		}
							// 収束しない場合
		else {
			for (i1 = 0; i1 < n; i1++)
				v1[i1] = v2[i1];
			l1 = l2;
			k++;
		}
	}

	return ind;
}