ガウスの消去法

今日は
勉強会の宿題 - malibu-bulldogの日記
の宿題提出のはずが…すんません，ガウスの消去法しか終わってませんでしたって感じです。
（でも勉強会の時間は昼休み30分のみなので充分だったりする…）

で，プログラムの方を印刷して確認してました。

プログラムはこんな感じです。ちなみにピボット操作とかしてないんでダメダメです。

#include <iostream>
#include <cmath>
using namespace std;

namespace fem{
  namespace solver{
	class Matrix;
	class Vector;
	class GaussianElimination;
  }
}
using namespace  fem::solver;

class fem::solver::Matrix
{
private:
  int		_size;
  double**	   _c;

public:
  Matrix(int size);
  ~Matrix();
  
public:
  int	get_size(){return _size;}
  double** get_c(){return _c;}
};

class fem::solver::Vector
{
private:
  int  _size;
  double* _v;

public:
  Vector(int size);
  ~Vector();

public:
  int	get_size(){return _size;}
  double* get_v(){return _v;}
};

class fem::solver::GaussianElimination
{
private:
  Matrix*	_coefficient;
  Vector*	   _constant;
  Vector*           _unknown;
public:
  GaussianElimination(Matrix* coefficient, Vector* constant, Vector* unknown);
  ~GaussianElimination();
  int solve(void);
  void debugPrint();
};

fem::solver::Matrix::Matrix(int size) :
  _size(size)
{
  _c = new double*[_size];
  for( int i = 0; i < _size; i++ ){
	_c[i] = new double[_size];
	memset(_c[i], 0, sizeof(double)*_size);
  }
}

fem::solver::Matrix::~Matrix()
{
  for( int i = 0; i < _size; i++ ){
	delete [] _c[i];
  }
  delete [] _c;
}

fem::solver::Vector::Vector(int size) :
  _size(size)
{
  _v = new double[_size];
  memset(_v, 0, sizeof(double)*_size);
}

fem::solver::Vector::~Vector()
{
  delete [] _v;
}

fem::solver::GaussianElimination::GaussianElimination(Matrix* coefficient,Vector* constant,Vector* unknown) :
  _coefficient(coefficient),
  _constant(constant),
  _unknown(unknown)
{
}


fem::solver::GaussianElimination::~GaussianElimination()
{
}


void fem::solver::GaussianElimination::debugPrint(void)
{

  for( int i = 0; i < _coefficient->get_size(); i++ ){
	for( int j = 0; j < _coefficient->get_size(); j++ ){
	  cout.width(5);
	  cout << _coefficient->get_c()[i][j] << " ";
	}
	cout.width(5);
	cout << _unknown->get_v()[i] << " = ";
	cout.width(5);
	cout << _constant->get_v()[i] << endl;
  }	
}


int fem::solver::GaussianElimination::solve(void)
{
  cout << "input" << endl;
  debugPrint();
  int n = _coefficient->get_size();
  double** c = _coefficient->get_c();
  double* b = _constant->get_v();

  // m is number of the line eliminate terms.
  for( int m = 1; m < n; m++ ){

	// i is number of the eliminated term .
	for( int i = 0; i < m; i++ ){

	  // j is after the liminated term.
	  for( int j = i+1; j < n; j++ ){
		c[m][j] = c[m][j] - c[m][i]*(c[i][j]/c[i][i]);
	  }
	  b[m] = b[m] - c[m][i]*(b[i]/c[i][i]);

	  // Targeted liminated term are zero.
	  c[m][i] = 0;
	}
  }
  cout << "After forward elimination" << endl;
  debugPrint();

  for( int m = n-1; m >= 0; m-- ){
	double tmp = b[m];
	for( int i = m+1; i < n; i++ ){
	  tmp -= c[m][i]*_unknown->get_v()[i];
	}
	_unknown->get_v()[m] = tmp/c[m][m];
  }

  cout << "After backward subsitution" << endl;
  debugPrint();
}


int main(int argc, char* argv[])
{  
  Matrix* mat = new Matrix(3);
  Vector* con = new Vector(3);
  Vector* unk = new Vector(3);

  mat->get_c()[0][0] = 1;
  mat->get_c()[0][1] = -1;
  mat->get_c()[0][2] = 1;
  mat->get_c()[1][0] = 1;
  mat->get_c()[1][1] = 1;
  mat->get_c()[1][2] = -1;
  mat->get_c()[2][0] = -1;
  mat->get_c()[2][1] = 1;
  mat->get_c()[2][2] = 1;

  con->get_v()[0] = 0;
  con->get_v()[1] = -2;
  con->get_v()[2] = 2;

  GaussianElimination* method = new GaussianElimination(mat, con, unk);
  method->solve();
  return 0;
}