Click the blue text above to subscribe!
In scientific computing and engineering applications, matrix operations are fundamental and critical operations. This article will introduce a complete implementation of a C++ matrix class that supports common matrix operations, including addition, multiplication, transposition, determinant calculation, and inverse matrix solving. We will delve into design details, implementation techniques, and performance optimizations.
Mathematical Foundations of the Matrix Class
A matrix is a core concept in linear algebra, and a matrix can be represented as:
Matrix operations follow specific algebraic rules:
- Addition: where
- Multiplication: where
- Transposition: where
- Determinant: defined only for square matrices, denoted as
- Inverse Matrix: satisfies
Complete Implementation of the Matrix Class
Matrix.h Header File
#pragma once
#include<vector>
#include<stdexcept>
#include<cmath>
#include<algorithm>
#include<utility>
#include<iostream>
#include<iomanip>
#include<initializer_list>
class Matrix {
public:
// Constructor
Matrix(int rows = 0, int cols = 0);
Matrix(const std::vector<std::vector<double>>& data);
Matrix(std::initializer_list<std::initializer_list<double>> list);
// Five Rule Implementations
Matrix(const Matrix& other);
Matrix(Matrix&& other) noexcept;
Matrix& operator=(const Matrix& other);
Matrix& operator=(Matrix&& other) noexcept;
~Matrix() = default;
// Element Access Operators
double& operator()(int i, int j);
const double& operator()(int i, int j) const;
std::vector<double>& operator[](int i);
const std::vector<double>& operator[](int i) const;
// Matrix Operations
Matrix operator+(const Matrix& other) const;
Matrix operator-(const Matrix& other) const;
Matrix operator*(const Matrix& other) const;
Matrix operator*(double scalar) const;
friend Matrix operator*(double scalar, const Matrix& matrix);
// Special Operations
operator std::vector<std::vector<double>>() const { return _data; }
Matrix transpose() const;
double determinant() const;
Matrix inverse() const;
// Utility Functions
int rows() const { return _rows; }
int cols() const { return _cols; }
bool isSquare() const { return _rows == _cols; }
void resize(int rows, int cols);
void zero();
void identity();
double norm() const;
void print(const std::string& title = "") const;
// Static Utility Functions
static Matrix zeros(int rows, int cols);
static Matrix ones(int rows, int cols);
static Matrix eye(int n);
private:
std::vector<std::vector<double>> _data;
int _rows, _cols;
// Determinant Calculation Helper Function
double determinantHelper(const Matrix& mat) const;
};
// Stream Output Operator
std::ostream& operator<<(std::ostream& os, const Matrix& matrix);
Matrix.cpp Implementation File
#include "Matrix.h"
#include<stdexcept>
#include<cmath>
#include<iomanip>
#include<string>
// Constructor: Specify rows and columns
Matrix::Matrix(int rows, int cols) : _rows(rows), _cols(cols) {
if (rows < 0 || cols < 0) {
throw std::invalid_argument("Matrix dimensions must be non-negative");
}
_data.resize(static_cast<size_t>(_rows),
std::vector<double>(static_cast<size_t>(_cols), 0.0));
}
// Constructor: 2D Vector
Matrix::Matrix(const std::vector<std::vector<double>>& data) {
if (data.empty()) {
_rows = _cols = 0;
return;
}
_rows = static_cast<int>(data.size());
_cols = static_cast<int>(data[0].size());
for (int i = 0; i < _rows; ++i) {
if (static_cast<int>(data[i].size()) != _cols) {
throw std::invalid_argument("All rows of the matrix must have the same size");
}
}
_data = data;
}
// Constructor: Initializer List
Matrix::Matrix(std::initializer_list<std::initializer_list<double>> list) {
if (list.size() == 0) {
_rows = _cols = 0;
return;
}
_rows = static_cast<int>(list.size());
_cols = static_cast<int>(list.begin()->size());
_data.reserve(static_cast<size_t>(_rows));
for (const auto& rowList : list) {
if (static_cast<int>(rowList.size()) != _cols) {
throw std::invalid_argument("All rows of the matrix must have the same size");
}
_data.push_back(std::vector<double>(rowList));
}
}
// Copy Constructor
Matrix::Matrix(const Matrix& other)
: _rows(other._rows), _cols(other._cols), _data(other._data) {}
// Move Constructor
Matrix::Matrix(Matrix&& other) noexcept
: _rows(other._rows), _cols(other._cols), _data(std::move(other._data)) {
other._rows = 0;
other._cols = 0;
}
// Copy Assignment Operator
Matrix& Matrix::operator=(const Matrix& other) {
if (this != &other) {
_rows = other._rows;
_cols = other._cols;
_data = other._data;
}
return *this;
}
// Move Assignment Operator
Matrix& Matrix::operator=(Matrix&& other) noexcept {
if (this != &other) {
_rows = other._rows;
_cols = other._cols;
_data = std::move(other._data);
other._rows = 0;
other._cols = 0;
}
return *this;
}
// Element Access (non-const)
double& Matrix::operator()(int i, int j) {
if (i < 0 || i >= _rows || j < 0 || j >= _cols) {
throw std::out_of_range("Matrix index out of range");
}
return _data[static_cast<size_t>(i)][static_cast<size_t>(j)];
}
// Element Access (const)
const double& Matrix::operator()(int i, int j) const {
if (i < 0 || i >= _rows || j < 0 || j >= _cols) {
throw std::out_of_range("Matrix index out of range");
}
return _data[static_cast<size_t>(i)][static_cast<size_t>(j)];
}
// Row Access (non-const)
std::vector<double>& Matrix::operator[](int i) {
if (i < 0 || i >= _rows) {
throw std::out_of_range("Matrix row index out of range");
}
return _data[static_cast<size_t>(i)];
}
// Row Access (const)
const std::vector<double>& Matrix::operator[](int i) const {
if (i < 0 || i >= _rows) {
throw std::out_of_range("Matrix row index out of range");
}
return _data[static_cast<size_t>(i)];
}
// Matrix Addition
Matrix Matrix::operator+(const Matrix& other) const {
if (_rows != other._rows || _cols != other._cols) {
throw std::invalid_argument("Matrix dimensions must match for addition");
}
Matrix result(_rows, _cols);
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < _cols; ++j) {
result(i, j) = _data[static_cast<size_t>(i)][static_cast<size_t>(j)] + other(i, j);
}
}
return result;
}
// Matrix Subtraction
Matrix Matrix::operator-(const Matrix& other) const {
if (_rows != other._rows || _cols != other._cols) {
throw std::invalid_argument("Matrix dimensions must match for subtraction");
}
Matrix result(_rows, _cols);
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < _cols; ++j) {
result(i, j) = _data[static_cast<size_t>(i)][static_cast<size_t>(j)] - other(i, j);
}
}
return result;
}
// Matrix Multiplication
Matrix Matrix::operator*(const Matrix& other) const {
if (_cols != other._rows) {
throw std::invalid_argument(
"Matrix dimensions incompatible for multiplication: " +
std::to_string(_rows) + "x" + std::to_string(_cols) + " * " +
std::to_string(other._rows) + "x" + std::to_string(other._cols));
}
Matrix result(_rows, other._cols);
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < other._cols; ++j) {
double sum = 0.0;
for (int k = 0; k < _cols; ++k) {
sum += _data[static_cast<size_t>(i)][static_cast<size_t>(k)] * other(k, j);
}
result(i, j) = sum;
}
}
return result;
}
// Scalar Multiplication (Right)
Matrix Matrix::operator*(double scalar) const {
Matrix result(_rows, _cols);
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < _cols; ++j) {
result(i, j) = _data[static_cast<size_t>(i)][static_cast<size_t>(j)] * scalar;
}
}
return result;
}
// Scalar Multiplication (Left)
Matrix operator*(double scalar, const Matrix& matrix) {
return matrix * scalar;
}
// Matrix Transposition
Matrix Matrix::transpose() const {
Matrix result(_cols, _rows);
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < _cols; ++j) {
result(j, i) = _data[static_cast<size_t>(i)][static_cast<size_t>(j)];
}
}
return result;
}
// Determinant Calculation Helper Function (Recursive)
double Matrix::determinantHelper(const Matrix& mat) const {
const int n = mat.rows();
if (n == 1) return mat(0, 0);
if (n == 2) return mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0);
double det = 0.0;
int sign = 1;
for (int j = 0; j < n; ++j) {
// Create Minor Matrix
Matrix minor(n-1, n-1);
for (int r = 1; r < n; ++r) {
int colIdx = 0;
for (int c = 0; c < n; ++c) {
if (c == j) continue;
minor(r-1, colIdx++) = mat(r, c);
}
}
// Add signed minor determinant
det += sign * mat(0, j) * minor.determinant();
sign = -sign;
}
return det;
}
// Determinant Calculation
double Matrix::determinant() const {
if (!isSquare()) {
throw std::invalid_argument("Determinant is only defined for square matrices");
}
// Special matrix handling
if (_rows == 0) return 1.0; // Determinant of 0x0 matrix is 1 (by convention)
if (_rows == 1) return _data[0][0];
if (_rows == 2) return _data[0][0] * _data[1][1] - _data[0][1] * _data[1][0];
return determinantHelper(*this);
}
// Inverse Matrix Calculation
Matrix Matrix::inverse() const {
if (!isSquare()) {
throw std::invalid_argument("Inverse is only defined for square matrices");
}
const int n = _rows;
// Special matrix handling
if (n == 0) return Matrix(0, 0);
if (n == 1) {
Matrix result(1, 1);
result(0, 0) = 1.0 / _data[0][0];
return result;
}
const double det = determinant();
if (std::abs(det) < 1e-12) {
throw std::runtime_error("Matrix is singular (non-invertible)");
}
// Calculate Cofactor Matrix
Matrix cofactors(n, n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
// Create Minor Matrix
Matrix minor(n-1, n-1);
for (int r = 0, ri = 0; r < n; ++r) {
if (r == i) continue;
for (int c = 0, cj = 0; c < n; ++c) {
if (c == j) continue;
minor(ri, cj++) = _data[static_cast<size_t>(r)][static_cast<size_t>(c)];
}
ri++;
}
// Algebraic Cofactor = Sign * Minor Determinant
double cofactor = minor.determinant();
int sign = ((i+j) % 2 == 0) ? 1 : -1;
cofactors(i, j) = sign * cofactor;
}
}
// Inverse Matrix = (1/Determinant) * Transpose of Cofactor Matrix
return cofactors.transpose() * (1.0 / det);
}
// Resize Matrix
void Matrix::resize(int rows, int cols) {
if (rows < 0 || cols < 0) {
throw std::invalid_argument("Matrix dimensions must be non-negative");
}
_rows = rows;
_cols = cols;
_data.resize(static_cast<size_t>(_rows));
for (auto& row : _data) {
row.resize(static_cast<size_t>(_cols), 0.0);
}
}
// Zero Matrix
void Matrix::zero() {
for (auto& row : _data) {
std::fill(row.begin(), row.end(), 0.0);
}
}
// Identity Matrix
void Matrix::identity() {
if (!isSquare()) {
throw std::invalid_argument("Identity matrix must be square");
}
zero();
for (int i = 0; i < _rows; ++i) {
_data[static_cast<size_t>(i)][static_cast<size_t>(i)] = 1.0;
}
}
// Matrix Norm (Frobenius Norm)
double Matrix::norm() const {
double sum = 0.0;
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < _cols; ++j) {
double val = _data[static_cast<size_t>(i)][static_cast<size_t>(j)];
sum += val * val;
}
}
return std::sqrt(sum);
}
// Static Function: Zero Matrix
Matrix Matrix::zeros(int rows, int cols) {
return Matrix(rows, cols);
}
// Static Function: Ones Matrix
Matrix Matrix::ones(int rows, int cols) {
Matrix result(rows, cols);
for (int i = 0; i < rows; ++i) {
for (int j = 0; j < cols; ++j) {
result(i, j) = 1.0;
}
}
return result;
}
// Static Function: Identity Matrix
Matrix Matrix::eye(int n) {
Matrix result(n, n);
for (int i = 0; i < n; ++i) {
result(i, i) = 1.0;
}
return result;
}
// Matrix Print
void Matrix::print(const std::string& title) const {
if (!title.empty()) {
std::cout << title << ":\n";
}
for (int i = 0; i < _rows; ++i) {
for (int j = 0; j < _cols; ++j) {
std::cout << std::setw(10) << std::setprecision(4)
<< _data[static_cast<size_t>(i)][static_cast<size_t>(j)] << " ";
}
std::cout << "\n";
}
}
// Stream Output Operator
std::ostream& operator<<(std::ostream& os, const Matrix& matrix) {
for (int i = 0; i < matrix.rows(); ++i) {
for (int j = 0; j < matrix.cols(); ++j) {
os << std::setw(10) << std::setprecision(4) << matrix(i, j) << " ";
}
os << "\n";
}
return os;
}
Testing
1. Test Main Function
#include "Matrix.h"
#include<iostream>
#include<cmath>
#include<iomanip>
void testConstructors() {
std::cout << "==== Testing Constructors ====" << \n";
// Default constructor
Matrix m1;
std::cout << "Default constructor (0x0):\n" << m1 << std::endl;
// Size constructor
Matrix m2(3, 4);
std::cout << "Size constructor (3x4):\n" << m2 << std::endl;
// Initializer list constructor
Matrix m3 = { {1, 2, 3}, {4, 5, 6}, {7, 8, 9} };
std::cout << "Initializer list constructor:\n" << m3 << std::endl;
// Vector constructor
std::vector<std::vector<double>> data = {
{1.5, 2.5},
{3.5, 4.5}
};
Matrix m4(data);
std::cout << "Vector constructor:\n" << m4 << std::endl;
// Copy constructor
Matrix m5 = m3;
std::cout << "Copy constructor:\n" << m5 << std::endl;
// Move constructor
Matrix m6 = Matrix(2, 2);
m6(0, 0) = 1; m6(0, 1) = 2;
m6(1, 0) = 3; m6(1, 1) = 4;
std::cout << "Move constructor (2x2):\n" << m6 << std::endl;
}
void testAccessors() {
std::cout << "\n==== Testing Accessors ====" << \n";
Matrix m(3, 3);
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
m(i, j) = i * 3 + j + 1;
}
}
std::cout << "Original matrix:\n" << m << std::endl;
// Test element access
std::cout << "Element (1,1): " << m(1, 1) << std::endl;
m(1, 1) = 99;
std::cout << "After setting (1,1) to 99:\n" << m << std::endl;
// Test row access
std::vector<double> row = m[0];
std::cout << "Row 0: ";
for (double val : row) std::cout << val << " ";
std::cout << std::endl;
// Test const access
const Matrix& cm = m;
std::cout << "Const element (0,0): " << cm(0, 0) << std::endl;
// Test out-of-bound access (should throw)
try {
std::cout << "Trying to access (3,3): ";
std::cout << m(3, 3) << std::endl;
}
catch (const std::out_of_range& e) {
std::cout << "Caught exception: " << e.what() << std::endl;
}
}
void testOperations() {
std::cout << "\n==== Testing Operations ====" << \n";
Matrix A = { {1, 2}, {3, 4} };
Matrix B = { {5, 6}, {7, 8} };
std::cout << "Matrix A:\n" << A;
std::cout << "Matrix B:\n" << B;
// Addition
Matrix C = A + B;
std::cout << "A + B:\n" << C;
// Subtraction
Matrix D = A - B;
std::cout << "A - B:\n" << D;
// Multiplication
Matrix E = A * B;
std::cout << "A * B:\n" << E;
// Scalar multiplication
Matrix F = A * 2.5;
std::cout << "A * 2.5:\n" << F;
// Left scalar multiplication
Matrix G = 3.0 * B;
std::cout << "3.0 * B:\n" << G;
// Transpose
Matrix H = A.transpose();
std::cout << "Transpose of A:\n" << H;
}
void testAdvancedOperations() {
std::cout << "\n==== Testing Advanced Operations ====" << \n";
Matrix A = { {4, 7}, {2, 6} };
std::cout << "Matrix A:\n" << A;
// Determinant
double det = A.determinant();
std::cout << "Determinant of A: " << det << " (expected: 10)\n";
// Inverse
Matrix inv = A.inverse();
std::cout << "Inverse of A:\n" << inv;
// Verify A * inv ≈ identity
Matrix identity = A * inv;
identity.print("A * inv (should be identity):");
for (int i = 0; i < identity.rows(); i++) {
for (int j = 0; j < identity.cols(); j++) {
double val = identity(i, j);
if (i == j) {
if (std::abs(val - 1.0) > 1e-6) {
std::cout << "Warning: diagonal element " << i << "," << j
<< " = " << val << " (expected 1.0)" << std::endl;
}
}
else {
if (std::abs(val) > 1e-6) {
std::cout << "Warning: off-diagonal element " << i << "," << j
<< " = " << val << " (expected 0.0)" << std::endl;
}
}
}
}
// Test with singular matrix
Matrix singular = { {1, 2}, {2, 4} };
std::cout << "Singular matrix:\n" << singular;
try {
double detSingular = singular.determinant();
std::cout << "Determinant of singular matrix: " << detSingular << " (expected: 0)\n";
Matrix invSingular = singular.inverse();
}
catch (const std::exception& e) {
std::cout << "Caught expected exception: " << e.what() << std::endl;
}
}
void testUtilityFunctions() {
std::cout << "\n==== Testing Utility Functions ====" << \n";
// Resize
Matrix m(2, 3);
m(0, 0) = 1; m(0, 1) = 2; m(0, 2) = 3;
m(1, 0) = 4; m(1, 1) = 5; m(1, 2) = 6;
std::cout << "Original matrix (2x3):\n" << m;
m.resize(3, 2);
std::cout << "After resize to 3x2:\n" << m;
// Zero
m.zero();
std::cout << "After zero():\n" << m;
// Identity
m.resize(3, 3);
m.identity();
std::cout << "Identity matrix:\n" << m;
// Norm
Matrix n = { {1, 2}, {3, 4} };
double norm = n.norm();
std::cout << "Norm of [[1,2],[3,4]]: " << norm
<< " (expected: " << std::sqrt(1 + 4 + 9 + 16) << ")\n";
}
void testStaticFunctions() {
std::cout << "\n==== Testing Static Functions ====" << \n";
// Zeros
Matrix z = Matrix::zeros(3, 2);
std::cout << "Zeros (3x2):\n" << z;
// Ones
Matrix o = Matrix::ones(2, 3);
std::cout << "Ones (2x3):\n" << o;
// Eye
Matrix e = Matrix::eye(4);
std::cout << "Eye (4x4):\n" << e;
}
void testLargeMatrixOperations() {
std::cout << "\n==== Testing Large Matrix Operations ====" << \n";
// Create 4x4 matrix instead of 100x100 for practical testing
const int N = 4;
Matrix large = Matrix::zeros(N, N);
// Fill with identity pattern
for (int i = 0; i < N; i++) {
large(i, i) = 2.0; // Diagonal
if (i > 0) large(i, i - 1) = -1.0; // Lower diagonal
if (i < N - 1) large(i, i + 1) = -1.0; // Upper diagonal
}
std::cout << "Large matrix (4x4):\n" << large;
// Test determinant
try {
double det = large.determinant();
std::cout << "Determinant of 4x4 matrix: " << det << std::endl;
}
catch (const std::exception& e) {
std::cout << "Error calculating determinant: " << e.what() << std::endl;
}
// Test inverse
try {
Matrix inv = large.inverse();
std::cout << "Successfully computed inverse of 4x4 matrix\n";
// Verify product is approximately identity
Matrix product = large * inv;
product.print("Product (should be identity):");
}
catch (const std::exception& e) {
std::cout << "Error computing inverse: " << e.what() << std::endl;
}
}
int main() {
try {
testConstructors();
testAccessors();
testOperations();
testAdvancedOperations();
testUtilityFunctions();
testStaticFunctions();
testLargeMatrixOperations();
}
catch (const std::exception& e) {
std::cerr << "Uncaught exception: " << e.what() << std::endl;
return 1;
}
std::cout << "\nAll tests completed successfully!" << std::endl;
return 0;
}
2. Test Results
==== Testing Constructors ====
Default constructor (0x0):
Size constructor (3x4):
0 0 0 0
0 0 0 0
0 0 0 0
Initializer list constructor:
1 2 3
4 5 6
7 8 9
Vector constructor:
1.5 2.5
3.5 4.5
Copy constructor:
1 2 3
4 5 6
7 8 9
Move constructor (2x2):
1 2
3 4
==== Testing Accessors ====
Original matrix:
1 2 3
4 5 6
7 8 9
Element (1,1): 5
After setting (1,1) to 99:
1 2 3
4 99 6
7 8 9
Row 0: 1 2 3
Const element (0,0): 1
Trying to access (3,3): Caught exception: Matrix indices out of range
==== Testing Operations ====
Matrix A:
1 2
3 4
Matrix B:
5 6
7 8
A + B:
6 8
10 12
A - B:
-4 -4
-4 -4
A * B:
19 22
43 50
A * 2.5:
2.5 5
7.5 10
3.0 * B:
15 18
21 24
Transpose of A:
1 3
2 4
==== Testing Advanced Operations ====
Matrix A:
4 7
2 6
Determinant of A: 10 (expected: 10)
Inverse of A:
0.6 -0.7
-0.2 0.4
A * inv (should be identity)::
1 0
0 1
Singular matrix:
1 2
2 4
Determinant of singular matrix: 0 (expected: 0)
Caught expected exception: Matrix is singular (non-invertible)
==== Testing Utility Functions ====
Original matrix (2x3):
1 2 3
4 5 6
After resize to 3x2:
1 2
4 5
0 0
After zero():
0 0
0 0
0 0
Identity matrix:
1 0 0
0 1 0
0 0 1
Norm of [[1,2],[3,4]]: 5.477 (expected: 5.477)
==== Testing Static Functions ====
Zeros (3x2):
0 0
0 0
0 0
Ones (2x3):
1 1 1
1 1 1
Eye (4x4):
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
==== Testing Large Matrix Operations ====
Large matrix (4x4):
2 -1 0 0
-1 2 -1 0
0 -1 2 -1
0 0 -1 2
Determinant of 4x4 matrix: 5
Successfully computed inverse of 4x4 matrix
Product (should be identity)::
1 0 0 0
1.11e-16 1 0 0
-5.551e-17 -1.11e-16 1 1.11e-16
0 0 0 1
All tests completed successfully!
Key Design and Implementation Details
1. Robust Memory Management
The matrix class uses<span>std::vector<std::vector<double>></span> as the underlying storage, automatically managing memory. We have implemented the completeFive Rule (copy constructor, move constructor, copy assignment, move assignment, and destructor) to ensure proper resource management:
// Move constructor example
Matrix::Matrix(Matrix&& other) noexcept
: _rows(other._rows), _cols(other._cols), _data(std::move(other._data)) {
other._rows = 0;
other._cols = 0;
}
2. Multi-dimensional Access Interface
Provides two ways to access elements:
<span>operator()</span>for row-column index access<span>operator[]</span>for row vector access
// Element access example
double& Matrix::operator()(int i, int j) {
if (i < 0 || i >= _rows || j < 0 || j >= _cols) {
throw std::out_of_range("Matrix index out of range");
}
return _data[static_cast<size_t>(i)][static_cast<size_t>(j)];
}
3. Determinant Calculation Optimization
The determinant calculation uses a recursive algorithm, with a time complexity of. It is efficient for small matrices, but for large matrices, consider more efficient algorithms like LU decomposition:
double Matrix::determinantHelper(const Matrix& mat) const {
const int n = mat.rows();
// Base case handling
if (n == 1) return mat(0, 0);
if (n == 2) return mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0);
// Recursive determinant calculation
double det = 0.0;
int sign = 1;
for (int j = 0; j < n; ++j) {
// Create minor matrix
Matrix minor(n-1, n-1);
// ... (fill minor)
det += sign * mat(0, j) * minor.determinant();
sign = -sign;
}
return det;
}
4. Numerical Stability of Inverse Matrix
The inverse matrix calculation uses the adjugate matrix method and includes a singularity check:
Matrix Matrix::inverse() const {
// ... (check for square matrix)
const double det = determinant();
if (std::abs(det) < 1e-12) {
throw std::runtime_error("Matrix is singular (non-invertible)");
}
// ... (calculate adjugate matrix)
return cofactors.transpose() * (1.0 / det);
}
5. Type Safety and Boundary Checks
All operations include dimension checks and boundary checks to prevent illegal operations:
Matrix Matrix::operator*(const Matrix& other) const {
if (_cols != other._rows) {
throw std::invalid_argument(
"Matrix dimensions incompatible for multiplication: " +
std::to_string(_rows) + "x" + std::to_string(_cols) + " * " +
std::to_string(other._rows) + "x" + std::to_string(other._cols));
}
// ... (multiplication implementation)
}
Performance Analysis and Optimization Directions
-
Storage Efficiency: Currently using nested vectors for storage, which may lead to non-contiguous memory. An optimization direction is to use a single contiguous memory block to store matrix elements.
-
Algorithm Optimization:
- Determinant calculation: Use LU decomposition instead of the recursive algorithm
- Inverse matrix calculation: Use Gaussian-Jordan elimination instead of the adjugate matrix method
Parallel Computing: Utilize modern CPU’s SIMD instructions and parallel capabilities to accelerate matrix operations
Expression Templates: Implement expression templates to avoid temporary matrix creation, optimizing performance for complex expressions
Conclusion
This article implements a fully functional and robust C++ matrix class that supports basic matrix operations and advanced linear algebra operations. Through reasonable interface design and robust error handling, this class is suitable for both educational purposes and practical scientific computing applications. For scenarios with higher performance requirements, algorithm optimizations and storage improvements can be made based on this implementation.
Matrix operations are ubiquitous in scientific computing, machine learning, and engineering simulations, and an efficient and reliable matrix implementation is a foundational tool in these fields. This implementation balances functionality, code clarity, and performance, laying a solid foundation for more complex numerical computation tasks.
• end • 
Companionship is the longest confession of love
We push the most practical information for you
Scan the QR code to follow us