选主元三角分解Gauss消去法的C++程序,急!!

要计算全主元三角分解Gauss消去法和列主元三角分解的Gauss消去法的C++程序,分开来写,谢谢啦

第1个回答  2009-12-30
全主元C语言程序,在这基础上修改一下就可以成为C++程序。列主元的自己修改一下也可以出来。你参考一下吧

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

//file <RGAUSS.C>全选主元高斯消去法

int rgauss(n,a,b)
//函数返回整型标志值。若返回0,表示原方程组的系数矩阵奇异;不为0,则表示正常返回
int n; //方程组的阶数
double a[],b[]; //a[]存放方程组的系数矩阵,返回时被破坏;
//b[]存放方程组右端常数向量,返回方程组的解向量;
{
int *js,l,k,i,j,is,p,q;
double d,t;
js=malloc(n*sizeof(int)); //开辟用于记忆列交换信息的动态空间
l=1; //置非奇异标志
for (k=0;k<=n-2;k++)
{ d=0.0;
for (i=k;i<=n-1;i++) //全选主元
for (j=k;j<=n-1;j++)
{ t=fabs(a[i*n+j]);
if (t>d) { d=t; js[k]=j; is=i;} //记忆行、列交换信息
}
if (d+1.0==1.0) l=0; //置奇异标志
else
{ if (js[k]!=k)
for (i=0;i<=n-1;i++) //列交换
{ p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if (is!=k)
{ for (j=k;j<=n-1;j++) //行交换
{ p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if (l==0) //奇异返回
{ free(js); printf("fail\n");
return(0);
}
d=a[k*n+k];
for (j=k+1;j<=n-1;j++) //归一化
{ p=k*n+j; a[p]=a[p]/d;}
b[k]=b[k]/d;
for (i=k+1;i<=n-1;i++) //消元
{ for (j=k+1;j<=n-1;j++)
{ p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if (fabs(d)+1.0==1.0) //奇异返回
{ free(js); printf("fail\n");
return(0);
}
b[n-1]=b[n-1]/d; //计算解向量的最后一个分量
for (i=n-2;i>=0;i--) //回代
{ t=0.0;
for (j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for (k=n-1;k>=0;k--) //恢复解向量
if (js[k]!=k)
{ t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;} //解向量行交换
free(js); //释放动态空间
return(1); //返回正常标志
}

main()
{

int i;
static double a[2][2]=
{ {5,1},{1,-1} };

static double c[4][4]=
{ {1,3,2,13},
{7,2,1,-2},
{9,15,3,-2},
{-2,-2,11,5} };

static double b[2]={200,1000};
static double d[4]={0,4,7,-1};

// static double M[3][3] = {};
// static double f[3] = {1,2,3} ;
if (rgauss(4,c,d)!=0)
for (i=0;i< 4;i++)
printf("x(%d1)=%e\n",i+1,d[i]);
printf("\n");

/* if (rgauss(4,M,f)!=0)
{
for (i=0;i<19;i++)
printf("x(%d)=%8.6f\n",i+1,f[i]);
}
*/}
相似回答