计算方法之高斯消去

#include "stdio.h"
#include "math.h"
#include "stdlib.h"
typedef double arrtype;
//高斯消去法
void GuassScale(arrtype coff[],int n,arrtype *x = NULL,int m = 1)
{
int col = n + m;
arrtype *s = new arrtype
;
int *nrow = new int
;
arrtype *base = coff;
for(int i = 0; i<n;i++)
{
s[i] = 0;
for(int j = 0;j<n;j++)
{
arrtype t = fabs(coff[i*col + j]);
if(s[i] < t)
s[i] = t;
}
nrow[i] = i;
}
for(i=0;i<n-1;i++)
{
int p = i;
arrtype t = fabs(coff[nrow[p]*col+i])/s[nrow[p]];
for(int j= i+1;j<n;j++)
{
arrtype temp= fabs(coff[nrow[j]*col+i])/s[nrow[j]];
if(t<temp)
{
t=temp;
p=j;
}
}
if(coff[nrow[p]*col+i]==0)
return;
if(nrow[p]!=nrow[i])
{
int ncopy = nrow[i];
nrow[i] = nrow[p];
nrow[p] =ncopy;
}
for(j= i+1;j<n;j++)
{
arrtype m = coff[nrow[j]*col+i]/coff[nrow[i]*col+i];
for(int k=i+1;k<col;k++)
{
coff[nrow[j]*col+k]-=coff[nrow[i]*col+k]*m;
}
}
}
if(coff[nrow[n-1]*col+n-1]==0)
return;
for(i=n;i<col;i++)
{
x[(n-1)*m+(i-n)]=coff[nrow[n-1]*col+i]/coff[nrow[n-1]*col+n-1];
}
for(i=n-2;i>=0;i--)
{
for(int j = n;j<col;j++)
{
for(int k=n-1;k>i;k--)
coff[nrow[i]*col+j]-=coff[nrow[i]*col+k]*x[k*m+(j-n)];
x[i*m+(j-n)]=coff[nrow[i]*col+j]/coff[nrow[i]*col+i];
}
}
}
void main(int argc,char* argv[])
{
FILE *fpi,*fpo;
if(argc > 1)
fpi=fopen(argv[1],"r");
else
fpi=stdin;
if(!fpi)return;
if(argc>2)
fpo = fopen(argv[2],"w");
else
fpo = stdout;
if(!fpo)return;
int m,n;
printf("请输入ｎ和ｍ的值：/n");
fscanf(fpi,"%d",&n);
fscanf(fpi,"%d",&m);
arrtype *co;
co = new arrtype[n*(n+m)];
for(int i = 0;i<n*(n+m);i++)
{
float t;
fscanf(fpi,"%f",&t);
co[i]=t;
}
fclose(fpi);
arrtype *x = new arrtype[m*n];
GuassScale(co,n,x,m);
delete[] co;
for(i=0;i<n;i++)
{
fprintf(fpo,"x[%d]=",i);
for(int j=0;j<m;j++)
fprintf(fpo,"%f",x[i*m+j]);
fprintf(fpo,"/n");
}
delete[] x;
fclose(fpo);
}