解线性方程组的迭代法--大型稀疏矩阵---雅可比迭代法C++实现

//
//    Jacobi  迭代求解线性方程组   //
//
#include<math.h>
#define//设置方程的最大维数
#define
#define//求解精度
#include<iostream>
using namespace std;

int main()
{
  int n;
  int i, j, k;
  double err;
  static double a[MAX_N][MAX_N], b[MAX_N][MAX_N], c[MAX_N], g[MAX_N];
  static double x[MAX_N], nx[MAX_N];
  printf("\nInput n value(dim of AX=C):");//输入方程维度n*n
  cin>>n; 
  if (n > MAX_N) {
    cout << "The input n is larger than Max_N,please redefine the Max_N." << endl;
    return 1;
  }
  else if (n <= 0) {
    printf("Please input a number between 1 and %d\n", MAX_N);
    return 1;
  }
  //初始化工作：
  //任务一：开始输入AX=C的系数矩阵A
  printf("Now input the matrix a(i,j),i,j=0,...,%d:\n", n - 1);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      cin >> a[i][j];
  //任务二：输入C矩阵
  printf("Now input the matrix b(i),i=0,...,%d\n", n - 1);
  for (i = 0; i < n; i++)  cin >> c[i];
  //任务三：生成  x^(k+1) = b * x^(k)+ g  迭代矩阵B 和 g [ x^k表示第k次获得的x]
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      if (a[i][i] == 0)
        cerr << "暂时不支持系数矩阵对角线元素为0" << endl;
      b[i][j] = -a[i][j] / a[i][i];
      g[i] = c[i] / a[i][i];      //先假设a[i][i]!=0,学习中暂时用不到
    }
    b[i][i] = 0;
  }
  for (i = 0; i < MAXREPT; i++) {//最大迭代周期数

    for (j = 0; j < n; j++)
      nx[j] = g[j];
    for (j = 0; j < n; j++) {
      for (k = 0; k < n; k++) {
        if (j == k) continue;
        nx[j] += b[j][k] * x[k];    // x^(k+1) = b * x^(k)+ g  迭代
      }
    }

    err = 0;      //x^(k+1)-x^k的误差
    for (j = 0; j < n; j++) {
      if (err < fabs(nx[j] - x[j])) {
        err = fabs(nx[j] - x[j]);    //求最大差即最大范式
      }
    }

    for (j = 0; j < n; j++)  x[j] = nx[j];

    if (err < epsilon) {    //控制误差
      cout << "Solve ... x_i=" << endl;
      for (j = 0; j < n; j++)
        cout << x[j] << endl;
      return 0;
    }
  }
  printf("After %d repeat ,no result\n", MAXREPT);
  return 1;
}

/*
方程组：
  64x1 - 3x2 - x3 =14
  2x1 -90 x2+x3 = -5
  x1 + x2 + 40x3 = 20
  输入样例：
  3
  64 -3 -1 2 -90 1 1 1 40
  14 -5 20
*/