Fix Issue #9: Add partial pivoting to LU Factorization

src/LU_factorisation.cpp

@@ -1,17 +1,24 @@
 // lu_decomposition.cpp
-#include "lu_decomposition.h"
+#include "../include/lu_decomposition.h"
 #include <omp.h>
 #include <iostream>
 #include <cstdio>
 #include <cstdlib>
+#include <cmath>
+
 using namespace std;
 // LU Decomposition function
-void l_u_d(float** a, float** l, float** u, int size)
+void l_u_d(float** a, float** l, float** u, int* p , int size)  // Added pivot array p
 {
     // Initialize a simple lock for parallel region
     omp_lock_t lock;
     omp_init_lock(&lock);
 
+    // Initialize permutation array p
+    for(int i = 0 ; i < size ; i++){
+        p[i] = i ;
+    }
+
     // Initialize L and U matrices
     for (int i = 0; i < size; i++) {
         for (int j = 0; j < size; j++) {
@@ -29,9 +36,30 @@ void l_u_d(float** a, float** l, float** u, int size)
     }
 
     // Parallelize the LU decomposition
-    #pragma omp parallel shared(a, l, u)
+    #pragma omp parallel shared(a, l, u , p)  // Add p to shared
     {
         for (int k = 0; k < size; k++) {
+            // Find pivot row
+            int pivot = k ;
+            float max_val = fabs(a[k][k]) ;
+            for(int i = k+1 ; i < size ; i++){
+                if(fabs(a[i][k] > max_val)){
+                    pivot = i ;
+                    max_val = fabs(a[i][k]) ;
+                }
+            }
+
+            // Row swapping
+            if(pivot != k){
+                omp_set_lock(&lock) ;
+                swap(p[k] , pivot) ; // Swapping in permutation array
+                for(int j = 0 ; j < size ; j++){
+                    if(j < k) swap(a[k][j] , l[pivot][j]) ;
+                }
+                omp_unset_lock(&lock);
+            }
+
+
             // Update U matrix
             #pragma omp for schedule(static)
             for (int j = k; j < size; j++) {
@@ -62,11 +90,13 @@ void l_u_d(float** a, float** l, float** u, int size)
 int main(int argc, char *argv[]) {
     int size = 2;
     float **a, **l, **u;
+    int *p ;  // Added permutation array
 
     // Allocate memory for the 2D arrays
     a = (float **)malloc(size * sizeof(float *));
     l = (float **)malloc(size * sizeof(float *));
     u = (float **)malloc(size * sizeof(float *));
+    p = (int *)malloc(size * sizeof(int));  // Allocated memory for pivot array
     for (int i = 0; i < size; i++) {
         a[i] = (float *)malloc(size * sizeof(float));
         l[i] = (float *)malloc(size * sizeof(float));
@@ -85,7 +115,7 @@ int main(int argc, char *argv[]) {
     }
 
     // Perform LU decomposition
-    l_u_d(a, l, u, size);
+    l_u_d(a , l , u , p , size) ; // Passed pivot array
 
     // Print L matrix
     printf("L Matrix:\n");