#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
int main()
{
char filename1[256],filename2[256];
cout<<"Enter the path to input file:";
cin>>filename1;
cout<<"Enter the path to output file:";
cin>>filename2;
ifstream infile(filename1,ios::in|ios::nocreate);
ofstream outfile(filename2);
if(!infile)
{
cout<<"Can not open input file:"<<filename1<<’
’;
exit(1);
}
int n,w,k,i,j;
infile>>n;
w=n+1;
double *a,*x,*b,*g;
double t;
a=new double[n*n];
g=new double[w];
x=new double[w];
b=new double[w];
for(i=0;i<n*n;i++)
infile>>a[i];
for(i=1;i<w;i++)
infile>>b[i];
if(a[0]<1e-8) //Step1
{
cout<<"Method failed!";
exit(2);
}
/*其實這裡可以定義：#define a(int u,int v) a[u*n-n+v-1]
這樣引用起來比較直觀!但是運算是一樣的.*/
g[1]=a[n]; //Step2
a[n]=g[1]/a[0];
a[n+1]-=g[1]*a[n];
if(a[n+1]<1e-8) //Step3
{
cout<<"Method failed!";
exit(3);
}
for(i=3;i<w;i++) //Step4
{
g[1]=a[i*n-n]; //Step5
for(j=2;j<i;j++) //Step6
{
for(t=0,k=1;k<j;k++)
t+=g[k]*a[j*n-n+k-1];
g[j]=a[i*n-n+j-1]-t;
}
a[i*n-n]=g[1]/a[0]; //Step7
for(j=2;j<i;j++) //Step8
a[i*n-n+j-1]=g[j]/a[j*n-n+j-1];
for(t=0,k=1;k<i;k++) //Step9
t+=g[k]*a[i*n-n+k-1];
a[i*n-n+i-1]-=t;
if(a[i*n-n+i-1]<1e-8) //Step10
{
cout<<"Method failed!";
exit(4);
}
}
for(i=2;i<w;i++) //Step11
{
for(t=0,k=1;k<i;k++)
t+=a[i*n-n+k-1]*b[k];
b[i]-=t;
}
x[n]=b[n]/a[n*n-1]; //Step12
for(i=n-1;i>0;i--) //Step13
{
for(t=0,k=i+1;k<w;k++)
t+=a[k*n-n+i-1]*x[k];
x[i]=b[i]/a[i*n-n+i-1]-t;
}
for(i=1;i<w;i++)
outfile<<"g("<<i<<")="<<g[i]<<’
’;
outfile<<"y=(";
for(i=1;i<n;i++)
outfile<<b[i]<<’,’;
outfile<<b[n]<<");
";
outfile<<"diag(D)=(";
for(i=1;i<n;i++)
outfile<<a[i*n-n+i-1]<<’,’;
outfile<<a[n*n-1]<<");
";
outfile<<"L=
";
for(i=1;i<w;i++)
{ for(j=1;j<w;j++)
{ if(j<=i)
outfile<<’ ’<<a[i*n-n+j-1];
else outfile<<’ ’<<’0’;
}
outfile<<’
’;
}
outfile<<"The roots are:("; //Step14
for(i=1;i<n;i++)
outfile<<x[i]<<’,’;
outfile<<x[n]<<");";
delete a;
delete g;
delete b;
delete x;
infile.close();
outfile.close();
return 0;
}
input=============
10
1 1 1 1 1 1 1 1 1 1
1 2 2 2 2 2 2 2 2 2
1 2 3 3 3 3 3 3 3 3
1 2 3 4 4 4 4 4 4 4
1 2 3 4 5 5 5 5 5 5
1 2 3 4 5 6 6 6 6 6
1 2 3 4 5 6 7 7 7 7
1 2 3 4 5 6 7 8 8 8
1 2 3 4 5 6 7 8 9 9
1 2 3 4 5 6 7 8 9 10
8 2 -12 4 -45 6 -5 9 10 30
output=============
g(1)=1
g(2)=1
g(3)=1
g(4)=1
g(5)=1
g(6)=1
g(7)=1
g(8)=1
g(9)=1
g(10)=-6.27744e+066
y=(8,-6,-14,16,-49,51,-11,14,1,20);
diag(D)=(1,1,1,1,1,1,1,1,1,1);
L=
1 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0
1 1 1 0 0 0 0 0 0 0
1 1 1 1 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0
1 1 1 1 1 1 0 0 0 0
1 1 1 1 1 1 1 0 0 0
1 1 1 1 1 1 1 1 0 0
1 1 1 1 1 1 1 1 1 0
1 1 1 1 1 1 1 1 1 1
The roots are:(14,8,-30,65,-100,62,-25,13,-19,20);
