#code(C){{
/*!
* ガウス・ザイデル法による連立方程式Ax = bの計算
* @param[in] a n×n正値対称行列
* @param[in] b n次元定数ベクトル
* @param[inout] x 引数として初期値を入れておき、ルーチン内で解を代入
* @param[inout] max_iter 最大反復回数
* @retval true 解を求めるのに成功
* @retval false 解なし
*/
static bool GaussSeidel(vector< vector<double> > &A, vector<double> &b, vector<double> &x, const int n, int max_iter = 1000)
{
vector<double> x_old;
x_old.resize(n);
double denom, dev, max_dev;
int i, j, iter = 0;
// 行列AとベクトルBの正規化
for(i = 0; i < n; i++){
denom = A[i][i];
if(denom < RX_FEQ_EPS) return false;
b[i] /= denom;
for(j = 0; j < n; j++) A[i][j] /= denom;
}
// ガウス・ザイデル反復の実行
while(true){
for(i = 0; i < n; i++){
x_old[i] = x[i];
x[i] = 0;
for(j = 0; j < n; j++){
if(j != i){
x[i] -= A[i][j]*x[j];
}
}
x[i] += b[i];
}
// 収束判定
max_dev = fabs(x_old[0]-x[0])/x[0];
for(i = 1; i < n; i++){
dev = fabs(x_old[i]-x[i])/x[i];
max_dev = (dev > max_dev) ? dev : max_dev;
}
if(max_dev <= RX_FEQ_EPS){
return true;
}
else{
iter++;
if(iter > max_iter){
return false;
}
}
}
}
}}