Househoder Method and Jacobi Rotaion

/* C code to tridiagonalize a Matrix using Householder's Method and then finding the eigen values using Jacobi Rotation Method*/

#include<stdio.h>
#include<math.h>
#define ERROR 0.0000001
//Function to read input matrix
void read_matrix(int order ,float matrix[][order])
{
    int i, j;
    for (i=0; i<order;i++)
        for (j=0;j<order;j++)
            scanf("%f",&matrix[i][j]);
}

//Function to print the matrix
void print_matrix(int order , float matrix[][order])
{
    int i, j;
    for (i=0; i<order;i++)
    {
        for (j=0;j<order;j++)
            printf("%.4f\t",matrix[i][j]);
        printf("\n");
    }
}

int check_symmetric(int n , float a[][n])
{
    int i, j, tag=1;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            if (a[i][j]!=a[j][i])
            {
                tag=0;
                break;
            }   
        }
       
        if (tag==0)
            break;
    }

    return tag;
}

//Function to multiply two matrices "a" and "b" and put the result in matrix "c"
float multiply(int n , float a[][n], float b[][n] , float c[][n] )
{
    int i,j,k ;
    float sum;

    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
        {
            sum=0;
            for(k=0;k<n;k++)
                sum=sum+(a[i][k])*(b[k][j]);
            c[i][j]=sum;
        }
}   

/*Function to get an identity matrix*/
void identity(int n , float S[][n])
{
    int i , j;
    for (i=0; i<n;i++)
        for (j=0;j<n;j++)
        {
            if (i==j)
                S[i][i]=1;
            else
                S[i][j]=0;   
        }
}

/*Function to copy one matrix into another matrix*/
void copy(int n, float a[][n] , float b[][n])
{
    int i , j;
    for (i=0;i<n;i++)
        for(j=0;j<n;j++)
            a[i][j]=b[i][j];
}

void transpose_vector(int n, float vector[] , float matrix[][n])
{
    int i,j;
    float sum=0;
   
    for(i=0;i<n;i++)
    {
        for (j=0;j<n;j++)
        {
            matrix[i][j] = vector[i]*vector[j];
        }
    }
}

void print_vector(int n, float a[])
{
    int i;
    for (i=0;i<n;i++)
    {
        printf("%12.4f  ",a[i]);
    }
    printf("\n");   
}

void swap(float a[], int i, int j)
{
    float temp;
    temp=a[i];
    a[i]=a[j];
    a[j]=temp;
}

void bubblesort( int n, float a[])
{
    int i , j;
     for(j=n-2;j>=1;j--)
         for(i=0;i<=j;i++)
             if(a[i]>a[i+1])
                 swap(a,i,i+1);        
}

void jacobi_rotation(int n , float a[][n])
{
    int i,j;
    float pi = 3.1415926;
    float S[n][n] , invS[n][n] ,temp1[n][n] ,temp2[n][n], theta;   

    identity(n,S);
    identity(n,invS);

    int u,v,tag=1;
    while(1)
    {
        tag=1;
        for (i=0;i<n;i++)
        {
           
            for (j=i+1 ; j <n;j++)
            {
                if ((a[i][j])!=0)
                {
                    if (a[i][i]!=a[j][j])
                        theta = (atan(2*a[i][j]/(a[i][i]-a[j][j])))/2;

                    else
                    {
                        if (a[i][j]>0)
                            theta = pi/4;
                       
                        else if (a[i][j]<0)
                            theta = -pi/4;
                    }
                                           
                   
                    identity(n,S);
                    identity(n,invS);

                    S[i][i]= invS[i][i]= invS[j][j]=S[j][j]=cos(theta);
                    S[i][j]= invS[j][i]= -sin(theta);
                    S[j][i]= invS[i][j]=  sin(theta);

                    multiply(n,invS,a,temp1);
                    multiply(n,temp1,S,temp2);
   
                    copy(n,a,temp2);
                   
                }
                else continue;
            }   
        }   
   
        for (u=0;u<n;u++)
        {
            for (v=u+1;v<n;v++)
            {
                if (abs(a[u][v]) > ERROR)
                {
                    tag=0;
                    break;
                }   
            }
            if (tag==0)
                break;
        }
   
        if (tag==1)
            break;
   
    }   
}


/*main function*/
int main()
{
    int n,i,j,sign,k;
    scanf("%d",&n);

    float a[n][n] , P[n][n] ,temp1[n][n] ,temp2[n][n] , I[n][n] , w[n], transpose[n][n], x , sum,s,temp;
    read_matrix(n,a);

    int check;
    check = check_symmetric(n,a);
    if (check==0)
    {
        printf("Given matrix is not symmetric\n");
        return 0;
    }
   
    for(i=0;i<n-2;i++)
    {
        j=0;
        sum=0;
        while(j<=i)
        {
            w[j++]=0;
        }
       
        for(k=i+1;k<n;k++)
            sum=sum+a[i][k]*a[i][k];   
       
        s= sqrt(sum);
        if (s==0)
            continue;
        //printf("s : %f\n",s);
        if (j==i+1)
        {
            if (a[i][j]<0)
                sign=-1;
            else
                sign=1;
               
            temp = (1+a[i][j]*sign/s)/2;
            x=sqrt(temp);
            w[j]=x;
            j++;
        }
       
        for (j=i+2;j<n;j++)
        {
            w[j]=(a[i][j]*sign)/(2*s*x);
        }
       
        identity(n, I);
        transpose_vector(n,w,transpose);
       
        for(k=0;k<n;k++)
        {
            for(j=0;j<n;j++)
            {
                P[k][j]=I[k][j]-2*transpose[k][j];
            }
        }

        //print_matrix(n, P);

        multiply(n,a,P,temp1);
       
        multiply(n,P,temp1,temp2);       
        copy(n,a,temp2);
       
        //print_matrix(n, a);
    }
   
    //printf("\nTridiagonal matrix using Householder Method is : \n\n");
    //printf("\n");
    print_matrix(n,a);
   
    //printf("\nEigne values of Tridiagonal  : \n\n");
    jacobi_rotation(n, a);

    float eigen[n];
    for (i=0;i<n;i++)
        eigen[i]=a[i][i];
       
    bubblesort(n,eigen);   
    //printf("\n");
    for (i=n-1;i>=0;i--)
        printf("%.4f\n",eigen[i]);
   
    return 0;   
}