高斯-塞德尔方法解方程组

#include <iostream>
#include <cmath>
using namespace std;
int main()
{
    double a[3][3]= {{9,-2,1},{1,-8,1},{2,-1,-8}};//系数矩阵
    double b[3]= {6,-8,9};//方程组结果矩阵
    double chg[3][3];//变换后的系数矩阵
    double emax=0.0001;//精确度,根据要求更改


    //变换系数矩阵，用以初始化chg[][].
    for(int i=0; i<3; i++)
    {
        for(int j=0; j<2; j++)
        {

            if(j>=i)
                chg[i][j]=-a[i][j+1]/a[i][i];
            else
                chg[i][j]=-a[i][j]/a[i][i];

        }
        chg[i][2]=b[i]/a[i][i];
    }

    //为了节省空间，将迭代结果存入方程组结果矩阵b[]中
    double tem[3];//临时矩阵用来计算迭代的结果b[].
    double comp[3];//存储b[i-1],用来与b[]做差达到要求精度

    for(int i=0; i<3; i++)//初始化b[]
        b[i]=0;
    int se;
    for( se=0; se<1000; se++)//开始迭代求结果，为了防止出现不收敛的情况，设迭代次数上限1000次
    {
        for(int i=0; i<3; i++)//初始化tem[]和comp[]
        {
            tem[i]=0;
            comp[i]=b[i];
        }

        //以下二重循环用来求迭代一次的结果存入b[]
        for(int i=0; i<3; i++)
        {
            for(int j=0; j<2; j++)
            {
                if(i>j)
                    tem[i]+=chg[i][j]*b[j];
                else     //j>=i
                    tem[i]+=chg[i][j]*b[j+1];
            }
            tem[i]+=chg[i][2];
            b[i]=tem[i];

        }

        int cal=0;
        for(int i=0; i<3; i++)//比较是否达到要求精度emax
        {
            if(fabs(comp[i]-b[i])<emax)
                cal++;
        }

        if(cal>=3)//即b[3]中每个值都满足精度要求
            break;

    }
    if(se==1000) cout<<"迭代1000次还未达到指定精度要求，可能是不收敛"<<endl<<"此时的结果是：";
    for(int i=0; i<3; i++)//输出显示
        cout<<"x["<<i+1<<"]="<<b[i]<<endl;

    return 0;
}