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;
}
到此已经实现了矩阵的一些基本功能，在接下来的博客中，我将完善其他矩阵的功能，并且逐步实现一些古老的数值代数算法。