ブルーシフト

大学院生の幸せな日々

Eigenを触ってみる

*勉強メモです

 

昨日導入した行列演算ライブラリ「Eigen」で遊んだのでメモ。

と言っても、やったことはエルミート行列(Hermitian matrix)の対角化だけ。最初は3行3列の実数行列で試したけど複素数行列はどうすんのかなーと思って作った。

C++を触るのは初めてだったのでいい練習にもなったかなと言う印象。

 

問題に使った行列は某マセ○の参考書に載っていた↓を用いた。

f:id:pa_yoyoyo:20160708141507p:plain固有値固有ベクトルを求めよ。

 

結論だけ述べるとこんくらいシンプルなコードで完成↓

  #include <iostream>

#include <Eigen/Dense> //今回使うライブラリ

   
  using namespace std;
  using namespace Eigen;
  int main(){
  complex<double> a;
  a=complex<double>(0,1);
  Matrix<complex<double>, 3, 3>A; //複素数行列Aの型指定
  A << 1.0,sqrt(2.0)*a,0,-sqrt(2.0)*a,2.0,-2.0*a,0,2.0*a,1.0; //行列要素の代入
   
  ComplexEigenSolver<Matrix<complex<double>,3,3>> es(A);
  if (es.info() != Success) abort();
  cout << "固有値\n"
  << es.eigenvalues() << endl;
  cout << "固有ベクトル\n"
  << es.eigenvectors() << endl;
  getchar();
  return 0;
  }

最初#include<cmath>やら#include<complex>やらを置いていたけど消しても実行できたから#include<Eigen/Dence>に入っているんだと思う(知らんが)。

System("pause")じゃなくてgetchar()と言う便利なものも覚えてほくほくしている。

 

コーディングについて

実数の行列だと次元の指定

Matrix<complex<double>, 3, 3>A;

の部分を

Matrix3d A;

と表記しても問題なかった。しかし複素数を扱いたい場合はどうも

Matrix<complex<double>, 3, 3>A;

と直接指定してやる必要があるらしい。こうしないとエラる。

同じようなことだが実数だと

ComplexEigenSolver<Matrix<complex<double>,3,3>> es(A);

EigenSolver<Matrix<complex<double>,3,3>> es(A);

としていたが、複素数ではここも変えなきゃいけないっぽい。

 

計算結果↓ (a,b)のaが実部,bが虚部。

f:id:pa_yoyoyo:20160708143135p:plain

 

思ったんだが、平方根で結果を表示できないんだろうか。1行1列の-√10/√15を-0.816497とか書かれても微妙に分かり辛くないか、と思うんだが・・・。

 

 

また、このブログにコードを張っ付けたのには↓のGistと言うサイトを使った。

gist.github.com

 「Filename including extension...」に「ファイル名(適当).cpp」と打ち込んでCreateすると色付きで出力される、パワポやブログに張るときに便利だと思う。

 

C++にかなり手を焼くと思っていたが、Cの記述を進化させただけで後は同じようなものなのかな?案外すぐ完成した。Eigenを使ってブイブイState-of-the-artしまくろうと思う。

 

 

ではまた。