ルギア君の戯言

雑多な記事。

行列

まえに言っていた行列計算用のライブラリのうち、かけ算だけできた。


というか、それが先週の木曜日に課題としてあったので、超ラッキーと思って、取り組んでいたのだが、バグだらけで結局、今日終わった(ぁ


ちょっと長いけど解析したければ、頑張って読んでみてくれ。


ちなみに、行列はまず派生元の vector を行ベクトルを繋げたものとして、x_size によって分割するという考えのもと、作成されたので、ちょっとややっこしいと言えば、ややっこしい。
もし、もっと本格的に作るなら、vector のもとになっている allocator を調べて、行列用にSTLコンテナをつくる必要があるかもね。
また、たぶん、このままだと、もしかしたら、x_size の管理が杜撰(ずさん)だっていうかも知れないけどその場合は、どうすればちゃんと管理が行き届くかアドバイスもいっしょにくれると有り難いです。
コピーコンストラクタはありませんが、コピーの場合は x_size はちゃんとコピーしてくれるようですが・・・(まあ、当然


では、まずは、ヘッダコード。

#ifndef _MATRIX_H_
#define _MATRIX_H_

#include <vector>

using namespace std;

template <typename T, typename A = allocator<T> >
class matrix : public vector<T, A> {
	typedef matrix<T, A> matrix_type;
	typedef vector<T, A> mat_vec_type;

	int x_size;
    public:

	matrix() : mat_vec_type() { x_size = 0; }
	matrix(int a) : mat_vec_type(a) { x_size = a; }
	matrix(int x, int y) : mat_vec_type(x * y) { x_size = x; }

	template<typename Iterator>
	matrix(Iterator a, Iterator b, int x) : mat_vec_type(a, b) { x_size = x; }

	int getXsize() {
		return x_size;
	}

	int getYsize() {
		return this->size() / x_size;
	}

	void setXsize(int x) {
		x_size = x;
	}
	
	T& operator () (int x, int y) {
		return *(this->begin() + x + getXsize() * y);
	}

	matrix_type operator *(matrix_type& b) {
		matrix c(b.getXsize(), this->getYsize());
		if(this->getXsize() != b.getYsize()) {
			return matrix_type();
		}
		
		for(int i = 0; i < this->getYsize(); i++) {
			for(int j = 0; j < b.getXsize(); j++) {
				c(j, i) = 0.0;
				for(int k = 0; k < b.getYsize(); k++) {
// 					cout << "c(" << j << ", " << i << ") += ";
// 					cout << "a(" << k << ", " << i << ")(" << (*this)(k, i) << ")";
// 					cout << " * ";
// 					cout << "b(" << j << ", " << k << ")(" << b(j, k) << ")";
// 					cout << " = " << (*this)(k, i) * b(j, k) << endl;
					c(j, i) = c(j, i) + (*this)(k, i) * b(j, k);
				}
			}
		}
		return c;
	}
};

#endif // _MATRIX_H_

途中のコメントアウトされている cout の羅列はデバッグ用に使っていたものです。参考のために残しておきました。
(サイズ不整合で)かけ算に失敗した場合は、空行列が帰ってくるので empty 関数(vector に付属)で調べてあげてください。

使用方法

  double mat_a[9] = {
    1, 0, 0,
    0, 1, 0,
    0, 0, 1,
  };
  matrix<double> A(mat_a, mat_a + 9, 3);   // 最後は x 方向のサイズ。配列からは取得できないので。
  matrix<double> B(3, 3);                  // 3 X 3 の行列
  matrix<double> C(3);                     // 大きさが 3 の行ベクトル
  
  B = A * A;                               // かけ算はこのように普通にできます。
  if(B.empty()) {
     cout << "内部エラー";
     return;
  }

  cout << B(1, 1) << endl;                 // (1, 1) 成分を表示。配列同様、最初の値は 0。
  cout << B[1 + B.getXsize() * 1] << endl; // (1, 1) 成分を表示。vector に付いてくる operator [] を使ったもの。
  cout << B.getYsize() << endl;            // Y 方向のサイズを表示。割り切れない状態でも例外スローはしません(直せ

サンプルプログラム

#include <iostream>
#include "matrix.h"

using namespace std;

double A[4][4] = {
	{0,1,0,0},
	{0,0,1,0},
	{0,0,0,1},
	{0,0,0,0},
};

double B[8] = {
	1,2,3, 5,
	3,1,2,-1,
};

double C[8] = {
	 0, 3,
	 1, 2,
	-3, 4,
	-2,-1,
};


void disp_matrix(matrix<double>& m);

int main() {
	matrix<double> a(&(A[0][0]), &(A[3][3])+1, 4);
	matrix<double> b;
	disp_matrix(a);
	b = a * a;	// a^2
	disp_matrix(b);
	b = b * a;	// a^3;
	disp_matrix(b);
	b = b * a;	// a^4;
	disp_matrix(b);

	matrix<double> c(&(B[0]), &(B[7])+1, 4);
	matrix<double> d(&(C[0]), &(C[7])+1, 2);
	disp_matrix(c);
	disp_matrix(d);
	c = c * d;	// c*d
	disp_matrix(c);
	b = c * c;	// (c*d)^2;
	disp_matrix(b);
	b = b * c;	// (c*d)^3;
	disp_matrix(b);
	b = b * c;	// (c*d)^4;
	disp_matrix(b);

	double D[9] = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0};
	matrix<double> e(&(D[0]), &(D[8])+1, 3);
	disp_matrix(e);
	e = e * e;
	disp_matrix(e);
	e = e * e;
	disp_matrix(e);
}

void disp_matrix(matrix<double>& m) {
	for(int i = 0; i < m.getYsize(); i++) {
		cout << "( ";
		for(int j = 0; j < m.getXsize(); j++) {
			cout << m(j, i);
			cout << ", ";
		}
		cout << "\b\b )" << endl;
	}
	cout << endl;
}


あとは、とりあえず、転置、エルミート行列、逆行列、正則判定、行列式基本変形、階数、連立一次方程式とみなして解く、上三角化、対角化、(ジョルダン)標準形変換、・・・



やはり、行列は多いなw