ヤコビ法

ラプラシアン行列を勉強するために、とりあえずヤコビ法のプログラムを作ってみました。

標準入力から最初に次元の数をもらって、その後各行を受け取るプログラムです。
反復回数は適当に１０回とかにしてみました（本当は非対角要素の大きさを見た方がいいかな）。

次回の日記ぐらいで、これを使ってラプラシアンマトリックスの固有値を解いたり、実際にスペクトル分解をやったものを日記にアップしていくつもり（多分、仕事がはっちゃけなければ…）。

//	Preprocessor
#include <iostream>
#include <sstream>
#include <cmath>
using namespace std;

//	Alias
#define ROOT_2 1.41421356
#define ITERATION_NUM 10

//	Define Template
template <typename TYPE> class Jacobi
{
public:
  TYPE**	m_eigenvector;
  TYPE**	  m_newmatrix;
  int			m_dim;

public:
  Jacobi(int dim);
  ~Jacobi();
  bool solve(TYPE** matrix) throw (const char* );
  bool copyMatrix(TYPE** buf1, TYPE** buf2);
  bool displayMatrix(TYPE** buf);
};


template<typename TYPE> 
Jacobi<TYPE>::Jacobi(int dim) :
  m_dim(dim)
{
  m_eigenvector=NULL;
  m_newmatrix=NULL;
  m_eigenvector = new TYPE*[m_dim];
  for( int i = 0; i < m_dim; i++ ) m_eigenvector[i] = new double[m_dim];
  m_newmatrix = new TYPE*[m_dim];
  for( int i = 0; i < m_dim; i++ ) m_newmatrix[i] = new double[m_dim];
  for( int i = 0; i < m_dim; i++ ){
	for( int j = 0; j < m_dim; j++ ){
	  i==j?m_eigenvector[i][j]=1:m_eigenvector[i][j]=0;
	}
  }
  displayMatrix(m_eigenvector);
}

template<typename TYPE> 
Jacobi<TYPE>::~Jacobi()
{
  if(m_eigenvector!=NULL ){
	for( int i = 0; i < m_dim; i++ ) delete [] m_eigenvector[i];
  }
}

template<typename TYPE> 
bool Jacobi<TYPE>::solve(TYPE** matrix) throw (const char* )
{
  int counter = ITERATION_NUM;
  while ( counter > 0 ){
	for( int p = 0; p < m_dim; p++ ){
	  for( int q = p+1; q < m_dim; q++ ){
		if( matrix[p][q] == 0.0 ) continue;
		copyMatrix(m_newmatrix,matrix);
		TYPE difference=(matrix[p][p]-matrix[q][q])/2.0;
		TYPE sum=(matrix[p][p]+matrix[q][q])/2.0;
		TYPE denominator=sqrt(difference*difference+matrix[p][q]*matrix[p][q]);
		TYPE c = 0;
		TYPE s = 0;
		if( difference > 0.0 ){
		  c = sqrt( 1.0+ difference/denominator)/ROOT_2;
		  s = matrix[p][q]/ (2.0*c*denominator);
		}else{
		  c = sqrt( 1.0- difference/denominator)/ROOT_2;
		  s = -1.0*matrix[p][q]/ (2.0*c*denominator);
		}

		// non diagram component
		for( int j = 0; j < m_dim; j++ ){
		  if( j == m_dim ) continue;
		  if( j == p ) continue;
		  if( j == q ) continue;
		  TYPE value1 = matrix[p][j] * c + matrix[q][j]*s;
		  m_newmatrix[p][j] = m_newmatrix[j][p] = value1;
		  
		  TYPE value2 = matrix[q][j]*c - matrix[p][j]*s;
		  m_newmatrix[j][q] = m_newmatrix[q][j] = value2;
		}

		// diagram component
		if( difference > 0 ){
		  m_newmatrix[p][p] = sum + denominator;
		  m_newmatrix[q][q] = sum - denominator;
		}else{
		  m_newmatrix[p][p] = sum - denominator;
		  m_newmatrix[q][q] = sum + denominator;
		}

		// target component
		m_newmatrix[p][q] = m_newmatrix[q][p] = 0;

		//displayMatrix(m_newmatrix);
		copyMatrix(matrix,m_newmatrix);

		// calculate vector
		for( int i = 0; i < m_dim; i++ ){
		  TYPE evip = m_eigenvector[i][p];
		  TYPE eviq = m_eigenvector[i][q];
		  m_eigenvector[i][p] = evip*c + eviq*s;
		  m_eigenvector[i][q] = eviq*c - evip*s;
		}
	  }
	}
	counter--;
  }
  return true;
}

template<typename TYPE> 
bool Jacobi<TYPE>::copyMatrix(TYPE** buf1, TYPE** buf2)
{
  for( int i = 0; i < m_dim; i ++ ){
	for( int j = 0; j < m_dim; j ++ ){
	  buf1[i][j] = buf2[i][j];
	}
  }
  return true;
}

template<typename TYPE> 
bool Jacobi<TYPE>::displayMatrix(TYPE** buf)
{
  for( int i = 0; i < m_dim; i ++ ){
	for( int j = 0; j < m_dim; j ++ ){
	  cout << buf[i][j] << " ";
	}
	cout << endl;
  }
  return true;
}

//	Define function
int main(int argc, char* argv[])
{
  char buffer[256];
  int dim;
  double** matrix;
  while( cin.getline(buffer, 256) ){
	string line(buffer);
	stringstream stout(line);
	stout >> dim;
	break;
  }
  matrix = new double*[dim];
  for( int i = 0;i < dim; i++) matrix[i] = new double[dim];
  int counter = 0;
  while( cin.getline(buffer, 256) ){
	string line(buffer);
	stringstream stout(line);
	if( counter >=  dim ){
	  cout << "Error line number = " 
		   << counter 
		   << ": dimention is " << dim << endl;
	  return 1;
	}
	for( int i = 0; i < dim; i++ ){
	  stout >> matrix[counter][i];
	}
	counter++;
  }
  Jacobi<double> jacobi(dim);

  jacobi.solve(matrix);

  cout << "eigen" << endl; 
  jacobi.displayMatrix(jacobi.m_newmatrix);
  cout << "associated vector" << endl;
  jacobi.displayMatrix(jacobi.m_eigenvector);
	
  for( int i = 0; i < dim; i++ ){
	delete [] matrix[i];
  }
  delete [] matrix;

  return 0;
}