使用高斯消元的模板中的矩阵行列式

Question

我在矩阵模板中实现矩阵的行列式时遇到问题。我试图通过计算减少的关联矩阵中主对角线的乘积来获得行列式。问题是它只适用于某些情况并且不可靠。这是行列式的代码

template <typename T>
T Matrix<T>::Det() const  {
    if (Rows != Cols) {
        cout << "Matrix must be square" << endl;
    }
    Matrix<T> r = Reduced();
    T Det = 1;
    for (int i=0; i<Rows; i++) {
        Det *= r.getValue(i, i);
    }
    return Det;
}

由于高斯消元有效，我为任何给定矩阵（符合要减少的标准）获得的减少矩阵很好，所以我认为不会有任何问题，但我没有能够看到我错过了什么。欢迎任何建议！

编辑： 正如有人指出的那样，这是一个最小的可重现示例。这将是模板

#ifndef MATRIX_H
#define MATRIX_H

#include <vector>
#include <iostream>
using namespace std;

template <typename T>
class Matrix {
private:
    unsigned int Rows;
    unsigned int Cols;
    T *Mat;
public:
    Matrix(unsigned int Dim);
    Matrix(unsigned int Rows, unsigned int Cols);
    Matrix(unsigned int Rows, unsigned int Cols, const std::vector<T>& Vec);
    Matrix(unsigned int Dim, const std::vector<T>& Vec);
    Matrix(const Matrix<T>& M);
    ~Matrix();
    T& getValue(unsigned int Row, unsigned int Col) const;
    T Det() const;
    Matrix<T> Reduced() const;
};

template <typename T>
Matrix<T>::Matrix(unsigned int Rows, unsigned int Cols) 
: Rows(Rows), Cols(Cols) {
    if (Rows<=0 || Cols<=0) {
        cout << "Número inválido de filas o columnas" << endl;
    }
    Mat = new T[Rows*Cols];
    for (int i=0; i<Rows*Cols; i++) {
        Mat[i] = 0;
    }
}

template <typename T>
Matrix<T>::Matrix(unsigned int Dim) 
: Rows(Dim), Cols(Dim) {
    if(Rows<=0 || Cols<=0) {
        cout << "Número inválido de filas o columnas" << endl;
    }
    Mat = new T[Rows*Cols];
    for (int i=0; i<Rows*Cols; i++) {
        Mat[i] = 0;
    }
}

template <typename T>
Matrix<T>::Matrix(unsigned int Rows, unsigned int Cols, const std::vector<T>& Vec) 
: Rows(Rows), Cols(Cols) {
    if (Rows<=0 || Cols<=0) {
        cout << "Número inválido de filas o columnas" << endl;
    }
    Mat = new T[Rows*Cols];
    if (Vec.size() != Rows*Cols) {
        cout << "Los tamaños de la matriz y el vector no son iguales" << endl;
    }
    for (int i=0; i<Rows*Cols; i++) {
        Mat[i] = Vec[i];
    }
}

template <typename T>
Matrix<T>::Matrix(unsigned int Dim, const std::vector<T>& Vec) 
: Rows(Dim), Cols(Dim) {
    if (Rows<=0 || Cols<=0) {
        cout << "Número inválido de filas o columnas" << endl;
    }
    Mat = new T[Rows*Cols];
    if (Vec.size() != Rows*Cols) {
        cout << "Número inválido de filas o columnas" << endl;
    }
    for (int i=0; i<Rows*Cols; i++) {
        Mat[i] = Vec[i];
    }
}

template <typename T>
Matrix<T>::Matrix(const Matrix<T>& M)
: Rows(M.Rows), Cols(M.Cols), Mat(new T[Rows * Cols]) {
    for (int i = 0; i < Rows * Cols; i++)
        Mat[i] = M.Mat[i];
}

template <typename T>
Matrix<T>::~Matrix() {
    delete[] Mat;
}

template <typename T>
T &Matrix<T>::getValue(unsigned int Row, unsigned int Col) const  {
    if (Row<0 || Row>=Rows || Col<0 || Col>=Cols) {
        cout << "Índice incorrecto" << endl;
    }
    return Mat[Row*Cols + Col];
}

template <typename T>
Matrix<T> Matrix<T>::Reduced() const {
    if (Rows != Cols) {
        cout << "La matriz debe ser cuadrada para reducirse" << endl;
    }
    Matrix<T> Tri(*this);
    
    int n = Rows;
    int m = 0;
    for (int k=0; k<n-1; k++) {
        if (Tri.Mat[k*Cols + k] == 0)
            cout << "La matriz es singular" << endl;
        
        for (int i = k+1; i<n; i++) {
            m = Tri.Mat[i*Cols + k]/Tri.Mat[k*Cols + k];
            for (int j = k+1; j<n; j++) {
                Tri.Mat[i*Cols + j] = Tri.Mat[i*Cols + j] - m*Tri.Mat[k*Cols + j];
            }
            Tri.Mat[i*Cols + k] = 0;
        }
    }
    return Tri;
}

template <typename T>
T Matrix<T>::Det() const  {
    if (Rows != Cols) {
        cout << "Matrix must be square" << endl;
    }
    Matrix<T> r = Reduced();
    T Det = 1;
    for (int i=0; i<Rows; i++) {
        Det *= r.getValue(i, i);
    }
    return Det;
}

template <class T>
ostream & operator<<(ostream &os, const Matrix<T> &Shw) {
    for (int i=0; i<Shw.getRows(); i++) {
        os << "| ";
        for (int j=0; j<Shw.getCols(); j++) {
            os << Shw.getValue(i,j) << " ";
        }
        os << "|\n";
    }
    return os;
}
#endif

这是 main.cc

上的实现

#include "matrix.h"
int main() {
    Matrix<double> M1(2, 2, {1, 3, 5, 7});
    cout << M1.Det() << endl;
    
    Matrix<double> M2(3, 3, {1, 3, 5, 7, 8, 9, 11, 2, 14});
    cout << M2.Det() << endl;
    
    return 0;
}

输出如下

-8    //This one is okay
-143  //This one is not

Answer 1

经过一些检查，我想我发现了你的问题。您在定义关键变量之一 m 时犯了一个简单的错误。如您所知，这是高斯消元法中的主要定标器，它的工作是将下一行的元素归零。与其将其定义为 int m = 0，不如将其定义为 float m = 0 或 double m = 0.
此外，在您计算 m 的行中，您需要将除法的操作数强制转换为适当的类型。所以在重新定义 m 之后，你应该改变这一行：

m = Tri.Mat[i*Cols + k]/Tri.Mat[k*Cols + k]

至此

m = Tri.Mat[i*Cols + k] / (double)Tri.Mat[k*Cols + k]

请记住，您的模板类型也应为 double，否则将无法正常工作。

此外，我建议您查看 here 以了解高斯消元的一些限制。在某些情况下算法无法有效工作。

使用高斯消元的模板中的矩阵行列式

Matrix determinant in template using gaussian elimination

c++

templates

matrix

linear-algebra