再来一个矩阵类myVector,这是最普遍意义上的矩阵,并不是特殊的对称矩阵、对角矩阵、稀疏矩阵……。成员变量只有四个,即行维数、列维数、总体元素数目、还有double指针用于动态申请释放内存来存储矩阵。构造函数可用行数列数和元素值构造(矩阵中元素都是同一个值),也可用另一个矩阵来做复制构造(Copy constructor),析构函数中释放内存空间。重载赋值符=,用于两矩阵间的赋值操作,重载操作符()用于引用矩阵的(i,j)位置的元素,重载输入流>>用于输入矩阵,重载输出流<<用于输出矩阵,重载运算符+,+=用于矩阵和标量或矩阵的加,自加操作,重载运算符*,*=用于矩阵和标量或矩阵相乘,矩阵和标量的自乘运算,友元函数determinant是求矩阵(方阵)的行列式,用的是最笨的方法,即行列式等于任意一行(这里选的是第1行)中的每一个元素与其代数余子式乘积的加和,显然,该方法要用到递归处理,即,函数的自我调用。
注意:矩阵的转置、求逆、求特征值、求解Ax=B等操作并没有实现。
矩阵类myMatrix的框架如下:
成员变量
- 行维数
- 列维数
- 总体元素数目
- double指针用于动态申请释放内存空间,来存储矩阵元素
成员函数
- 构造函数1,行维数、列维数、元素值构造,缺省值为0
- 构造函数2,由另一个矩阵构造
- 析构函数,释放内存空间
- 赋值运算符=,左矩阵 = 右矩阵,将右侧矩阵的值赋给给左侧矩阵
- 判断给定的行和列标识是否合规,即,是否在矩阵行列维数范围内
- 重载操作符(),返回矩阵i行j列位置元素的引用,即,可直接用A(i,j)来获得相应位置的元素
- 提取矩阵行维数、列维数、总体元素数目的getNumRow,getNumCol,getNumCmp函数
- 提取矩阵i行j列元素值的getAij函数
- 设置矩阵i行j列元素值的setAij函数
- 输出流<<,输出矩阵
- 输入流>>,输入矩阵
- 操作符+:矩阵+矩阵,矩阵+标量
- 操作符+=:矩阵+=矩阵,矩阵+=标量
- 操作符 *:矩阵 *标量,矩阵 * 矩阵
- 操作符*=:矩阵*=标量
- 求方阵行列式的函数determinant
写好的程序如下:
myMatrix.h 头文件
// 2020-06-17 class myMatrix - head file
// a general matrix class
// Author: Guanhua Mei
#ifndef myMatrix_H
#define myMatrix_H
#include <iostream>
using std::ostream;
using std::istream;
using std::cout;
using std::endl;
class myMatrix
{
private:
unsigned int m_nRow; // number of rows
unsigned int m_nCol; // number of columns
unsigned int m_nCmp; // number of total components(elements)
double * m_pdMatrix; // double pointer to allocate and free memory
public:
// constructors
// construct with number of rows and columns, and initial value, default 0
myMatrix(unsigned int nRow = 0, unsigned int nCol = 0, double dbl = 0);
// construct with another matrix, copy constructor
myMatrix(const myMatrix &);
// overload =, assignment function
myMatrix & operator=(const myMatrix &);
// destructor
~myMatrix();
// judge if a valid index, i.e., i < m_nRow and j < m_nCol
bool ifValidIdx(unsigned int, unsigned int) const;
// get reference to Aij use A(i, j), note Aij can be modified.
double & operator()(unsigned int, unsigned int);
// get number of Rows, Columns and Components
unsigned int getNumRow() const {return m_nRow;}
unsigned int getNumCol() const {return m_nCol;}
unsigned int getNumCmp() const {return m_nCmp;}
// get Aij, set Aij
double getAij(unsigned int, unsigned int) const;
void setAij(unsigned int, unsigned int, double);
// ostream <<, output matrix to screen
friend ostream & operator<<(ostream &, const myMatrix &);
// istream >>, read matrix by keyboard, row by row
friend istream & operator>>(istream &, myMatrix &);
// operator +, matrix + scalar, matrix + matrix
friend myMatrix operator+(const myMatrix &, double);
friend myMatrix operator+(const myMatrix &, const myMatrix &);
// operator +=, matrix += scalar, matrix += matrix
myMatrix & operator+=(double);
myMatrix & operator+=(const myMatrix &);
// operator -, matrix - scalar, matrix - matrix
// opeartor -=, matrix -= scalar, matrix -= matrix
// similar as +, +=
// operator *, matrix * scalar, matrix * matrix
friend myMatrix operator*(const myMatrix &, double);
friend myMatrix operator*(const myMatrix &, const myMatrix &);
// operator *=, matrix *= scalar
myMatrix & operator*=(double);
// determinant
friend double determinant(const myMatrix &);
// transpose matrix
// inverse matrix
// eigenvectors and eigenvalues
// solve Ax = b, Gauss Elimination, Gauss-Seidel, CG, PCG, BiCG, AMG, ...
// to be continued...
};
#endif
myMatrix.cpp文件
// 2020-06-17 class myMatrix - cpp file
#include "myMatrix.h"
// construct with input row and col number, default 0
myMatrix::myMatrix(unsigned int nRow, unsigned int nCol, double dbl)
{
m_nRow = nRow; // set matrix row dimension
m_nCol = nCol; // set matrix column dimension
m_nCmp = nRow * nCol; // set matrix components number
m_pdMatrix = NULL;
// allocate memory and initialize with given value
if (m_nCmp > 0)
{
m_pdMatrix = new double [m_nCmp];
for (unsigned i = 0; i < m_nCmp; i++)
m_pdMatrix[i] = dbl;
}
}
// construct with another matrix, copy constructor
myMatrix::myMatrix(const myMatrix & mtx)
{
m_nRow = mtx.m_nRow; // set matrix row dimension
m_nCol = mtx.m_nCol; // set matrix column dimension
m_nCmp = mtx.m_nCmp; // set matrix components number
m_pdMatrix = NULL;
// allocate memory and assign values
if (m_nCmp > 0)
{
m_pdMatrix = new double [m_nCmp];
for (unsigned i = 0; i < m_nCmp; i++)
m_pdMatrix[i] = mtx.m_pdMatrix[i];
}
}
// overload =, assignement function
myMatrix & myMatrix::operator=(const myMatrix & mtx)
{
if (this == & mtx) // prevent self copy
return *this;
m_nRow = mtx.m_nRow; // set matrix row dimension
m_nCol = mtx.m_nCol; // set matrix column dimension
m_nCmp = mtx.m_nCmp; // set matrix components number
// free old memory
if (m_pdMatrix != NULL)
delete [] m_pdMatrix;
m_pdMatrix = NULL;
// allocate new memory and assign values
if (m_nCmp > 0)
{
m_pdMatrix = new double [m_nCmp];
for (unsigned i = 0; i < m_nCmp; i++)
m_pdMatrix[i] = mtx.m_pdMatrix[i];
}
return *this;
}
// destructor
myMatrix::~myMatrix()
{
// free memory
if (m_pdMatrix != NULL)
delete [] m_pdMatrix;
m_pdMatrix = NULL;
}
// judge if i < m_nRow and j < m_nCol
bool myMatrix::ifValidIdx(unsigned int i, unsigned int j) const
{
if ((i < m_nRow) && (j < m_nCol))
return true;
else
{
cout << "Invalid idinces! Error!" << endl;
return false;
}
}
// get reference to Aij use A(i, j), note Aij can be modified.
double & myMatrix::operator()(unsigned int i, unsigned int j)
{
if (ifValidIdx(i, j))
return m_pdMatrix[i * m_nCol + j];
else
exit(1);
}
// get Aij
double myMatrix::getAij(unsigned int i, unsigned int j) const
{
if (ifValidIdx(i, j))
return m_pdMatrix[i * m_nCol + j];
else
exit(1);
}
// set Aij
void myMatrix::setAij(unsigned int i, unsigned int j, double dbl)
{
if (ifValidIdx(i, j))
m_pdMatrix[i * m_nCol + j] = dbl;
}
// ostream <<
ostream & operator<<(ostream & os, const myMatrix & mtx)
{
if (mtx.m_nCmp <= 0)
{
os << "Empty matrix! " << endl;
return os;
}
for (unsigned int i = 0; i < mtx.m_nRow; i++)
{
for (unsigned int j = 0; j < mtx.m_nCol; j++)
os << mtx.getAij(i, j) << " ";
os << endl;
}
return os;
}
// istream >>
istream & operator>>(istream & is, myMatrix & mtx)
{
for (unsigned i = 0; i < mtx.m_nRow; i++)
for (unsigned j = 0; j < mtx.m_nCol; j++)
is >> mtx(i, j);
return is;
}
// operator +, matrix + scalar, matrix + matrix
myMatrix operator+(const myMatrix & mtx, double dbl)
{
myMatrix mtxTmp(mtx);
for (unsigned int i = 0; i < mtxTmp.m_nCmp; i++)
mtxTmp.m_pdMatrix[i] += dbl;
return mtxTmp;
}
myMatrix operator+(const myMatrix & mtx1, const myMatrix & mtx2)
{
if ((mtx1.m_nRow != mtx2.m_nRow) || (mtx1.m_nCol != mtx2.m_nCol))
{
cout << "Two matrices dimensions do not match, cannot add! " << endl;
exit(1);
}
else
{
myMatrix mtxTmp(mtx1);
for (unsigned int i = 0; i < mtxTmp.m_nCmp; i++)
mtxTmp.m_pdMatrix[i] += mtx2.m_pdMatrix[i];
return mtxTmp;
}
}
// operator +=, matrix += scalar, matrix += matrix
myMatrix & myMatrix::operator+=(double dbl)
{
for (unsigned int i = 0; i < m_nCmp; i++)
m_pdMatrix[i] += dbl;
return *this;
}
myMatrix & myMatrix::operator+=(const myMatrix & mtx)
{
if ((m_nRow != mtx.m_nRow) || (m_nCol != mtx.m_nCol))
{
cout << "Two matrices dimensions do not match, add 0! " << endl;
return *this;
}
else
{
for (unsigned int i = 0; i < m_nCmp; i++)
m_pdMatrix[i] += mtx.m_pdMatrix[i];
return *this;
}
}
// operator *, matrix * scalar, matrix * matrix
myMatrix operator*(const myMatrix & mtx, double dbl)
{
myMatrix mtxTmp(mtx);
for (unsigned int i = 0; i < mtxTmp.m_nCmp; i++)
mtxTmp.m_pdMatrix[i] *= dbl;
return mtxTmp;
}
myMatrix operator*(const myMatrix & mtxLft, const myMatrix & mtxRgt)
{
if (mtxLft.getNumCol() != mtxRgt.getNumRow())
{
cout << "Left matrix column dimension do not equals right matrix row "
<< "dimension, cannot *, return an empty matrix! " << endl;
myMatrix mtxTmp;
return mtxTmp;
}
else
{
unsigned int i = mtxLft.getNumRow();
unsigned int j = mtxLft.getNumCol();
unsigned int k = mtxRgt.getNumCol();
myMatrix mtxTmp(i, k, 0);
double dSum(0); // sum(matrixLeft[i,k] * matrixRight[k,j], k=0,1,...)
for (unsigned int iRow = 0; iRow < i; iRow++)
for (unsigned int iCol = 0; iCol < k; iCol++)
{
dSum = 0;
for (unsigned int iMid = 0; iMid < j; iMid++)
{
dSum += mtxLft.getAij(iRow, iMid) * mtxRgt.getAij(iMid, iCol);
}
mtxTmp(iRow, iCol) = dSum;
}
return mtxTmp;
}
}
// operator *=, matrix *= scalar
myMatrix & myMatrix::operator*=(double dbl)
{
for (unsigned int i = 0; i < m_nCmp; i++)
m_pdMatrix[i] *= dbl;
return *this;
}
// determinant
double determinant(const myMatrix & mtx)
{
// if not a square matrix, can not calculate determinant
if (mtx.m_nRow != mtx.m_nCol)
{
cout << "Not a square matrix, no determinant. " << endl;
exit(1);
}
// if a 0*0 matrix
if (mtx.m_nRow == 0)
return 0;
// if a 1*1 matrix
if (mtx.m_nRow == 1)
return mtx.getAij(0, 0);
// if a 2*2 matrix
if (mtx.m_nRow == 2)
return (mtx.getAij(0, 0)*mtx.getAij(1, 1)
- mtx.getAij(0, 1)*mtx.getAij(1, 0));
// if a 3*3 matrix
if (mtx.m_nRow == 3)
{
return mtx.getAij(0, 0) *
(mtx.getAij(1, 1) * mtx.getAij(2, 2) -
mtx.getAij(2, 1) * mtx.getAij(1, 2))
- mtx.getAij(0, 1) *
(mtx.getAij(1, 0) * mtx.getAij(2, 2) -
mtx.getAij(2, 0) * mtx.getAij(1, 2)) +
+ mtx.getAij(0, 2) *
(mtx.getAij(1, 0) * mtx.getAij(2, 1) -
mtx.getAij(2, 0) * mtx.getAij(1, 1));
}
// if a n*n matrix, n > 3
// recursion to calculate determinant
double det(0);
for (unsigned j = 0; j < mtx.m_nCol; j++) // the jth column to be deleted
{
// a matrix do not contain first row and Jth col
myMatrix mtxTmp(mtx.m_nRow - 1, mtx.m_nCol - 1, 0.0);
if (j == 0) // if delete the 1st column
{
for (unsigned iRow = 0; iRow < mtx.m_nRow - 1; iRow++)
for (unsigned jCol = 0; jCol < mtx.m_nCol - 1; jCol++)
mtxTmp(iRow, jCol) = mtx.getAij(iRow + 1, jCol + 1);
}
if (j == mtx.m_nCol - 1) // if delete the last column
{
for (unsigned iRow = 0; iRow < mtx.m_nRow - 1; iRow++)
for (unsigned jCol = 0; jCol < mtx.m_nCol - 1; jCol++)
mtxTmp(iRow, jCol) = mtx.getAij(iRow + 1, jCol);
}
if ((0 < j) && (j < mtx.m_nCol - 1)) //if delete the middle column
{
// 0, 1, ..., j-1 retained
for (unsigned iRow = 0; iRow < mtx.m_nRow - 1; iRow++)
for (unsigned jCol = 0; jCol <= j - 1; jCol++)
mtxTmp(iRow, jCol) = mtx.getAij(iRow + 1, jCol);
// j+1, j+2, ..., nCol - 1 retained
for (unsigned iRow = 0; iRow < mtx.m_nRow - 1; iRow++)
for (unsigned jCol = j + 1; jCol <= mtx.m_nCol - 1; jCol++)
mtxTmp(iRow, jCol - 1) = mtx.getAij(iRow + 1, jCol);
}
// cout << mtxTmp << endl << endl;
// add (-1)^(i+j)*A[i, j] * det(matrix(without iRow and jCol)) to dSum
if (((2 + j) % 2) == 0) // odd
det += (mtx.getAij(0, j) * determinant(mtxTmp));
else // even
det -= (mtx.getAij(0, j) * determinant(mtxTmp));
}
return det;
}
测试主函数main.cpp
// 2020-06-17 test class myMatrix
// Author: Guanhua Mei
#include <iostream>
#include "myMatrix.h"
using namespace std;
int main()
{
myMatrix matrix1; // construct with none
cout << "After myMatrix matrix1, matrix1: " << endl;
cout << matrix1 << endl; // test <<
myMatrix matrix2(2); // construct with nRow
cout << "After myMatrix matrix2(2), matrix2: " << endl;
cout << matrix2 << endl;
myMatrix matrix3(2, 3, 1); // construct with nRow, nCol, val
cout << "After myMatrix matrix3(2, 3, 1), matrix3: " << endl;
cout << matrix3 << endl;
myMatrix matrix4(3, 2, 2); // construct with nRow, nCol, val
cout << "After myMatrix matrix4(3, 2, 2), matrix4: " << endl;
cout << matrix4 << endl;
myMatrix matrix5(matrix3); // construct with matrix
cout << "After myMatrix matrix5(matrix3), matrix5: " << endl;
cout << matrix5 << endl;
matrix5 = matrix4; // test =
cout << "After matrix5 = matrix4, matrix5: " << endl;
cout << matrix5 << endl;
matrix5.setAij(0, 0, 1); // test set
matrix5.setAij(0, 1, 2);
matrix5(1, 0) = 3; // test ()
matrix5(1, 1) = 4;
matrix5(2, 0) = 5;
matrix5(2, 1) = 6;
cout << "matrix5.setAij=X and matrix5(i,j)=X:" << endl;
cout << matrix5 << endl;
cout << "matrix5.getAij(1, 1):"; // test setAij
cout << matrix5.getAij(1, 1) << endl;
cout << "matrix5(1, 1):";
cout << matrix5(1, 1) << endl;
cout << "matrix5 + matrix5: " << endl;
cout << matrix5 + matrix5 << endl;
cout << "matrix5 + 1:" << endl;
cout << matrix5 + 1 << endl;
cout << "matrix5 * 2:" << endl;
cout << matrix5 * 2 << endl;
cout << "matrix5 * matrix3:" << endl;
cout << matrix5 * matrix3 << endl;
cout << "matrix5 += 1:" << endl;
matrix5 += 1;
cout << matrix5 << endl;
cout << "matrix5 *= 2:" << endl;
matrix5 *= 2;
cout << matrix5 << endl;
// test determinant
myMatrix matrix0(0, 0, 0);
cout << "det(emptyMatrix): " << endl;
cout << determinant(matrix0) << endl << endl;
myMatrix mtx1(1, 1, 1);
cout << "det(1): " << endl;
cout << determinant(mtx1) << endl << endl;
myMatrix mtx2(2, 2, 0);
mtx2(0, 0) = 1;
mtx2(0, 1) = 2;
mtx2(1, 0) = 3;
mtx2(1, 1) = 4;
cout << "det of matrix " << endl << mtx2 << "is: ";
cout << determinant(mtx2) << endl << endl;
myMatrix mtx3(3, 3, 0);
mtx3(0, 0) = 3;
mtx3(0, 1) = 1;
mtx3(0, 2) = 4;
mtx3(1, 0) = 5;
mtx3(1, 1) = 9;
mtx3(1, 2) = 2;
mtx3(2, 0) = 6;
mtx3(2, 1) = 8;
mtx3(2, 2) = 7;
cout << "det of matrix " << endl << mtx3 << "is: ";
cout << determinant(mtx3) << endl << endl;
myMatrix mtx4(4, 4, 0);
cout << "Please input a mtrix by keyboard, row by row. " << endl;
cin >> mtx4; // test >>
// type numbers by keyboard: 3 1 4 2 5 9 2 3 6 8 7 5 4 9 7 6
cout << "det of matrix " << endl << mtx4 << "is: ";
cout << determinant(mtx4) << endl << endl; // 45, correct!
return 0;
}
测试结果,表明所有函数均实现了其功能。