C++實現的線性代數矩陣計算

/**
* 線性代數矩陣計算
* 實現功能：行列向量獲取，子矩陣獲取，轉置矩陣獲取，
* 行列式計算，轉置伴隨矩陣獲取，逆矩陣計算
*
* Copyright 2011 Shi Y.M. All Rights Reserved
*/

#ifndef MATICAL_TMATRIX_H_INCLUDED
#define MATICAL_TMATRIX_H_INCLUDED

#include <memory.h>

namespace matical
{

template<typename T>
class TMatrix
{
private:
T * m_Elems;
int m_Rows;
int m_Cols;

public:
typedef T ElemType;

public:

TMatrix() : m_Elems(NULL), m_Rows(0), m_Cols(0) {}
TMatrix(int rows, int cols) : m_Rows(rows), m_Cols(cols) {
m_Elems = new ElemType[m_Rows * m_Cols];
memset(m_Elems, 0, sizeof(ElemType) * m_Rows * m_Cols);
}
// 生成n階方陣
TMatrix(int n) : m_Rows(n), m_Cols(n){
m_Elems = new ElemType[m_Rows * m_Cols];
memset(m_Elems, 0, sizeof(ElemType) * m_Rows * m_Cols);
}

TMatrix(const TMatrix& m) : m_Rows(m.Rows()), m_Cols(m.Cols()) {
m_Elems = new ElemType[m_Rows * m_Cols];
memcpy(m_Elems, m.m_Elems, sizeof(ElemType) * m_Rows * m_Cols);
}
virtual ~TMatrix() {
if ( m_Elems ) delete[] m_Elems;
m_Elems = NULL;
m_Rows = 0;
m_Cols = 0;
}

public:

int Rows() const { return m_Rows;}
int Cols() const { return m_Cols;}

// 是否為方陣
bool IsSquare() const { return ( m_Rows > 0 && (m_Rows == m_Cols)); }

// 取行向量
TMatrix<T> Row(int row) const {
TMatrix m(1, m_Cols);
for(int i = 0; i < m_Cols; i++) {
m(0, i) = (*this)(row, i);
}
return m;
}

// 取列向量
TMatrix<T> Col(int col) const {
TMatrix m(m_Rows, 1);
for(int i = 0; i < m_Rows; i++) {
m(i, 0) = (*this)(i, col);
}
return m;
}

public:

// 將當前矩陣設置為單位矩陣
void Identity() {
for(int r = 0; r < this->m_Rows; r++) {
for(int c = 0; c < this->m_Cols; c++) {
if ( r == c ) (*this)(r, c) = 1;
else (*this)(r, c) = 0;
}
}
return ;
}

// 當前矩陣設置為零矩陣
void Zero() {
int sz = m_Rows * m_Cols;
for (int i = 0; i < sz; i++ ) m_Elems[i] = 0;
return ;
}

// 獲取指定范圍的子矩陣，rrow>=lrow, rcol>=lcol
TMatrix<T> SubMatrix(int lrow, int lcol, int rrow, int rcol) const {
TMatrix m(rrow - lrow + 1, rcol - lcol + 1);
for ( int r = 0; r < m.m_Rows; r++) {
for (int c = 0; c < m.m_Cols; c++) {
m(r, c ) = (*this)(r + lrow, c + lcol);
}
}
return m;
}

// 獲取當前矩陣元素(row, col)的子矩陣（去掉元素（row,col）所在的行和列）
TMatrix<T> SubMatrixEx(int row, int col) const {
int rr = 0;
int cc = 0;
TMatrix m(this->m_Rows - 1, this->m_Cols - 1);
for ( int r = 0; r < m.m_Rows; r++, rr++) {
if ( r == row) rr++;
cc = 0;
for (int c = 0; c < m.m_Cols; c++, cc++) {
if ( c == col ) cc++;
m(r, c) = (*this)(rr, cc);
}
}
return m;
}

// 計算當前矩陣的轉置矩陣
TMatrix<T> Transpose() const {
printf("Transpose\n");
TMatrix m(this->m_Cols, this->m_Rows);
for (int r = 0; r < m.m_Rows; r++ ) {
for ( int c = 0; c < m.m_Cols; c++ ) {
m(r, c) = (*this)(c, r);
}
}
return m;
}

// 計算當前矩陣的行列式(遞歸方式), 計算行列式的矩陣，須是n階方陣
ElemType Det() const {
ElemType det = 0;
if ( m_Rows == 1) {
det = (*this)(0, 0);
} else if ( m_Rows == 2) {
det = (*this)(0, 0) * (*this)(1, 1) - (*this)(0 ,1) * (*this)(1, 0);
} else {
for (int r = 0; r < m_Rows; r++) {
TMatrix<ElemType> m = this->SubMatrixEx(r, 0);
if ( r % 2 ) det += (-1) * (*this)(r, 0) * m.Det();
else det += (*this)(r, 0) * m.Det();
}
}
return det;
}

// 計算當前矩陣的轉置伴隨矩陣，當前矩陣須為n階方陣
TMatrix<T> Adj() const {
TMatrix<T> m(m_Rows, m_Cols);
for (int r = 0; r < m_Rows; r++) {
for (int c = 0; c < m_Cols; c++) {
if ( (r + c) % 2 ) m(c, r) = -1 * this->SubMatrixEx(r, c).Det();
else m(c, r) = this->SubMatrixEx(r, c).Det();
}
}
return m;
}

// 計算當前矩陣的逆矩陣，當前矩陣須為n階方陣
TMatrix<T> Inv() const {
ElemType det = this->Det();
TMatrix<T> adj = this->Adj();
return adj * (1 / det);
}

public:

ElemType operator()(int row, int col) const { return m_Elems[row * m_Cols + col];};
ElemType& operator()(int row, int col) { return m_Elems[row * m_Cols + col]; }

public:

// 兩矩陣相加，必須具有相同的行數以及列數
TMatrix<T> operator+(const TMatrix& m) const {
TMatrix mr(this->m_Rows, this->m_Cols);
for(int r = 0; r < this->m_Rows; r++) {
for(int c = 0; c < this->m_Cols; c++) {
mr(r, c) = (*this)(r, c) + m(r, c);
}
}
return mr;
} // operator+(const TMatrix& m) const

// 兩矩陣相減，必須具有相同的行數以及列數
TMatrix<T> operator-(const TMatrix& m) const {
TMatrix mr(this->m_Rows, this->m_Cols);
for(int r = 0; r < this->m_Rows; r++) {
for(int c = 0; c < this->m_Cols; c++) {
mr(r, c) = (*this)(r, c) - m(r, c);
}
}
return mr;
} // operator-(const TMatrix& m) const

// 矩陣與常數相乘(數乘)
TMatrix<T> operator*(ElemType v) const {
TMatrix mr(this->m_Rows, this->m_Cols);
for(int r = 0; r < this->m_Rows; r++) {
for(int c = 0; c < this->m_Cols; c++) {
mr(r, c) = (*this)(r, c) * v;
}
}
return mr;
} // operator*(ElemType v) const

// 矩陣相乘(當前矩陣列數須等於參數矩陣的列數)，
// 結果矩陣行數等於當前矩陣，列數等於參數矩陣的列數
TMatrix<T> operator*(const TMatrix& m) const {
TMatrix mr(this->m_Rows, m.m_Cols);
for (int r = 0; r < this->m_Rows; r++) {
for (int c = 0; c < m.m_Cols; c++) {
for (int i = 0; i < this->m_Cols; i++) {
mr(r, c) += ( (*this)(r, i) * m(i, c) );
printf("(%d, %d)\n", r, c);
}
}
}
return mr;
} // operator*(const TMatrix& m) const
}; // class TMatrix

} // namespace matical

#endif // MATICAL_TMATRIX_H_INCLUDED