C++实现的线性代数矩阵计算

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