C++矩阵计算器

Question

我编写了一个能够用矩阵执行计算的程序：

main.cpp (仅为测试目的)

#include <iostream> 
#include <vector>
#include "matrix.hpp"

int main() { 
    std::cout << "Width (matrix 1):\n";
    int width = getDimension();

    std::cout << "Height (matrix 1):\n";
    int height = getDimension();

    std::vector<std::vector<double>> matrix(height, std::vector<double> (width));

    //Now, the user has to enter the matrix line by line, seperated by commas
    for(int i = 1; i <= height; i++) {
        getUserInput(matrix, i, width);
    }

    //Output
    printMatrix(matrix);
    std::cout << "\n";
    std::cout << "Determinant is " << getDeterminant(matrix) << "\n";
    std::cout << "\nMatrix * Matrix =\n"; 
    printMatrix(getMultiplication(matrix, matrix));
    std::cout << "\nTransposed Matrix:\n"; 
    printMatrix(getTranspose(matrix));
    std::cout << "\nCofactor-Matrix:\n"; 
    printMatrix(getCofactor(matrix));
    std::cout << "\nInverse-Matrix:\n"; 
    printMatrix(getInverse(matrix));
    return 0;
} 

matrix.cpp (负责实际计算)

#include <iostream>
#include <vector>
#include <math.h>
#include <iomanip>
#include <stdexcept>
#include "matrix.hpp"

//Bareiss-Algorithm with O(n^3); not Laplace-expansion with O(n!)
double getDeterminant(std::vector<std::vector<double>> vect) {
    for(int i = 0; i < vect.size(); i++) {
        for(int j = i + 1; j < vect.size(); j++) {
            for(int k = i + 1; k < vect.size(); k++) {
                vect[j][k] = vect[j][k] * vect[i][i] - vect[j][i] * vect[i][k];
                if(i != 0) {
                    vect[j][k] /= vect[i - 1][i - 1];
                }
            }
        }
    }
    return vect[vect.size() - 1][vect.size() - 1];
}

//O(n^3) - much faster isn't possible
std::vector<std::vector<double>> getMultiplication(const std::vector<std::vector<double>> matrix1, const std::vector<std::vector<double>> matrix2) {

    if(matrix1[0].size() != matrix2.size()) {
        throw std::runtime_error("Matrices cannot be multiplicated");
    }
    int height1 = matrix1.size();

    int width2 = matrix2[0].size();
    int height2 = matrix2.size();

    std::vector<std::vector<double>> solution(height1, std::vector<double> (width2));

    //Formula for matrix-multiplication
    for(int i = 0; i < height1; i++) {
        for(int k = 0; k < width2; k++) {
            for(int j = 0; j < height2; j++) {
                solution[i][k] = matrix1[i][j] * matrix2[j][k]; 
            }
        }
    }

    return solution;
}

//O(n^2) - faster is not possible; 
std::vector<std::vector<double>> getTranspose(const std::vector<std::vector<double>> matrix1) {

    //Transpose-matrix: height = width(matrix), width = height(matrix)
    std::vector<std::vector<double>> solution(matrix1[0].size(), std::vector<double> (matrix1.size()));

    //Filling solution-matrix
    for(size_t i = 0; i < matrix1.size(); i++) {
        for(size_t j = 0; j < matrix1[0].size(); j++) {
            solution[j][i] = matrix1[i][j];
        }
    }
    return solution;
}

std::vector<std::vector<double>> getCofactor(const std::vector<std::vector<double>> vect) {
    if(vect.size() != vect[0].size()) {
        throw std::runtime_error("Matrix is not quadratic");
    } 

    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));
    std::vector<std::vector<double>> subVect(vect.size() - 1, std::vector<double> (vect.size() - 1));

    for(std::size_t i = 0; i < vect.size(); i++) {
        for(std::size_t j = 0; j < vect[0].size(); j++) {

            int p = 0;
            for(size_t x = 0; x < vect.size(); x++) {
                if(x == i) {
                    continue;
                }
                int q = 0;

                for(size_t y = 0; y < vect.size(); y++) {
                    if(y == j) {
                        continue;
                    }

                    subVect[p][q] = vect[x][y];
                    q++;
                }
                p++;
            }
            solution[i][j] = pow(-1, i + j) * getDeterminant(subVect);
        }
    }
    return solution;
}

std::vector<std::vector<double>> getInverse(const std::vector<std::vector<double>> vect) {
    double det = getDeterminant(vect);
    if(det == 0) {
        throw std::runtime_error("Determinant is 0");
    } 
    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));

    solution = getTranspose(getCofactor(vect));

    for(size_t i = 0; i < vect.size(); i++) {
        for(size_t j = 0; j < vect.size(); j++) {
            solution[i][j] *= (1/det);
        }
    }

    return solution;
}

int getDimension() {
    int dimension;
    std::cout << "Please enter dimension of Matrix: ";
    std::cin >> dimension;
    std::cout << "\n";

    if(dimension < 0 || std::cin.fail()) {
        std::cin.clear();
        std::cin.ignore();
        std::cout << "ERROR: Dimension cannot be < 0.\n";
        return getDimension();
    }

    return dimension;
}

void printMatrix(const std::vector<std::vector<double>> vect) {
    for(std::size_t i = 0; i < vect.size(); i++) {
        for(std::size_t j = 0; j < vect[0].size(); j++) {
            std::cout << std::setw(8) << vect[i][j] << " ";
        }
        std::cout << "\n";
    }
}

void getUserInput(std::vector<std::vector<double>>& vect, int i, int dimension) {
    std::string str = "";
    std::cout << "Enter line " << i << " only seperated by commas: ";
    std::cin >> str;
    std::cout << "\n";
    str = str + ',';

    std::string number = "";
    int count = 0;

    for(std::size_t k = 0; k < str.length(); k++) {
        if(str[k] != ',') {
            number = number + str[k];
        }
        else if(count < dimension) {
            if(number.find_first_not_of("0123456789.-") != std::string::npos) {
                std::cout << "ERROR: Not only numbers entered.\n";
                getUserInput(vect, i, dimension);
                break;
            }
            else if(number.find_first_not_of("0123456789-.") == std::string::npos) {
                vect[i - 1][count] = std::stod(number);
                number = "";
                count++;
            }
            else {
                std::cout << "ERROR: Not enough numbers entered.\n";
                getUserInput(vect, i, dimension);
                break;
            }       
        }
        else {
            std::cout << "ERROR: Too many numbers entered.\n";
            getUserInput(vect, i, dimension);
            break;
        }
    }    
}

matrix.hpp

#ifndef MATRIX_HPP
#define MATRIX_HPP

#include <vector>

double getDeterminant(const std::vector<std::vector<double>> vect);
std::vector<std::vector<double>> getMultiplication(const std::vector<std::vector<double>> matrix1, const std::vector<std::vector<double>> matrix2);
std::vector<std::vector<double>> getTranspose(const std::vector<std::vector<double>> matrix1);
std::vector<std::vector<double>> getCofactor(const std::vector<std::vector<double>> vect);
std::vector<std::vector<double>> getInverse(const std::vector<std::vector<double>> vect);
int getDimension();
void printMatrix(const std::vector<std::vector<double>> vect);
void getUserInput(std::vector<std::vector<double>>& vect, int i, int dimension);
#endif

我的问题：

• getCofactor()和getInverse()的运行时(O-表示法)是什么？
• 是否有更有效的方法来完成这些任务？
• 如何对代码进行总体改进？

Answer 1

创建一个class Matrix

创建一个class Matrix并向其添加成员函数来操作矩阵，而不是有一个向量向量，并且具有操纵这些向量的全局函数。您可能还应该为算术操作创建重载，这样您就可以编写类似于auto matrix3 = matrix1 + matrix2;的东西了。

看看现有C++矩阵库，看看什么是可能的。

复杂性的getCofactor()与getInverse()

在getCofactor()中，您有四个嵌套的for-loops，它们仅在x == i或y == j情况下跳过迭代，因此将使其成为O(N^4)。但是，在大小为(N1)^2的子向量上的第二个最外层循环中调用getDeterminant()，这样就会使它成为O(N^2 * (N-1)^3) = O(N^5)。

在getInverse()中，复杂度主要取决于对getCofactor()的调用，因此它也是O(N^5)。这对大型矩阵确实不好，因为你应该在第一年的数学学习中学习的Gauss方法就是O(N^3)。

维基百科有一个矩阵代数算法列表复杂性，您可以检查看看什么是可能的。但是，请注意，具有奇怪指数的最佳算法可能很难实现，实际上，对于小矩阵而言，效率可能要比简单的算法低得多。因此，对您来说最好的算法取决于程序将要处理的矩阵的大小。