本系列博客将利用C++实现一系列数值算法。数值算法离不开矩阵，但是C++并未自带矩阵这一对象，直接使用数组又会带来诸多不便，因此我们需要做一些预备工作————编写一个矩阵类，实现矩阵的基本功能。一般来说，读者可以直接使用Eigen库进行矩阵计算，从头开始造轮子仅仅是为了满足笔者个人的需要。

一、成员组成

回顾矩阵的定义，我们仅需三个量就可以具体描述一个矩阵：行指标，列指标，对应位置的元素。因此我们在类Matrix（下文就如此称呼了，和代码保持一致）中定义三个数据成员：行指标，列指标，一个二重指针。

        typedef unsigned int Index;
        class Matrix{
            private:
                Index Number_of_row;//行数
                Index Number_of_column;//列数
                double **p_Matrix;//二重指针构造矩阵
}

二、基本功能的分析与实现

按一般类的定义，类Matrix需要有构造函数、析构函数和拷贝函数。构造函数生成矩阵时，矩阵的每一个位置都需要被赋值，最合适的默认值莫过于0，因此在用户未指定的情况下，默认每个值为零；如果用户指定了某个值a，则将每个位置赋值a。因此，如下创建构造函数：

        Matrix( Index num_of_row, Index num_of_column){ //一般矩阵，默认为全零矩阵。
            Number_of_row = num_of_row;
            Number_of_column = num_of_column;
            p_Matrix = new double*[num_of_row];
            for( int i = 0; i < num_of_row; i++){
                p_Matrix[i] = new double[num_of_column];
            }
            for( int i = 0; i < num_of_row; i++){
                for( int j = 0; j < num_of_column; j++){
                    p_Matrix[i][j] = 0;
                }
            }
        }
        Matrix( Index num_of_row, Index num_of_column, double value){ //一般矩阵，默认为全为value
            Number_of_row = num_of_row;
            Number_of_column = num_of_column;
            p_Matrix = new double*[num_of_row];
            for( int i = 0; i < num_of_row; i++){
                p_Matrix[i] = new double[num_of_column];
            }
            for( int i = 0; i < num_of_row; i++){
                for( int j = 0; j < num_of_column; j++){
                    p_Matrix[i][j] = value;
                }
            }
        }        

对应的析构函数和拷贝函数如下：

        //析构函数
        ~Matrix(){
            for( int i = 0; i < Number_of_row; i++){
                delete[] p_Matrix[i];
            }
            delete[] p_Matrix;
        }
        //拷贝函数
        Matrix( const Matrix &Copy_Matrix){
            Number_of_row = Copy_Matrix.Number_of_row;
            Number_of_column = Copy_Matrix.Number_of_column;
            for(int i = 0; i < Number_of_row; i++){
                p_Matrix[i] = new double[Number_of_column];
            }
            for( int i = 0; i < Number_of_row; i++){
                for( int j = 0; j < Number_of_column; j++){
                    p_Matrix[i][j] = Copy_Matrix.p_Matrix[i][j];
                }
            }
        }

对于类Matrix而言，它自然必须有能显示和改变元素值的功能，我们将这个需求交给以下两个函数：

        //输出矩阵
        void Print_Matrix(){
            for( int i = 0; i < Number_of_row; i++){
                for( int j = 0; j < Number_of_column; j++){
                    cout << p_Matrix[i][j] << " ";
                }
                cout << endl;
            }            
        }
        //改变元素，或者赋值
        void Change_Value(Index m, Index n, double x){
            if( m >= Number_of_row || n >= Number_of_column){
                cout << "Invalid! The index exceeds the range!" << endl;
            }else{
                p_Matrix[m][n] = x;
            }   
        }

当然，这种赋值的方式极度不友好，我们将在未来改写这一赋值方式。此外，我们常常需要取出某一位置的某一值，我们将这一功能交给下面这个函数：

        //取出某一元素
        friend double Get_element(Matrix &A, int m, int n){
            if(m >= 0 && m < A.Number_of_row && n >= 0 && n < A.Number_of_column){
                return A.p_Matrix[m][n];                
            }
            else{
                cout<<"The index exceeds the range!"<<endl;
                return 0;
            }
        }

三、运算符的重载

考虑矩阵的基本运算：+、-、*（包含数乘与矩阵相乘）与=（赋值），我们需要重载运算符。值得注意的是，在矩阵乘法中我们采用的循环策略，这是考虑到C++从左往右扫描数组的特点，具体可以参考《Matrix Computation》。

        //重载运算符
        //重载加法
        Matrix operator+ (const Matrix &A){
            Matrix tempMatrix(A.Number_of_row, A.Number_of_column);

            if (A.Number_of_column != this->Number_of_column || A.Number_of_row != this->Number_of_row){
                cout << "Matrices are in different size!" << endl;
            }
            else{
                for(int i = 0; i < this->Number_of_row; i++){
                    for(int j = 0; j < A.Number_of_column; j++){
                        tempMatrix.p_Matrix[i][j] = this->p_Matrix[i][j] + A.p_Matrix[i][j];
                    }
                }  
            }
            return tempMatrix;
        }
        //重载赋值
        Matrix &operator=(const Matrix &A){
            if (A.Nu