ここではある３次実対称行列の対角化を行う。

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

#define N 3                     //実対象行列の次数
#define PI 3.1415926535897932385                //円周率

/*実対称行列とそれを対角化する直行行列を
格納する配列を受け取り対角化を行う関数。*/
void taikakuka(double a[N][N],double U[N][N]);  

int main()
{
        int i,j;
//配列aには実対称行列を入れる。また結果の直行行列を入れる配列Uを作っておく。
        double U[N][N],a[N][N]={{0.0,1.0,-1.0},{1.0,0.0,1.0},{-1.0,1.0,0.0}};

        taikakuka(a,U);         //対角化関数の呼び出し

//計算結果の対角行列の表示
        printf("対角行列は\n");
        for(i=0;i<N;i++) {
                for(j=0;j<N;j++) {
                        printf("%f ",a[i][j]);
                }
                printf("\n");
        }
//計算結果の直交行列の表示
        printf("直行行列は\n");
        for(i=0;i<N;i++) {
                for(j=0;j<N;j++) {
                        printf("%f ",U[i][j]);
                }
                printf("\n");
        }

        return 0;
}

void taikakuka(double a[N][N],double U[N][N])
{
        int i,j,l,m,p,q,count;
        double max,theta;
        double u[N][N],oldU[N][N],newa[N][N];

        FILE *output;
        output=fopen("output.dat","w");

/*計算結果としてだされた直行行列を
格納するための配列を単位行列に初期化しておく。*/
        for(p=0;p<N;p++) {
                for(q=0;q<N;q++) {
                        if(p==q) U[p][q]=1.0;
                        else U[p][q]=0.0;
                }
        }

        for(count=0;count<=10000;count++) {

/*配列olduは新たな対角化計算を行う前に
かけてきた直行行列を保持する。*/
                for(p=0;p<N;p++) {
                        for(q=0;q<N;q++) {
                                oldU[p][q]=U[p][q];
/*対角化するときの個々の直交行列を入れる配列uの初期化。(単位行列にする。)*/
                                if(p==q) u[p][q]=1.0;
                                else u[p][q]=0.0;
                        }
                }

/*非対角要素の中から絶対値の最大のものを見つける*/
                max=0.0;
                for(p=0;p<N;p++) {
                        for(q=0;q<N;q++) {
                                if(max<fabs(a[p][q]) && p!=q) {
                                        max=fabs(a[p][q]);
/*その最大のものの成分の行と列にあたる数を記憶しておく。*/
                                        i=p;
                                        j=q;
                                }
                        }
                }
/*先ほど選んだ最大のものが指定の値より小さければ対角化終了*/
                if(fabs(a[i][j]) < 1.0e-10) {
                        break;
                }
/*条件によってシータの値を決める*/
                if(a[i][i]==a[j][j]) theta=PI/4.0;
                else theta=atan(-2*a[i][j]/(a[i][i]-a[j][j]))/2.0;

/*ここでこのときに実対称行列にかける個々の直行行列uが決まるが
特にここでの計算の意味はない。(する必要はない。)*/
                u[i][j]=sin(theta);
                u[j][i]=-sin(theta);
                u[i][i]=u[j][j]=cos(theta);

/*ここでいままで実対称行列にかけてきた直行行列を配列Uに入れる。*/
                for(p=0;p<N;p++) {
                        U[p][i]=oldU[p][i]*cos(theta)-oldU[p][j]*sin(theta);
                        U[p][j]=oldU[p][i]*sin(theta)+oldU[p][j]*cos(theta);
                }

//対角化計算によってでた新たな実対称行列の成分を配列newaに入れる。
                newa[i][i]=a[i][i]*cos(theta)*cos(theta)
                        +a[j][j]*sin(theta)*sin(theta)-2.0*a[i][j]*sin(theta)*cos(theta);
                newa[j][j]=a[i][i]*sin(theta)*sin(theta)
                        +a[j][j]*cos(theta)*cos(theta)+2.0*a[i][j]*sin(theta)*cos(theta);
                newa[i][j]=newa[j][i]=0.0;
                for(l=0;l<N;l++) {
                        if(l!=i && l!=j) {
                                newa[i][l]=a[i][l]*cos(theta)-a[j][l]*sin(theta);
                                newa[l][i]=newa[i][l];
                                newa[j][l]=a[i][l]*sin(theta)+a[j][l]*cos(theta);
                                newa[l][j]=newa[j][l];
                        }
                }
                for(l=0;l<N;l++) {
                        for(m=0;m<N;m++) {
                                if(l!=i && l!=j && m!=i && m!=j) newa[l][m]=a[l][m];
                        }
                }

//次の対角化計算を行う行列の成分を配列aへ上書きする。
                for(p=0;p<N;p++) {
                        for(q=0;q<N;q++) {
                                a[p][q]=newa[p][q];
                        }
                }

                if(count==10000) {
                        printf("対角化するためにはまだ作業を繰り返す必要があります.\n");
                }

//対角化の途中計算をファイルに記録する。
                fprintf(output,"選択した非対角要素は %d 行 %d 列です。\n",i+1,j+1);
                fprintf(output,"このときのシータの値は %f です。\n",theta);
                fprintf(output,"また計算後の実対称行列は\n");
                for(i=0;i<N;i++) {
                        for(j=0;j<N;j++) {
                                fprintf(output,"%f ",newa[i][j]);
                        }
                        fprintf(output,"\n");
                }
                fprintf(output,"\nこのときまでに計算してきた直行行列は\n");
                for(i=0;i<N;i++) {
                        for(j=0;j<N;j++) {
                                fprintf(output,"%f ",U[i][j]);
                        }
                        fprintf(output,"\n");
                }
                fprintf(output,"\n\n");
        }
        
        fclose(output);
}

計算結果は

対角行列は
-2.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
直行行列は
0.577350 0.707107 -0.408248
-0.577350 0.707107 0.408248
0.577350 0.000000 0.816497
Press any key to continue

である。よって固有値は縮退ありの -2 と 1 である。