Untitled

#include <iostream>
#include <numa.h>
#include <omp.h>
#include <iostream>
#include <vector>
#include <algorithm> 
#include <cmath> 
#include <cstdlib>
#include <sys/time.h>
#include <iomanip>
using namespace std;

typedef vector<std::vector<double>> Matrix;

void PLU_Decomposition(Matrix& A, Matrix& P, Matrix& L, Matrix& U) {
    int n = A.size();

    // Initialize L and U
    L = Matrix(n, vector<double>(n, 0.0));
    U = Matrix(n, vector<double>(n, 0.0));
    P = Matrix(n, vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) {
        P[i][i] = 1.0; // Initialize P as an identity matrix
        L[i][i] = 1.0; // Diagonal of L is 1
    }

    for (int k = 0; k < n; ++k) {
        double max = 0.0;
        int maxIndex = k;
        // Find the pivot row.
        #pragma omp parallel for reduction(max:max)
        for (int i = k; i < n; ++i) {
            if (fabs(A[i][k]) > max) {
                max = fabs(A[i][k]);
                #pragma omp critical
                {
                    maxIndex = i;
                }
            }
        }


        if (maxIndex == 0.0) {
            cerr << "Matrix is singular!" << endl;
            return;
        }

        // Swap rows in P, A_copy
        swap(P[k], P[maxIndex]);
        swap(A[k], A[maxIndex]);

        // Perform the elimination process.
        #pragma omp parallel for
        for (int i = k + 1; i < n; ++i) {
            double factor = A[i][k] / A[k][k];
            for (int j = k; j < n; ++j) {
                A[i][j] -= factor * A[k][j];
            }
            L[i][k] = factor; // Update the L matrix.
        }
    }

    // Copying the upper triangle of A to U.
    for (int i = 0; i < n; ++i)
        for (int j = i; j < n; ++j)
            U[i][j] = A[i][j];
}

// Function to compute the matrix product
Matrix Multiply(const Matrix& A, const Matrix& B) {
    int n = A.size();
    Matrix result(n, vector<double>(n, 0.0));

    // #pragma omp parallel for collapse(2)
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            double sum = 0.0;
            for (int k = 0; k < n; ++k) {
                result[i][j] += A[i][k] * B[k][j];
            }
            result[i][j] = sum;   // No data race because each thread have unique[i][j]
        }
    }
    return result;
}


int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: " << argv[0] << " <matrix size>" << std::endl;
        return -1;
    }

    // Parse the matrix size from the command line arguments
    int n = atoi(argv[1]);
    int nthreads = atoi(argv[2]);

    omp_set_num_threads(nthreads);

    // Seed the random number generator
    srand48(time(nullptr));

    Matrix A(n, vector<double>(n));
    // Initialize A with random values
    // #pragma omp parallel for collapse(2)
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            A[i][j] = drand48() * 100; // Random double between 0 and 100
        }
    }

    Matrix A_original, P, L, U;
    A_original = A;     // Preserve the original matrix A
    
   struct timeval start, end;

   gettimeofday(&start,NULL);

   PLU_Decomposition(A, P, L, U);

   gettimeofday(&end, NULL);

   double diff = (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1000000.0;
   cout << "Time: " << fixed << setprecision(2) << diff <<  "s" << endl;
    // cout << "P:" << std::endl;
    // for (auto& row : P) {
    //     for (auto& element : row) {
    //         cout << element << " ";
    //     }
    //     cout << endl;
    // }

    // cout << "L:" << std::endl;
    // for (auto& row : L) {
    //     for (auto& element : row) {
    //         cout << element << " ";
    //     }
    //     cout << endl;
    // }

    // cout << "U:" << endl;
    // for (auto& row : U) {
    //     for (auto& element : row) {
    //         cout << element << " ";
    //     }
    //     cout << endl;
    // }

    // // Step 1: Compute PA
    // Matrix PA = Multiply(P, A_original);

    // // Step 2: Compute LU
    // Matrix LU = Multiply(L, U);

    // // Step 3: Calculate the residual matrix R = PA - LU
    // Matrix R(n, vector<double>(n, 0.0));
    // // #pragma omp parallel for collapse(2)
    // for (int i = 0; i < n; ++i) {
    //     for (int j = 0; j < n; ++j) {
    //         R[i][j] = PA[i][j] - LU[i][j];
    //     }
    // }

    // // Step 4: Compute the Euclidean norm of each column of R
    // // Step 5: Sum the norms to get the L2,1 norm
    // double L21_norm = 0.0;
    // // #pragma omp parallel for reduction(+:L21_norm)
    // for (int j = 0; j < n; ++j) {
    //     double col_norm = 0.0;
    //     for (int i = 0; i < n; ++i) {
    //         col_norm += R[i][j] * R[i][j];
    //     }
    //     L21_norm += sqrt(col_norm);
    // }

    // // Print the L2,1 norm of the residual
    // cout << "L2,1 norm of the residual: " << L21_norm << endl;


    return 0;
}