написанные программы / программирование метода гаусса и гаусса-зейделя / Gauss_old / CZ2

/*--------------------Gauss 1.1
----------------g204s07------------*/
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>

/*============global var===============*/
const int les = 10;
const double E=1E-20;
double M[les][les],Mc[les][les];
double root[les];
int count,d;
FILE *in;

int Z(int type);
double Gauss(int mov, int *d);
void change(double *aa,double *bb);
int BackMatrix(void);
/*===========main===================*/
void main(void)
{
register int _i,_j,_k;
 /*-----------------------file-----*/
 if((in = fopen("D:\\gIn.txt","rb"))==NULL)
 {printf("Error on open file. \n");getch();exit(1);}

 fscanf(in,"%d",&count);

 for(_i=0;_i<count;_i++)
  for(_j=0;_j<=count;_j++) fscanf(in,"%d",&M[_i][_j]);

 rewind(in);
 fscanf(in,"%d",&count);

 for(_i=0;_i<count;_i++)
  for(_j=0;_j<=count;_j++) fscanf(in,"%d",&Mc[_i][_j]);



/*printf("Sourse Matrix: \n");
for(_i=0;_i<count;_i++){
	 for(_j=0;_j<count;_j++) {printf("%d ",Cx[_j]);}
	 printf("\n");
							  */

  switch((int)Gauss(0,&d))
	{
	 case 0:{
	  printf("error: the quation has no roots ");
	  break;
	 }
	 case -1:
	 {
	  printf("error: the quation has infinite set of roots ");
	  break;
	 }
	 default:
	 {
	  root[count-1]=M[count-1][count]/M[count-1][count-1];
	  for(_k=1;_k<count;_k++)
	  {
		float curr=M[count-_k-1][count];
		 for(_i=count-1;_i>=count-_k;_i--)
		  curr-=M[count-_k-1][_i]*root[_i];
		root[count-_k-1]=curr;
	  }
	 }
  }

  printf("Roots of Gauss: \n");
  for(_i=0;_i<count;_i++) {printf("%f",root[_i]);printf("\n");}

  if( Z(2)) printf("error: matrix is upborn \n");
  else
  {
	 printf("Roots of Gauss+Zeidel: \n");
	 for(_i=0;_i<count;_i++) {printf("%f",root[_i]);printf("\n");}
	 printf("Epsilon: %e",E);printf("\n");
  }

  getch();
  fclose(in);
}
/*===========================Gauss===========*/
double Gauss(int mov, int *d)
{
 *d=0;
 double det=1;
 for(int k=0;k<count-1;k++)
 {
  if(!M[k][k])
  {
	int i=k;
	for(i;(i<=count)&&(!M[i][k]);i++);
	if(i>count)
	{
	 det=0;
	 if(!mov)
	 {
	  int j=k+1;
	  for(j;(j<=count)&&(!M[k][j]);j++);
	  if(j>=count)
	  {
		if(j>count)
		 *d=1;
		return det;
	  }
	  for(int i=0;i<count;i++)
		change(&M[i][k],&M[i][j]);
	 }
	 else
	  return det;
	}
	else
	{
	 for(int j=0;j<=count+mov;j++)
	  change(&M[i][k],&M[i][j]);
	 if(mov==-1)
	  det*=-1;
	}
  }
  if(mov==-1)
	det*=M[k][k];
  for(int j=count+mov;j>k-1;j--)
	M[k][j]=M[k][j]/M[k][k];
  for(int i=k+1;i<count;i++)
	for( j=count+mov;j>k-1;j--)
	 M[i][j]=M[i][j]-M[k][j]*M[i][k];
 }
 switch(mov)
 {
  case -1:
  {
	det*=M[count-1][count-1];
	break;
  }
  case 0:
  {
	if(!M[count-1][count-1])
	 if(!M[count-1][count])
	 det=-1;
	else
	 det=0;
	break;
  }
  default:
  {
	if(!M[count-1][count])
	 det=0;
	break;
  }
 }
 return det;
}
/*================zeidel=================================*/
int BackMatrix(void)
{ register int i,j,quan,k;
 int T[les][les];
 int sk=0;
 for(i=0;i<count;i++)
  for(j=count;j<count+count;j++)
	if(i==j-count)
	 M[i][j]=1;
	else
	 M[i][j]=0;
 if(!Gauss(count-1,&d))
  sk=1;
 else
 {
  for(quan=0;quan<count;quan++)
  {
	T[count-1][quan]=M[count-1][count+quan]/M[count-1][count-1];
	for(k=1;k<count;k++)
	{
	 double curr=M[count-k-1][count+quan];
	 for(i=count-1;i>=count-k;i--)
	  curr-=M[count-k-1][i]*T[i][quan];
	 T[count-k-1][quan]=curr;
	}
  }
  for(i=0;i<count;i++)
	for(j=0;j<count;M[i][j]=T[i][j++]);
 }
return sk;
}
void change(int *aa,int *bb)
{
 double tmp;
 tmp=*aa;
 *aa=*bb;
 *bb=tmp;
}
int Z(int type)
{ int i,j;
double X[10][2],C[10][10],curr,Ecurr=.1,a,b,A,B;
 do
 {
  for(int i=0;i<count;i++)
	for(int j=0;j<count;M[i][j++]=Mc[j][i]);
  for(j=0;j<count;j++)
  {
	M[count][j]=0;
	for(int i=0;i<count;M[count][j]+=M[i++][j]);
  }
  if(BackMatrix())
	return 1;
  for( i=0;i<count;i++)
  {
	curr=1-M[count][i]*Ecurr;
	for(int j=0,mj=-1;j<count;(j++,mj++))
	 if(!j)
	 {
	  C[i][0]=0;
	  for(int k=0;k<count;k++)
		 C[i][0]+=(M[i][k]-Ecurr)*Mc[k][count];
	  C[i][0]/=curr;
	  }
	 else
	 {
	  if(i==j-1)
		mj++;
	  C[i][j]=-M[count][mj]*Ecurr/curr;
	 }
  }
  A=B=0;
  for(i=0;i<count;i++)
  {
	a=b=0;
	for(int j=1;j<count;j++)
	{
	 a+=fabs(C[i][j]);
	 b+=fabs(C[j-1][i+1]);
	}
	if(A<a)
	 A=a;
	if(B<b)
	 B=b;
  }
  Ecurr*=.1;
 }
 while(A>=1&&B>=1);
 for( i=0;i<count;i++)
 {
  X[i][0]=C[i][0];
  if(type==2)
	X[i][1]=0;
 }
 int pos=0;
 do
 {
  if(type!=2)
  {
	for(int i=0;i<count;i++)
	{
	 if(type)
	  curr=X[i][pos];
	 X[i][pos]=C[i][0];
	 for(int j=1,mj=0;j<count;(j++,mj++))
	 {
	  if(i==j-1)
		mj++;
	  if(!type)
		X[i][pos]-=X[mj][pos?0:1]*C[i][j];
	  else
		X[i][pos]-=X[mj][pos]*C[i][j];
	 }
	 if(!type)
	  curr=X[i][pos]-X[i][pos?0:1];
	 else
	  curr=X[i][pos]-curr;
	}
  }
  else
  {
	int num;
	curr=0;
	for(int i=0;i<count;i++)
	{
	 if(fabs(curr)<fabs(X[i][0]))
	 {
	  curr=X[i][0];
	  num=i;
	 }
	}
	X[num][1]+=curr;
	for(i=0;i<count;i++)
	{
	 if(i==num)
	  X[i][0]=0;
	 else
	  X[i][0]-=C[i][i<num?num:num+1]*curr;
	}
  }

 }
 while(fabs(curr)>E);

 if(type==2)
  pos=1;
 for( i=0;i<count;i++)
  root[i]=X[i][pos];
 return 0;
}