jacobi

[b7505002 (...)]

// Jacobi method. //
#include<stdio.h>
#include<process.h>
#include<math.h>
#define dim 10
#define mistake 1.0e-5
double a[dim+1][dim+1], b[dim+1], old_x[dim+1], new_x[dim+1];
void mprintf(int m, int n, double a[][dim+1], double b[dim+1]);
void exchange(int n, int r1, int r2, double a[][dim+1], double b[dim+1]);
int main(void)
{
 int m,i,j,k,ipvt_store,flag,index;
 double temp;
 char ch;
 FILE *fp;

 fp=fopen("jacobi.dat","r");
 if ( fp==NULL )
 {
  printf("Can't open file.\n");
  exit(1);
 }

 // Loading... //
 fscanf(fp,"%d",&m);
 for (i=1; i<=m; i++)
 {
  for (j=1; j<=m; j++)
   fscanf(fp,"%lf ",&a[i][j]);
  fscanf(fp,"%lf\n",&b[i]);
 }
 for (i=1; i<=m; i++)
  fscanf(fp,"%lf ",&old_x[i]);

 // Pivoting... //
 for (k=1; k<m; k++)
 {
  ipvt_store=k;
  temp=a[k][k];
  for (i=k; i<=m; i++)
   if ( a[i+1][k] > temp )
   {
    ipvt_store=i+1;
    temp=a[i+1][k];
   }
  if ( ipvt_store != k )
   exchange(m, ipvt_store, k, a, b);
 }

 for (i=1; i<=m; i++)
 {
  temp=a[i][i];
  for (j=1; j<=m; j++)
   a[i][j]/=temp;
  a[i][i]=temp;
  b[i]/=temp;
 }

 index=0;
 do
 {
  index++;
  flag=0;
  for (i=1; i<=m; i++)
   new_x[i]=b[i];

  for (i=1; i<=m; i++)
  {
   for (j=1; j<=m; j++)
   {
    if ( i==j )
     continue;
    new_x[i] = new_x[i] - a[i][j]*old_x[j];
   }

   if ( fabs(new_x[i]-old_x[i]) > mistake )
    flag=1;
  }
  for (j=1; j<=m; j++)
   old_x[j]=new_x[j];

  printf("#%d\n",index);
  mprintf(m,m,a,new_x);
 }while ( flag==1 );

 for (i=1; i<=m; i++)
  printf("x%d = %10.5lf\n",i,new_x[i]);
 return 0;
}

void mprintf(int m, int n, double a[][dim+1], double b[dim+1])
{
 int i,j;
 printf("----------\n");
 for (i=1; i<=m; i++)
 {
  for (j=1; j<=n; j++)
   printf("%10.5lf ",a[i][j]);
  printf("%10.5lf\n",b[i]);
 }
}

void exchange(int n, int r1, int r2, double a[][dim+1], double b[dim+1])
{
 int j;
 double temp;
 for (j=1; j<=n; j++)
 {
  temp=a[r1][j];
  a[r1][j]=a[r2][j];
  a[r2][j]=temp;
 }
 temp=b[r1];
 b[r1]=b[r2];
 b[r2]=temp;
}

--
※ 發信站: 批踢踢實業坊(ptt.twbbs.org) 
◆ From: t195-250.dialup.seed.net.tw