0

I'm implementing PLU decomposition algorithm based on CLRS on c++, however, the permutation matrix is really strange:

class Matrix{
    public:
        Matrix(int);
        ~Matrix();
        Matrix PLU_decomposition(Matrix&, Matrix&);
        void exchange_rows(int, int);
        void print_data();
        int size;
        double** data;
};

Matrix::Matrix(int n){
    data = new double* [n];
    for(int i=0; i<n; i++){
        data[i] = new double[n];
    }
    for(int i=0;i<n; i++){
        for(int j=0; j<n; j++){
            data[i][j] = 0;
        }
    }
    size = n;
}

Matrix::~Matrix(){
    for(int i=0; i<size; i++){
        delete[] data[i];
    }
    delete[] data;
}

void Matrix::exchange_rows(int a, int b){
    for(int i=0; i<size; i++){
        double temp = data[a][i];
        data[a][i] = data[b][i];
        data[b][i] = temp;
    }
    return;
}

void Matrix::print_data(){
    for(int i=0; i<size; i++){
        for(int j=0; j<size; j++){
            cout << data[i][j] << " ";
        }
        cout << "\n";
    }
}
Matrix Matrix::PLU_decomposition(Matrix& L, Matrix& U){ //PA=LU
    //L U will be modified after this, and outputs a permutation P
    int n=size; 
    Matrix P(n);
    for(int i=0; i<n; i++){
        P.data[i][i] = 1;
    }
    for(int k=0; k<n; k++){
        double p=0;
        int k_prime = k;
        for(int i=k; i<n; i++){
            if(abs(data[i][k])>p){
                p = abs(data[i][k]);
                k_prime = i;
            }
        }
        if(p==0) cout << "Error in PLU_Decomposition: singular matrix\n";
        else{
            P.exchange_rows(k, k_prime);
            exchange_rows(k, k_prime);  
        }
        for(int i=k+1; i<n; i++){
            data[i][k] = data[i][k] /data[k][k];
            for(int j=k+1; j<n; j++){
                data[i][j] = data[i][j] - data[i][k]*data[k][j];
            } 
        }
    }
    for(int i=0; i<n;i++){
        for(int j=0; j<n; j++){
            if(i>j){
                L.data[i][j] = data[i][j] ;
            }else{
                if(i==j){
                    L.data[i][j] = 1;
                }
                U.data[i][j] = data[i][j] ;
            }
        }
    }
    return P;
}

int main(){
    Matrix m(3), P(3), L(3), U(3);
    m.data[0][0]=1;
    m.data[0][1]=2;
    m.data[0][2]=0;
    m.data[1][0]=3;
    m.data[1][1]=4;
    m.data[1][2]=4;
    m.data[2][0]=5;
    m.data[2][1]=6;
    m.data[2][2]=3;

    P = m.PLU_decomposition(L, U);
    P.print_data();
}

The L U are all correct, but the problem is P. I initialize it with some 1s only, however, at the end, the output is

4.68746e-310 4.68746e-310 0
1.97626e-323 0 0
4.68746e-310 4.68746e-310 0
free(): double free detected in tcache 2
Aborted

Can anybody tell me why there's something wrong with deleting the object in the destructor?

jerry
  • 37
  • 5
  • It may not be easy to understand how the duplicate questions helps at first, but read them through. In short: You need to implement a copy constructor, a move constructor, a copy assignment operator and a move assignment operator. Also read [The rule of three/five/zero](https://en.cppreference.com/w/cpp/language/rule_of_three) – Ted Lyngmo Jun 12 '22 at 07:06
  • You could also avoid having to implement all those constructors and operators by using a `std::vector>` for `data` instead of manually allocating the memory. [Example](https://godbolt.org/z/rh1YfTs5v). You should almost never actually use `new` / `delete` and `new[]`/`delete[]` in modern C++. – Ted Lyngmo Jun 12 '22 at 07:18

0 Answers0