C++数值计算——矩阵类的实现(一)

本系列博客将利用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.Number_of_column != this->Number_of_column || A.Number_of_row != this->Number_of_row){
                for(int i = 0; i < this->Number_of_row; i++){
                    delete[] p_Matrix[i];
                }
                delete[] p_Matrix;
                p_Matrix = new double*[A.Number_of_row];
                for( int i = 0; i < A.Number_of_row; i++){
                    p_Matrix[i] = new double[A.Number_of_column]; 
                }


                this->Number_of_row = A.Number_of_row;
                this->Number_of_column = A.Number_of_column; 
                for(int i = 0; i < this->Number_of_row; i++){
                    for(int j = 0; j < this->Number_of_column; j++){
                        this->p_Matrix[i][j] = A.p_Matrix[i][j];
                    }
                }
            }
            else{
                for(int i = 0; i < this->Number_of_row; i++){
                    for(int j = 0; j < this->Number_of_column; j++){
                        this->p_Matrix[i][j] = A.p_Matrix[i][j];
                    }
                }                              
            }      
            return *this;      
        }
        //重载减法
        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*(double value){
            Matrix tempMatrix(this->Number_of_row, this->Number_of_column);

            for(int i = 0; i < this->Number_of_row; i++){
                for(int j = 0; j < this->Number_of_column; j++){
                    tempMatrix.p_Matrix[i][j] = value*this->p_Matrix[i][j];
                }
            }
            return tempMatrix;
        }
        friend Matrix operator*(double value, const Matrix &A){
            Matrix tempMatrix(A.Number_of_row, A.Number_of_column);

            for(int i = 0; i < A.Number_of_row; i++){
                for(int j = 0; j < A.Number_of_column; j++){
                    tempMatrix.p_Matrix[i][j] = value*A.p_Matrix[i][j];
                }
            }
            return tempMatrix;            
        }
        //矩阵相乘
        friend Matrix operator*(Matrix &A, Matrix &B){
            Matrix tempMatrix(A.Number_of_row, B.Number_of_column);    
            if(A.Number_of_column != B.Number_of_row){
                cout<<"Invalid!"<<endl;
            }
            else{
                double s;
        
                for(int i = 0; i < A.Number_of_row; i++){
                    for(int k = 0; k < A.Number_of_column; k++){
                        s=A.p_Matrix[i][k];
                        for(int j = 0; j < B.Number_of_column; j++){
                            tempMatrix.p_Matrix[i][j] += s*B.p_Matrix[k][j];
                        }
                    }
                }
            }
            return tempMatrix;
        }

到此已经实现了矩阵的一些基本功能,在接下来的博客中,我将完善其他矩阵的功能,并且逐步实现一些古老的数值代数算法。

热门相关:总裁别再玩了   闺范   万古至尊   神算大小姐   闺范