(((silence)))

NA::VecDP jacobiIteration ( NA::MatDP const& a, // 系数矩阵
                            NA::VecDP const& b, // 右端向量
                            double    const  eps,       // 精度
                            int       const  M )        // 最多迭代次数
{
    NA::VecDP y(0.0, b.size());
    NA::VecDP x(0.0, b.size());
    double  tmp, maxXY;
    for ( int m=0; m<M ; ++m ) {
        maxXY = 0;
        for ( int i=0; i<b.size(); ++i ) {
            tmp = 0;
            for ( int j=0; j<b.size(); ++j ) {
                tmp+=x[j]*a[i][j];
            }
            tmp-=x[i]*a[i][i];
            y[i] = (b[i]-tmp)/a[i][i];

            tmp = fabs(x[i]-y[i]);
            if (tmp>maxXY)
                maxXY=tmp;
        }

        if (maxXY<eps)
            return y;
        x=y;
    }
    return NA::VecDP(0);
}

// 高斯-赛德尔迭代法
NA::VecDP gaussSeidelIteration (
                            NA::MatDP const& a,    // 系数矩阵
                            NA::VecDP const& b,    // 右端向量
                            double    const  eps,  // 精度
                            int       const  M )   // 最多迭代次数
{
    NA::VecDP y(0.0, b.size());
    NA::VecDP x(0.0, b.size());
    double  tmp, maxXY;
    for ( int m=0; m<M ; ++m ) {
        maxXY = 0;
        for ( int i=0; i<b.size(); ++i ) {
            tmp = 0;
            for ( int j=0; j<b.size(); ++j ) {
                tmp+=x[j]*a[i][j];
            }
            tmp-=x[i]*a[i][i];
            x[i] = (b[i]-tmp)/a[i][i];   // 这里直接用新的值取代了旧的，而不是先放在向量 y 中

            tmp = fabs(x[i]-y[i]);
            if (tmp>maxXY)
                maxXY=tmp;
        }

        if (maxXY<eps)
            return y;
        y=x; // 为了判断循环终止条件，还是保留了 y
    }
    return NA::VecDP(0);
}