高斯消元
July 11, 2013
Math
输入线性方程组的维数,然后随即生成一定有解的线性方程组的增广矩阵,求出解,然后输出时间和方程组的解,以及和标准解的误差(其实就是和标准解的方差)。
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <cstring>
#include <cmath>
#include <ctime>
using namespace std;
const double eps=1e-9;
const int MAX=20000;
double ans[MAX];
int n;
double **inputdata; double *result;
/*
输出标准解
*/
void printresult() {
for(int i=0;i<n;i++)
printf("%.4lf ",result[i]);
printf("\n");
}
/*
打印出生成的的增广矩阵
*/
void printb() {
for(int i=0;i<n;i++)
printf("%.4lf ",inputdata[i][n]);
printf("\n");
}
/*
随机生成增广矩阵
*/
void pro_inputdata()
{
int i,j,k;
double res=0;
printf("please input the number of the element:\n");
scanf("%d",&n);
srand(time(0));
result=(double*)malloc(n*sizeof(double));
inputdata=(double**)malloc(n*sizeof(double));
for(i=0;i<n;i++)
result[i]=(double)(rand()%10000)/1000.0;
for(i=0;i<n;i++) {
inputdata[i]=(double*)malloc((n+1)*sizeof(double));
res=0;
for(j=0;j<n;j++) {
k=(rand()%2==0)?1:(-1);
inputdata[i][j]=k*(double)(rand()%10000)/1000.0;
res=res+inputdata[i][j]*result[j];
}
inputdata[i][n]=res;
}
}
/*
高斯消元函数
*/
inline void solve(double **a, double ans[], int n) {
for (int i = 0; i < n; ++i) {
for (int j = i; j < n; ++j)
if (fabs(a[j][i]) > eps) {
for (int k = i; k <= n; ++k) swap(a[j][k], a[i][k]);
break;
}
for (int j = i+1; j < n; ++j) {
if (fabs(a[j][i]) > eps) {
double tmp = a[j][i]/a[i][i];
for (int k = i; k <= n; ++k) a[j][k] -= tmp * a[i][k];
}
}
}
for (int i = n - 1; i >= 0; --i) {
double sum = 0.0;
for (int j = n - 1; j > i; --j)
sum += ans[j] * a[i][j];
ans[i] = (a[i][n] - sum) / a[i][i];
}
for (int i = 0; i < n; ++i) printf("%.4lf ", ans[i]);
printf("\n");
return;
}
int main(void) {
freopen("in.txt", "r", stdin);
clock_t begin, end;
begin = clock();
pro_inputdata();
solve(inputdata, ans, n);
end = clock();
printf("Time is: %lf seconds\n",(double)(end-begin)/CLOCKS_PER_SEC);
double sum=0;
for (int i = 0; i < n; ++i) {
sum += (ans[i]-result[i])*(ans[i]-result[i]);
}
printf("Error is %.4lf\n",sqrt(sum));
return 0;
}
囧,等下在贴代码。== 刚才系统剪切板里面的内容不能粘贴到浏览器里面了,注销了一下就可以了=_=