#include <stdio.h>
#include <math.h>
#define Nmax 10
void JACOBI(int n,double a[Nmax][Nmax],double ev[Nmax],double vec[Nmax][Nmax]);
/*--------------------------------------------------------------------*/

int main()
{
    /* n：固有方程式の次数(配列宣言はn) */
    /* ev[k]：固有値(大きい順に出力)    */
    /* vec[i=0〜n-1][k]：固有ベクトル   */

    int i,j,k;
    int n=3;
    double a[Nmax][Nmax];
    double ev[Nmax];
    double vec[Nmax][Nmax];

    /*検算用*/
    double ar[Nmax][Nmax];
    double c[Nmax][Nmax];

    a[0][0]=4; a[0][1]=2; a[0][2]=3;
    a[1][0]=2; a[1][1]=6; a[1][2]=1;
    a[2][0]=3; a[2][1]=1; a[2][2]=5;

    for(i=0;i<=n-1;i++){
        for(j=0;j<=n-1;j++){
            ar[i][j]=a[i][j];
        }
    }

    JACOBI(n,a,ev,vec);

    printf("***Eigenvalue***\n");
    printf("%16.6e  %16.6e  %16.6e\n",ev[0],ev[1],ev[2]);
    printf("***Eigenvector***\n");
    for(i=0;i<=n-1;i++){
        printf("%16.6e  %16.6e  %16.6e\n",vec[i][0],vec[i][1],vec[i][2]);
    }

    /* 出力結果                                      */
    /* *** Eigen value ***                           */
    /* 9.000000e+000   4.732051e+000   1.267949e+000 */
    /* *** Eigen vector ***                          */
    /* 5.773503e-001   2.113249e-001   7.886751e-001 */
    /* 5.773503e-001  -7.886751e-001  -2.113249e-001 */
    /* 5.773503e-001   5.773503e-001  -5.773503e-001 */

    for(i=0;i<=n-1;i++){
        for(j=0;j<=n-1;j++){
            c[i][j]=0.0;
            for(k=0;k<=n-1;k++){
                c[i][j]=c[i][j]+ar[i][k]*vec[k][j];
            }
        }
    }

    /*検算  [A]{x}=λ{x} */
    printf("***検算結果 [A]{x}***\n");
    for(i=0;i<=n-1;i++){
        printf("%16.6e  %16.6e  %16.6e\n",c[i][0],c[i][1],c[i][2]);
    }
    printf("***検算結果 λ{x}***\n");
    for(i=0;i<=n-1;i++){
        printf("%16.6e  %16.6e  %16.6e\n",ev[0]*vec[i][0],ev[1]*vec[i][1],ev[2]*vec[i][2]);
    }

    printf("リターンキーを押して下さい");
    getchar();
    return 0;

}
/*--------------------------------------------------------------------*/

void JACOBI(int n,double a[Nmax][Nmax],double ev[Nmax],double vec[Nmax][Nmax])
{
    /* Jacobi法による固有値・固有ベクトル計算（奥村先生）*/
    /* 固有値ev[k],固有ベクトルvec[i][k]                 */

    int i,j,k,im=0;
    double s,c,aoff,eps,t,aa=0.0,aj,ak,vj,vk,vmax,bb;

    for(i=0;i<=n-1;i++){
        for(j=0;j<=n-1;j++){
            vec[i][j]=0.0;
        }
        vec[i][i]=1.0;
    }
    aoff=0.0;
    for(j=0;j<=n-2;j++){
        for(k=j+1;k<=n-1;k++){
            aoff=aoff+a[j][k]*a[j][k];
        }
    }
    eps=aoff*2.0;
    for(j=0;j<=n-1;j++){
        eps=eps+a[j][j]*a[j][j];
    }
    eps=eps*0.00000000005;

    while(aoff>eps){
        for(j=0;j<=n-2;j++){
            for(k=j+1;k<=n-1;k++){
                if(a[j][k]==0.0){
                    t=0.0;
                }else{
                    aa=(a[k][k]-a[j][j])/(2.0*a[j][k]);
                }
                if(aa>=0){
                    t=1.0/(aa+sqrt(1.0+aa*aa));
                }else{
                    t=1.0/(aa-sqrt(1.0+aa*aa));
                }
                c=1.0/sqrt(1.0+t*t);
                s=t*c;
                a[j][j]=a[j][j]-t*a[j][k];
                a[k][k]=a[k][k]+t*a[j][k];
                a[j][k]=0.0;
                for(i=0;i<=j-1;i++){
                    aj=a[i][j];
                    ak=a[i][k];
                    a[i][j]=aj*c-ak*s;
                    a[i][k]=aj*s+ak*c;
                }
                for(i=j+1;i<=k-1;i++){
                    aj=a[j][i];
                    ak=a[i][k];
                    a[j][i]=aj*c-ak*s;
                    a[i][k]=aj*s+ak*c;
                }
                for(i=k+1;i<=n-1;i++){
                    aj=a[j][i];
                    ak=a[k][i];
                    a[j][i]=aj*c-ak*s;
                    a[k][i]=aj*s+ak*c;
                }
                for(i=0;i<=n-1;i++){
                    vj=vec[i][j];
                    vk=vec[i][k];
                    vec[i][j]=vj*c-vk*s;
                    vec[i][k]=vj*s+vk*c;
                }
            }
        }
        aoff=0.0;
        for(j=0;j<=n-2;j++){
            for(k=j+1;k<=n-1;k++){
                aoff=aoff+a[j][k]*a[j][k];
            }
        }
    }
            
    for(k=0;k<=n-1;k++){
        ev[k]=a[k][k];
    }
    
    /*ソート（大きい順）*/
    for(k=0;k<=n-2;k++){
        vmax=-1.0E+30;
        for(i=k;i<=n-1;i++){
            if(ev[i]>=vmax){
                im=i;
                vmax=ev[i];
            }
        }
        aa=ev[k];ev[k]=ev[im];ev[im]=aa;
        for (i=0;i<=n-1;i++){
            bb=vec[i][k];
            vec[i][k]=vec[i][im];
            vec[i][im]=bb;
        }
    }

}
/*--------------------------------------------------------------------*/