Powered By Blogger

Wednesday 27 April 2011

Gauss Seidel method using C


#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
int dominat(float**);
float compute(float**,float,float,int);
int main(void)
{
float **arr,*x,*y,*z,e=0.001;
int i=0,j=0,dom,c=0;

//Dynamic allocation
x=(float*)malloc(1000*sizeof(float));
y=(float*)malloc(1000*sizeof(float));
z=(float*)malloc(1000*sizeof(float));
arr=(float**)malloc(3*sizeof(float*));
for(i=0;i<3;i++)
arr[i]=(float*)malloc(4*sizeof(float));
//Allocation complete


//User input
printf("\n***** JACOBI METHOD ****");
printf("\n\nEnter the coefficient matrix row wise: \n");
for(i=0;i<3;i++)
for(j=0;j<3;j++)
scanf("%f",&arr[i][j]);

printf("\n\nEnter the right hand side vector: \n");
for(i=0;i<3;i++)
scanf("%f",&arr[i][j]);
//Matrix updated with user input values


//Check whether given matrix is dominat or not
dom=dominat(arr);
if(dom==0)
{
printf("\nThe coefficient matrix is not in dominat form... Program will terminate !!!\n");
getch();
exit(1);
}
//checking complete

//Updation of matrix
printf("\nThe updated matris is:: \n\n");
for(i=0;i<3;i++)
{
for(j=0;j<4;j++)
{
if(i!=j&&j!=3)
arr[i][j]=(-arr[i][j]/arr[i][i]);
if(j==3)
arr[i][j]=(arr[i][j]/arr[i][i]); //printf("%f ",arr[i][j]);
} //printf("\n");
}
//End of updation

x[0]=0;y[0]=0;z[0]=0;i=0;
//X, Y, Z, i initialized

//Calculate x,y,z by using previous values
do
{
i=i+1;
x[i]=compute(arr,y[i-1],z[i-1],c);
y[i]=compute(arr,x[i],z[i-1],c+1);
z[i]=compute(arr,x[i],y[i],c+2);
}while((abs(x[i]-x[i-1])>=e)||(abs(y[i]-y[i-1])>=e)||(abs(y[i]-y[i-1])>=e));
//End of calculation

//Print result
printf("\nThe required solution is:: x-> %f,   y-> %f,   z-> %f\n",x[i],y[i],z[i]);
getch();
return 0;
}


int dominat(float **matrix)
{
int r=0;
if(abs(matrix[0][0])>=abs(matrix[0][1]+matrix[0][2])&&abs(matrix[1][1])>=abs(matrix[1][0]+matrix[1][2])&&abs(matrix[2][2])>=abs(matrix[2][0]+matrix[2][1]))
{
r++;
}
return r;
}


float compute(float **matrix,float a,float b,int c)
{
float result=0;
if(c==0)
result=((matrix[0][1]*a)+(matrix[0][2]*b)+matrix[0][3]);
if(c==1)
result=((matrix[1][0]*a)+(matrix[1][2]*b)+matrix[1][3]);
if(c==2)
result=((matrix[2][0]*a)+(matrix[2][1]*b)+matrix[2][3]);
return result;
}

No comments:

Post a Comment