Tridiagonalization using Given's Method

#include<stdio.h>
#define max 1000
#include<math.h>

//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("%f\t",matrix[i][j]);
        printf("\n");
    }
}

//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];
}

/*main function*/
int main()
{
    int n,i,j;

    printf("Enter the order of matrix : ");
    scanf("%d",&n);
    float a[n][n] , S[n][n] , invS[n][n] ,temp1[n][n] ,temp2[n][n];   
    read_matrix(n,a);

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

    for (i=0;i<=n-3;i++)
    {
        for (j=i+2 ; j <= n-1;j++)
        {
            if ((a[i][j])!=0)
            {
                float theta = atan(a[i][j]/a[i][i+1]);

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

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

                multiply(n,invS,a,temp1);
                multiply(n,temp1,S,temp2);
   
                copy(n,a,temp2);
            }
            else continue;
        }   
    }   
    print_matrix(n,a);
    return 0;   
}