高斯消元

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;
}

  囧,等下在贴代码。== 刚才系统剪切板里面的内容不能粘贴到浏览器里面了,注销了一下就可以了=_=

comments powered by Disqus