【数値計算】Gauss-Seidel法（ガウス・ザイデル法）のC++コード

現在「【数値計算】Gauss-Seidel法（ガウス・ザイデル法）の説明　実践編」を準備中ですが、それに先立ってGauss-Seidel法のC++コードを公開します。私はコーディングの専門家ではないのであまり綺麗ではありませんが悪しからず。皆様の勉強に役立てて頂ければ幸いです。解いている問題は

で、答えは

となります。実際に動かして遊んでみて下さい。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//inner product//
inline double ip(vector<double> &a, vector<double> &b, int &size)
{
	double value=0.;
	int i;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//Gauss-Seidel
inline void GS(vector<vector<double> > &A,vector<double> &x,vector<double> &b, int &size, double &eps)
{
	int i,k;
	double in,sum,max;
	
	for(k=0;k<1000000;k++)
	{
		max=0.0;
		for(i=0;i<size;i++)
		{
			in=x[i];
			x[i]=(b[i]-(ip(A[i],x,size)-A[i][i]*x[i]))/A[i][i];
			if(max<abs(in-x[i])) max=abs(in-x[i]);
		}
		if(eps>max) break;
	}
}

int main()
{
	int size=3;
	int i;
	double eps=pow(2.0,-50);
	double in,max;
	vector<double> b(size,0.);
	vector<double> x(size,0.);
	vector<vector<double> > A(size,(vector<double>(size,0.)));
	
	A[0][0]=3;A[0][1]=2;A[0][2]=1;
	A[1][0]=1;A[1][1]=4;A[1][2]=1;
	A[2][0]=2;A[2][1]=2;A[2][2]=5;
	b[0]=10;b[1]=12;b[2]=21;
	
	GS(A,x,b,size,eps);

	for(i=0;i<size;i++) cout<<x[i]<<endl;

        return 0;
}