本系列博客將利用C++實現一系列數值演算法。數值演算法離不開矩陣，但是C++並未自帶矩陣這一對象，直接使用數組又會帶來諸多不便，因此我們需要做一些預備工作————編寫一個矩陣類，實現矩陣的基本功能。一般來說，讀者可以直接使用Eigen庫進行矩陣計算，從頭開始造輪子僅僅是為了滿足筆者個人的需要。 #一、成 ...
本系列博客將利用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;
}
到此已經實現了矩陣的一些基本功能，在接下來的博客中，我將完善其他矩陣的功能，並且逐步實現一些古老的數值代數演算法。
