/*  ヤコビ法による行列固有値の計算　*/

/*  ---- 近似計算の途中を表示 ----  */

#include<stdio.h>

#include<math.h>

#define  N        4

#define  EPSILON  1.0E-4

double  a[N][N]={{ 8.00, -0.01,  1.05,  0.02 },

                {-0.01,  3.40,  0.24, -0.10 },

                { 1.05,  0.24,  3.50, -0.03 },

                { 0.02, -0.10, -0.03, -1.50 }};

int    k;

double  t[N][N],u[N][N],x[N][N];

void main()

{

    int i,j;

    void jacobi(void);

    printf("\n ================= given matrix ================= \n\n");

    for(i=0;i<N;i++)

    {

        for(j=0;j<N;j++)

            printf(" %8.4f",a[i][j]);

        printf("\n");

    }

    jacobi();

    printf("\n ======== characteristic value of matrix ======== \n\n");

    for(i=0;i<N;i++)

        printf(" %8.4f ",a[i][i]);

    printf("\n\n ======== characteristic vector of matrix ======= \n\n");

    for(j=0;j<N;j++)

    {

        printf("  u%d = (",j+1);

        for(i=0;i<N-1;i++)

            printf(" %8.4f,",u[i][j]);

        printf(" %8.4f )\n\n",u[N-1][j]);

    }

}

void jacobi()

{

    int   i,j,m,p,q;

    double xx,aa,amax,sn,cs;

    void  max_det(int *,int *,double *);

    void  calc(int,int,double *,double *);

    for(i=0;i<N;i++)

        for(j=0;j<N;j++)

            u[i][j] = 0.0;

    for(i=0;i<N;i++)

            u[i][i] = 1.0;

    max_det(&p,&q,&amax);

    printf("\n\n **** iteration for approximation \n\n");

    k = 0;

    do

    {

        calc(p,q,&sn,&cs);

        for(i=0;i<N;i++)

            for(j=0;j<N;j++)

                t[i][j] = 0.0;

        for(i=0;i<N;i++)

                t[i][i] = 1.0;

        t[p][q] = sn; t[q][p] = -sn; t[p][p] = t[q][q] = cs;

        for(i=0;i<N;i++)

            for(j=0;j<N;j++)

           {

                xx = 0.0;

                for(m=0;m<N;m++)

                      xx = xx+a[i][m]*t[m][j];

                x[i][j]= xx;

            }

        for(i=0;i<N;i++)

            for(j=0;j<N;j++)

            {

                aa = 0.0;

                for(m=0;m<N;m++)

                    aa = aa+t[m][i]*x[m][j];

                a[i][j] = aa;

            }

        for(i=0;i<N;i++)

            for(j=0;j<N;j++)

                x[i][j] = u[i][j];

        for(i=0;i<N;i++)

            for(j=0;j<N;j++)

            {

                xx = 0.0;

                for(m=0;m<N;m++)

                    xx = xx+x[i][m]*t[m][j];

                u[i][j] = xx;

            }

        k = k+1;

        printf("   #%d   ",k);

        for(i=0;i<N;i++)

        {

            for(j=0;j<N;j++)

                printf(" %8.5f ",a[i][j]);

            printf("\n        ");

        }

        printf("\n");

        max_det(&p,&q,&amax);

    }

    while(amax > EPSILON);

}

void calc(int p,int q,double *sn,double *cs)

{

    double t,tt,absn;

    if(a[p][p] != a[q][q])

    {

        t = 2.0*a[p][q]/(a[q][q]-a[p][p]);

        tt = sqrt(t*t+1.0);

        absn = sqrt(0.5*(tt-1.0)/tt);

        if(t>0.0)

            *sn =  absn;

        else

            *sn = -absn;

        *cs = sqrt(0.5*(tt+1.0)/tt);

    }

    else

        *sn = *cs = 0.5*sqrt(2.0);

}

void max_det(int *p,int *q,double *amax)

{

    int i,j;

    double absa;

    *amax = 0.0;

    for(i=0;i<N-1;i++)

        for(j=i+1

        ;j<N;j++)

        {

            absa = fabs(a[i][j]);

            if(absa > *amax)

            {

                *p = i; *q = j;

                *amax = absa;

            }

        }

}