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

/*--------------------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;
int M[les][les],Mc[les][les];
double root[les];
int count,d;
FILE *in;

int Z(int type);
void change(double *a,double *b);
double Gauss(int mov, int *d);
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",&Mc[_i][_j]);
  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===========*/
void change(int *a,int *b)
{
 double tmp;
 tmp=*a;
 *a=*b;
 *b=tmp;
}

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()
{
 int T[les][les];
 int d=0;
 for(int i=0;i<count;i++)
  for(int j=count;j<count+count;j++)
   if(i==j-count)
    M[i][j]=1;
   else
    M[i][j]=0;
 if(!Gauss(count-1,&d))
  d=1;
 else
 {
  for(int quan=0;quan<count;quan++)
  {
   T[count-1][quan]=M[count-1][count+quan]/M[count-1][count-1];
   for(int k=1;k<count;k++)
   {
    double curr=M[count-k-1][count+quan];
    for(int i=count-1;i>=count-k;i--)
     curr-=M[count-k-1][i]*T[i][quan];
    T[count-k-1][quan]=curr;
   }
  }
  for(int i=0;i<count;i++)
   for(int j=0;j<count;M[i][j]=T[i][j++]);
 }
return d;
}


int Z(int type)
{
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(int j=0;j<count;j++)
  {
   M[count][j]=0;
   for(int i=0;i<count;M[count][j]+=M[i++][j]);
  }
  if(BackMatrix()==1)
   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(int 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;
}