Matrix.cpp

#include "Matrix.h"
#include <stdexcept>
#include <cmath>
#include <algorithm>

//构造函数和析构函数
Matrix::Matrix(): rows(0), cols(0) {}

Matrix::Matrix(int rows, int cols) : Matrix(rows, cols, 0.0) {}

Matrix::Matrix(int rows, int cols, double initValue) : rows(rows), cols(cols) {
    if (rows <= 0 || cols <= 0) {
        throw std::invalid_argument("Matrix dimensions must be positive");
    }
    data.resize(rows, std::vector<double>(cols, initValue));
}

Matrix::Matrix(const std::vector<std::vector<double>>& values) {
    if (values.empty() || values[0].empty()) {
        rows = 0; cols = 0; return;
    }
    rows = int(values.size());
    cols = int(values[0].size());
    for (int i = 1; i < rows; i++) {
        if (values[i].size() != cols) {
            throw std::invalid_argument("All rows must have the same number of columns in a matrix");
        }
    }
    data = values;
}

Matrix::Matrix(const Matrix& other) : rows(other.rows), cols(other.cols), data(other.data) {}

Matrix::~Matrix() {}

//私有辅助方法
inline void Matrix::avoidEmpty() const {
    if (empty()) throw std::domain_error("The matrix is empty");
}

void Matrix::swapRows(int row1, int row2) {
    if (row1 < 0 || row1 >= rows || row2 < 0 || row2 >= rows) {
        throw std::out_of_range("Row index out of range");
    }
    if (row1 == row2) return;
    std::swap(data[row1], data[row2]);
}

void Matrix::scaleRow(int row, double scalar) {
    if (row < 0 || row >= rows) {
        throw std::out_of_range("Row index out of range");
    }
    for (int j = 0; j < cols; j++) {
        data[row][j] *= scalar;
    }
}

void Matrix::addScaledRow(int targetRow, int sourceRow, double scalar) {
    if (targetRow < 0 || targetRow >= rows || sourceRow < 0 || sourceRow >= rows) {
        throw std::out_of_range("Row index out of range");
    }
    for (int j = 0; j < cols; j++) {
        data[targetRow][j] += data[sourceRow][j] * scalar;
    }
}

//特解向量
Matrix Matrix::getSolutionVector(int len, const std::set<int>& fixedUnknowns) const {
    if (!isVector()) {
        throw std::domain_error("The matrix must be a vector");
    }
    Matrix result(len, 1);
    int row = 0;
    for (int i = 0; i < len; i++) {
        if (fixedUnknowns.find(i) != fixedUnknowns.end()) {
            if (row >= rows) {
                throw std::invalid_argument("The set of fixed unknowns is unreasonable");
            }
            result(i, 0) = data[row][0];
            row++;
            continue;
        }
        result(i, 0) = 0.0;
    }
    return result;
}

//解系向量
Matrix Matrix::getSolutionVector(int unknown, int len, const std::set<int>& fixedUnknowns) const {
    if (!isVector()) {
        throw std::domain_error("The matrix must be a vector");
    }
    if (unknown < 0 || unknown >= len) {
        throw std::invalid_argument("The column of unknown is out of range");
    }
    Matrix result(len, 1);
    int row = 0;
    for (int i = 0; i < len; i++) {
        if (fixedUnknowns.find(i) != fixedUnknowns.end()) {
            if (row >= rows || i == unknown) {
                throw std::invalid_argument("The set of fixed unknowns is unreasonable");
            }
            result(i, 0) = -data[row][0];
            row++;
            continue;
        }
        if (i == unknown) result(i, 0) = 1.0;
        else result(i, 0) = 0.0;
    }
    return result;
}

//矩阵操作
Matrix Matrix::getSubMatrix(int excludedRow, int excludedCol) const {
    if (excludedRow < 0 || excludedRow >= rows ||
        excludedCol < 0 || excludedCol >= cols) {
        throw std::out_of_range("Matrix index out of range");
    }
    if (rows == 1 && cols == 1) return Matrix{};
    Matrix sub(rows - 1, cols - 1);
    int subRow = 0;
    for (int i = 0; i < rows; i++) {
        if (i == excludedRow) continue;
        int subCol = 0;
        for (int j = 0; j < cols; j++) {
            if (j == excludedCol) continue;
            sub(subRow, subCol) = data[i][j];
            subCol++;
        }
        subRow++;
    }
    return sub;
}

Matrix Matrix::getSubMatrix(int startRow, int endRow, int startCol, int endCol) const {
    if (startRow < 0 || endRow >= rows || startRow > endRow ||
        startCol < 0 || endCol >= cols || startCol > endCol) {
        throw std::out_of_range("Matrix index out of range");
    }
    Matrix sub(endRow - startRow + 1, endCol - startCol + 1);
    int subRow = 0;
    for (int i = startRow; i <= endRow; i++) {
        int subCol = 0;
        for (int j = startCol; j <= endCol; j++) {
            sub(subRow, subCol) = data[i][j];
            subCol++;
        }
        subRow++;
    }
    return sub;
}

Matrix Matrix::combine(const Matrix& other) const {
    if (empty()) return other;
    if (rows != other.rows) {
        throw std::invalid_argument("The rows of matrices need to be identical");
    }
    Matrix result(rows, cols + other.cols);
    for (int i = 0; i < rows; i++) {
        int j = 0;
        for (int k = 0; k < cols; k++) {
            result(i, j++) = data[i][k];
        }
        for (int k = 0; k < other.cols; k++) {
            result(i, j++) = other.data[i][k];
        }
    }
    return result;
}

//基本访问器
int Matrix::getRows() const { return rows; }

int Matrix::getCols() const { return cols; }

double& Matrix::operator()(int row, int col) {
    if (row < 0 || row >= rows || col < 0 || col >= cols) {
        throw std::out_of_range("Matrix index out of range");
    }
    return data[row][col];
}

const double& Matrix::operator()(int row, int col) const {
    if (row < 0 || row >= rows || col < 0 || col >= cols) {
        throw std::out_of_range("Matrix index out of range");
    }
    return data[row][col];
}

Matrix::operator std::vector<std::vector<double>>() const {
    return data;
}

//赋值和算数运算
Matrix& Matrix::operator=(const Matrix& other) {
    if (this != &other) {
        rows = other.rows;
        cols = other.cols;
        data = other.data;
    }
    return *this;
}

Matrix Matrix::operator+(const Matrix& other) const {
    avoidEmpty();
    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[i][j] + other.data[i][j];
        }
    }
    return result;
}

Matrix Matrix::operator-(const Matrix& other) const {
    avoidEmpty();
    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[i][j] - other.data[i][j];
        }
    }
    return result;
}

Matrix Matrix::operator*(const Matrix& other) const {
    avoidEmpty();
    if (cols != other.rows) {
        throw std::invalid_argument("Matrix dimensions must match for multiplication");
    }
    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[i][k] * other.data[k][j];
            }
            result(i, j) = sum;
        }
    }
    return result;
}

Matrix Matrix::operator*(double scalar) const {
    avoidEmpty();
    Matrix result(rows, cols);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            result(i, j) = data[i][j] * scalar;
        }
    }
    return result;
}

Matrix operator*(double scalar, const Matrix& matrix) {
    return matrix * scalar;
}

//复合赋值运算
Matrix& Matrix::operator+=(const Matrix& other) {
    *this = *this + other;
    return *this;
}

Matrix& Matrix::operator-=(const Matrix& other) {
    *this = *this - other;
    return *this;
}

Matrix& Matrix::operator*=(const Matrix& other) {
    *this = *this * other;
    return *this;
}

Matrix& Matrix::operator*=(double scalar) {
    *this = *this * scalar;
    return *this;
}

// 比较运算
bool Matrix::operator==(const Matrix& other) const {
    avoidEmpty();
    if (rows != other.rows || cols != other.cols) {
        return false;
    }
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if (!equalsZero(data[i][j] - other.data[i][j])) {
                return false;
            }
        }
    }
    return true;
}

bool Matrix::operator!=(const Matrix& other) const {
    return !(*this == other);
}

//特殊矩阵
Matrix Matrix::identity(int size) {
    if (size <= 0) {
        throw std::invalid_argument("Invalid size for matrix");
    }
    Matrix result(size, size, 0.0);
    for (int i = 0; i < size; i++) {
        result(i, i) = 1.0;
    }
    return result;
}

Matrix Matrix::zeros(int rows, int cols) {
    if (rows <= 0 || cols <= 0) {
        throw std::invalid_argument("Invalid size for matrix");
    }
    Matrix result(rows, cols, 0.0);
    return result;
}

Matrix Matrix::ones(int rows, int cols) {
    if (rows <= 0 || cols <= 0) {
        throw std::invalid_argument("Invalid size for matrix");
    }
    Matrix result(rows, cols, 1.0);
    return result;
}

//行列式
double Matrix::determinant() const {
    if (!isSquare()) {
        throw std::domain_error("Determinant is only defined for square matrices");
    }
    int sign = 1;
    Matrix upperTri(*this);
    int topRow = 0;
    for (int j = 0; j < cols && topRow < rows; j++) {
        int pivotRow = -1;
        for (int i = topRow; i < rows; i++) {
            if (!equalsZero(upperTri(i, j))) {
                pivotRow = i;
                break;
            }
        }
        if (pivotRow == -1) continue;
        if (pivotRow != topRow) {
            upperTri.swapRows(topRow, pivotRow);
            sign = -sign;
        }
        for (int i = topRow + 1; i < rows; i++) {
            if (equalsZero(upperTri(i, j))) continue;
            double scalar = -upperTri(i, j) / upperTri(topRow, j);
            upperTri.addScaledRow(i, topRow, scalar);
        }
        topRow++;
    }
    double result = sign;
    for (int i = 0; i < rows; i++) {
        result *= upperTri.data[i][i];
    }
    return result;
}

//伴随矩阵
Matrix Matrix::adjugate() const {
    if (!isSquare()) {
        throw std::domain_error("Adjugate is only defined for square matrices");
    }
    Matrix result(rows, cols);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            Matrix sub(getSubMatrix(i, j));
            double algeCofactor = ((i + j) % 2 ? -1.0 : 1.0) * sub.determinant();
            result(j, i) = algeCofactor;
        }
    }
    return result;
}

//逆矩阵
Matrix Matrix::inverse() const {
    if (!isSquare()) {
        throw std::domain_error("Inverse is only defined for square matrices");
    }
    Matrix augmented(combine(identity(rows)));
    Matrix reduced = augmented.reducedEchelonForm();
    Matrix leftPart = reduced.getSubMatrix(0, rows - 1, 0, cols - 1);
    Matrix identityMatrix = identity(rows);
    if (leftPart != identityMatrix) {
        throw std::domain_error("Matrix is singular and cannot be inverted");
    }
    return reduced.getSubMatrix(0, rows - 1, cols, reduced.cols - 1);
}

//转置
Matrix Matrix::transpose() const {
    avoidEmpty();
    Matrix result(cols, rows);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) 
            result(j, i) = data[i][j];
    }
    return result;
}

//秩
int Matrix::rank() const {
    avoidEmpty();
    Matrix echelon(echelonForm());
    int r = 0;
    for (int i = 0; i < rows; i++) {
        bool isZeroRow = true;
        for (int j = 0; j < cols; j++) {
            if (!equalsZero(echelon(i, j))) {
                isZeroRow = false;
                break;
            }
        }
        if (isZeroRow) break;
        r++;
    }
    return r;
}

//方阵幂次
Matrix Matrix::pow(int exponent) const {
    if (!isSquare()) {
        throw std::domain_error("Power is only defined for square matrices");
    }
    if (exponent == 0) return identity(rows);
    Matrix base{};
    if (exponent < 0) {
        base = inverse();
        exponent = -exponent;
    } else base = *this;
    Matrix result(identity(rows));
    while (exponent) {
        if (exponent & 1) result *= base;
        exponent >>= 1;
        base *= base;
    }
    return result;
}

//内积
double Matrix::dotProduct(const Matrix& other) const {
    if (!isVector() || !other.isVector() || rows != other.rows) {
        throw std::invalid_argument("Matrices must be column vectors of the same length");
    }
    double result = 0.0;
    for (int i = 0; i < rows; i++) {
        result += data[i][0] * other.data[i][0];
    }
    return result;
}

//迹
double Matrix::trace() const {
    if (!isSquare()) {
        throw std::domain_error("Trace is only defined for square matrices");
    }
    double result = 0.0;
    for (int i = 0; i < rows; i++) {
        result += data[i][i];
    }
    return result;
}

//行阶梯型
Matrix Matrix::echelonForm() const {
    avoidEmpty();
    if (isEchelonForm()) return *this;
    Matrix result(*this);
    int topRow = 0;
    for (int j = 0; j < cols && topRow < rows; j++) {
        int pivotRow = -1;
        for (int i = topRow; i < rows; i++) {
            if (!equalsZero(result(i, j))) {
                pivotRow = i;
                break;
            }
        }
        if (pivotRow == -1) continue;
        if (pivotRow != topRow) result.swapRows(topRow, pivotRow);
        for (int i = topRow + 1; i < rows; i++) {
            if (equalsZero(result(i, j))) continue;
            double scalar = -result(i, j) / result(topRow, j);
            result.addScaledRow(i, topRow, scalar);
        }
        topRow++;
    }
    return result;
}

//简化行阶梯形
Matrix Matrix::reducedEchelonForm() const {
    avoidEmpty();
    if (isReducedEchelonForm()) return *this;
    Matrix result(echelonForm());
    for (int i = rows - 1; i >= 0; i--) {
        int pivotCol = -1;
        for (int j = 0; j < cols; j++) {
            if (!equalsZero(result(i, j))) {
                pivotCol = j;
                break;
            }
        }
        if (pivotCol == -1) continue;
        double pivotValue = result(i, pivotCol);
        if (!equalsZero(pivotValue - 1.0)) {
            result.scaleRow(i, 1.0 / pivotValue);
            result(i, pivotCol) = 1.0;
        }
        for (int r = 0; r < rows; r++) {
            double factor = result(r, pivotCol);
            if (r == i || equalsZero(factor)) continue;
            result.addScaledRow(r, i, -factor);
            if (equalsZero(result(r, pivotCol))) {
                result(r, pivotCol) = 0.0;
            }
        }
    }
    return result;
}

//齐次线性方程组
Matrix Matrix::solve() const {
    avoidEmpty();
    Matrix result(solve(Matrix(rows, 1, 0.0)));
    if (result.cols <= 1) return result;
    return result.getSubMatrix(0, result.rows - 1, 1, result.cols - 1);
}

//非齐次线性方程组
Matrix Matrix::solve(const Matrix& other) const {
    avoidEmpty();
    if (!other.isVector() || rows != other.rows) {
        throw std::invalid_argument("The argument matrix must be a vector with the same rows");
    }
    Matrix augmentedReduced(combine(other).reducedEchelonForm()); //增广的简化阶梯型
    int r = augmentedReduced.rank();
    Matrix reducedEchelon(augmentedReduced.getSubMatrix(0, rows - 1, 0, cols - 1));
    if (reducedEchelon.rank() < r) return Matrix{}; //无解
    Matrix b(augmentedReduced.getColVector(cols));
    if (r == cols) return b; //唯一解
    std::set<int> fixedUnknowns{}; //非自由未知量所在列集合
    int colIndex = -1;
    for (int i = 0; i < rows; i++) {
        bool isZeroRow = true;
        for (int j = colIndex + 1; j < cols; j++) {
            if (!equalsZero(reducedEchelon(i, j))) {
                fixedUnknowns.insert(colIndex = j);
                isZeroRow = false;
                break;
            }
        }
        if (isZeroRow) break;
    }
    Matrix result(b.getSolutionVector(cols, fixedUnknowns)); //特解
    for (int j = 0; j < cols; j++) {
        if (fixedUnknowns.find(j) != fixedUnknowns.end()) continue;
        Matrix colVector(reducedEchelon.getColVector(j));
        Matrix solutionVector(colVector.getSolutionVector(j, cols, fixedUnknowns));
        result = result.combine(solutionVector); //基础解系
    }
    return result;
}

//输入输出
std::ostream& operator<<(std::ostream& os, const Matrix& matrix) {
    for (int i = 0; i < matrix.rows; i++) {
        os << "[ ";
        for (int j = 0; j < matrix.cols; j++) {
            if (equalsZero(matrix.data[i][j])) os << 0.0;
            else os << matrix.data[i][j];
            if (j < matrix.cols - 1) os << ", ";
        }
        os << " ]" << std::endl;
    }
    return os;
}

std::istream& operator>>(std::istream& is, Matrix& matrix) {
    int rows, cols;
    is >> rows >> cols;
    matrix = Matrix(rows, cols);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            is >> matrix(i, j);
        }
    }
    return is;
}

//类型判断
inline bool Matrix::empty() const {
    return rows == 0 || cols == 0;
}

inline bool Matrix::isSquare() const {
    return !empty() && rows == cols;
}

//列向量
inline bool Matrix::isVector() const {
    return cols == 1;
}

//零矩阵
bool Matrix::isZeroMatrix() const {
    if (empty()) return false;
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if (!equalsZero(data[i][j])) return false;
        }
    }
    return true;
}

//奇异
bool Matrix::isSingular() const {
    if (!isSquare()) return false;
    return equalsZero(determinant());
}

//对称
bool Matrix::isSymmetric() const {
    if (!isSquare()) return false;
    for (int i = 0; i < rows; i++) {
        for (int j = i + 1; j < cols; j++) {
            if (!equalsZero(data[i][j] - data[j][i])) {
                return false;
            }
        }
    }
    return true;
}

//对角
bool Matrix::isDiagonal() const {
    if (!isSquare()) return false;
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if (j == i) continue;
            if (!equalsZero(data[i][j])) return false;
        }
    }
    return true;
}

//列满秩
bool Matrix::isFullRank() const {
    if (empty()) return false;
    return rank() == cols;
}

//正交
bool Matrix::isOrthogonal() const {
    if (!isSquare()) return false;
    Matrix multi(*this * transpose());
    Matrix identi(identity(rows));
    return multi == identi;
}

//正定
bool Matrix::isPositiveDefinite() const {
    if (!isSymmetric()) return false;
    Matrix sub(*this);
    for (int n = rows - 1; n >= 0; n--) {
        double d = sub.determinant();
        if (d < 0.0 || equalsZero(d)) return false;
        sub = sub.getSubMatrix(n, n);
    }
    return true;
}

//上三角
bool Matrix::isUpperTriangular() const {
    if (!isSquare()) return false;
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < i; j++) {
            if (!equalsZero(data[i][j])) return false;
        }
    }
    return true;
}

//阶梯型
bool Matrix::isEchelonForm() const {
    if (empty()) return false;
    int index = -1;
    for (int i = 0; i < rows; i++) {
        bool isZeroRow = true;
        for (int j = 0; j < cols; j++) {
            if (!equalsZero(data[i][j])) {
                if (j <= index) return false;
                index = j;
                isZeroRow = false;
                break;
            }
        }
        if (isZeroRow) index = cols;
    }
    return true;
}

//简化阶梯型
bool Matrix::isReducedEchelonForm() const {
    if (!isEchelonForm()) return false;
    int startCol = 0;
    for (int i = 0; i < rows; i++) {
        bool isZeroRow = true;
        for (int j = startCol; j < cols; j++) {
            if (!equalsZero(data[i][j])) {
                if (!equalsZero(data[i][j] - 1.0)) return false;
                for (int k = 0; k < i; k++) {
                    if (!equalsZero(data[k][j])) return false;
                }
                startCol = j + 1;
                isZeroRow = false;
                break;
            }
        }
        if (isZeroRow) break;
    }
    return true;
}

//工具函数
Matrix Matrix::getRowVector(int row) const {
    if (row < 0 || row >= rows) {
        throw std::out_of_range("Matrix index out of range");
    }
    Matrix result(1, cols);
    for (int i = 0; i < cols; i++) {
        result(0, i) = data[row][i];
    }
    return result;
}

Matrix Matrix::getColVector(int col) const {
    if (col < 0 || col >= cols) {
        throw std::out_of_range("Matrix index out of range");
    }
    Matrix result(rows, 1);
    for (int i = 0; i < rows; i++) {
        result(i, 0) = data[i][col];
    }
    return result;
}

void Matrix::print() const {
    std::cout << (*this);
}

inline bool equalsZero(double num) {
    return std::abs(num) < 1.0e-10;
}