[C++] 행렬 라이브러리 구현 (simple matrix library)

matrix library로 코드는 맨 아래 깃헙 링크에 공개 하였습니다.

 

행렬 연산을 위한 cpp 라이브러리 입니다.

 

먼저 constructor는 다음과 같습니다. 

KwanMat::KwanMat(const int rows, const int cols, float matrices[]) : numRows(rows), numCols(cols) 
{
	matNums = new float[numRows * numCols];
	if (matrices == 0)
	{
		for (unsigned int i = 0; i < numRows * numCols; i++)
			matNums[i] = 0;
	}
	else
	{
		for (unsigned int i = 0; i < numRows * numCols; i++)
			matNums[i] = matrices[i];
	}
};

 

2d vector를 사용하지 않고 1d array를 이용해 구현하였습니다.

 

그 이유는, 약간의 메모리(vector class로 구현하는 경우)와 약간의 속도(2d로 구현하는 경우 accesss time) 때문인데, 큰차이가 나지는 않습니다.

 

다음으로 구현한 내용은 uniform calculation과 memberwise calculation으로, uniform 연산의 경우 member function으로, memeberwise calc는 operation overloading으로 구현하였습니다. 자세한 내용은 생략하고 선언은 다음과 같습니다.

 

    //uniform calculation of float
    void add(const float value);
    void sub(const float value);
    void mul(const float value);
    void div(const float value);

    //memberwise calculation
    KwanMat& operator+(const KwanMat& ref);
    KwanMat& operator-(const KwanMat& ref);
    KwanMat& operator*(const KwanMat& ref);

 

다음으로 Trasnpose는 임시 array를 정적으로 할당하고 1024 이상의 사이즈는 고려하지 않았습니다.

좀 더 general 한 implemenation을 위해서는(1024 이상이 올 수 도 있는 경우) 동적할당 합니다.

KwanMat& KwanMat::T()
{
	//for general purpose
	//float* tempNums = new float[numRows*numCols];

	float tempNums[1024];

	//copy temporally
	for (unsigned int i = 0; i < numRows * numCols; i++)
		tempNums[i] = matNums[i];
	for (unsigned int i = 0; i < numRows; i++)
	{
		for (unsigned int j = 0; j < numCols; j++)
		{
			matNums[i*numCols+j] = tempNums[j * numCols + i];
		}
	}
	//delete[] tempNums;

	int tempN = numRows;
	numRows = numCols;
	numCols = tempN;

	return *this;
}

 

다음으로 행렬곱입니다.

행렬곱은 행렬1의 column 사이즈와 행렬2의 row 사이즈가 같아야 하므로 이를 확인하고, 계산하였습니다.

마찬가지로 1d array이므로 계산식은 

theMat.matNums[row*ref.numCols+col] += matNums[row*numCols+inner] * ref.matNums[inner*ref.numCols+col];

입니다.

KwanMat& KwanMat::MatMul(KwanMat& ref)
{
	if (numCols != ref.numRows)
	{
		throw std::out_of_range("Matrix multiplication should match the size");
	}

	KwanMat& theMat = Zero(numRows, ref.numCols);

	for (unsigned int row = 0; row < numRows; row++) {
		for (unsigned int col = 0; col < ref.numCols; col++) {
			for (unsigned int inner = 0; inner < numCols; inner++) {
				theMat.matNums[row*ref.numCols+col] += matNums[row*numCols+inner] * ref.matNums[inner*ref.numCols+col];
			}
		}
	}
	return theMat;
}

 

다음으로 determinant 연산입니다.

이때 부터 1d array로 implement한것을 굉장히 후회했습니다...

 

먼저 wrapper function에서 간단하게 square인지 확인하고 1x1 행렬인지 확인합니다.

다음으로 Laplace expansion 을 통해, 6x6보다 작은 행렬에 대한 det를 구합니다. 

recursive method이므로 행렬의 크기만 변경하면 얼마든지 더 큰 행렬에 대해서도 determinant를 구할 수 있습니다.

float KwanMat::RecursiveDet(const float mat[25], const int n)
{
	float det = 0;
	float subMat[25];
	if (n == 2)
		return mat[0] * mat[3] - mat[1] * mat[2];
	else {
		for (int x = 0; x < n; x++) {
			int subi = 0;
			for (int i = 1; i < n; i++) {
				int subj = 0;
				for (int j = 0; j < n; j++) {
					if (j == x)
						continue;
					subMat[subi*(n-1)+subj] = mat[i*n+j];
					subj++;
				}
				subi++;
			}
			det += x % 2 ? -mat[x] * RecursiveDet(subMat, n - 1) : mat[x] * RecursiveDet(subMat, n - 1);
			//det = det + (pow(-1, x) * mat[x] * RecursiveDet(submatrix, n - 1));
		}
	}
	return det;
}

먼저 2x2인 경우 a*d-b*c로 구하고, 그 이상인경우는 submatrix로 나누어 계산하게 됩니다.

 

예를들어 3x3인경우,

첫번째 row에서 x를 하나씩 뽑습니다.

하나씩 뽑은 x에 대해서 x와 row,column이 겹치지 않은 submatrix를 만들고 해당 submatrix의 determinant에 x를 곱하는 방식을 통해 임시 det값을 구합니다. 이때,  (-1)^(row+col)을 통해 부호를 결정하는데 여기서는

//fast
x % 2 ? -mat[x] * RecursiveDet(subMat, n - 1) : mat[x] * RecursiveDet(subMat, n - 1);

//slow
det = det + (pow(-1, x) * mat[x] * RecursiveDet(submatrix, n - 1));

위와 같은 방식으로 해당 부호결정을 진행하였습니다.

 

이제 역행렬을 함수를 구현하도록 하겠습니다.

수반행렬, 여인자를 이용해 구하는 방식으로 진행하였습니다.

KwanMat& KwanMat::Inverse()
{
	float det = KwanMat::DET();
	if (det == 0) {
		throw std::runtime_error("Determinant is 0");
	}
	float dInv = 1 / det;
	
	float adj[25]; 
	ADJ(matNums, adj);
	for (int i = 0; i < numRows* numRows; i++)
		matNums[i] = adj[i] * dInv;

	return *this;
}

먼저, det를 구함으로 역행렬이 존재하는지 확인합니다.

다음으로 adjoint 행렬을 구하고 여기에 1/det를 통해 inverse행렬을 구합니다.

 

adj행렬은 다음과 같이 구하게 됩니다.

void KwanMat::ADJ(float M[25], float adj[25])
{
	int s = 1;
	float t[25] = {0};
	for (int i = 0; i < numRows; i++) {
		for (int j = 0; j < numRows; j++) {
			//To get cofactor of M[i][j]
			Cofactor(M, t, i, j, numRows);
			s = ((i + j) % 2 == 0) ? 1 : -1; //sign of adj[j][i] positive if sum of row and column indexes is even.
			adj[j*numRows+i] = (s) * (RecursiveDet(t, numRows-1)); //Interchange rows and columns to get the transpose of the cofactor matrix
		}
	}
}

 

여기서 여인자는, 다음과 같습니다.

 

void KwanMat::Cofactor(const float M[25],float t[16],int p,int q, int n)
{
	int i = 0, j = 0;
	for (int r = 0; r < n; r++) {
		for (int c = 0; c < n; c++) //Copy only those elements which are not in given row r and column c
			if (r != p && c != q) {
				t[i*(n-1)+j++] = M[r*n+c]; //If row is filled increase r index and reset c index
				if (j == n - 1) {
					j = 0; i++;
				}
			}
	}
}

 

이외에도, symmetric 행렬인지, square matrix인지 파악하는 함수와

identity matrix를 생성하는 함수, zero matrix를 생성하는 함수 등을 만들었습니다.

 

 

https://github.com/kwanyun/KwanMatrix-libarary

GitHub - kwanyun/KwanMatrix-libarary: simple Matrix Library