From 7d3ffe1d2a6cf95fd9fee242f3f28a8efb217c17 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <francois.coppens@uvsq.fr>
Date: Fri, 29 Jan 2021 12:53:20 +0100
Subject: [PATCH] Added random matrix initialisation of arbitrary size. Added
 function to calculate determinant to test matrix invertability of A and A0.

---
 SM-MaponiA3.cpp | 66 +++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 56 insertions(+), 10 deletions(-)

diff --git a/SM-MaponiA3.cpp b/SM-MaponiA3.cpp
index 4d2c74f..d6c6346 100644
--- a/SM-MaponiA3.cpp
+++ b/SM-MaponiA3.cpp
@@ -3,6 +3,9 @@
 
 #include <iostream>
 #include <string>
+#include <cmath>
+#include <cstdlib>
+#include <ctime>
 using namespace std;
 
 uint getMaxIndex(double* arr, uint size);
@@ -12,9 +15,12 @@ template<typename T>void showMatrix(T** matrix, uint size, string name);
 template<typename T>void showMatrixT(T** matrix, uint size, string name);
 template<typename T>T** matMul(T** A, T** B, uint size);
 template<typename T1, typename T2>T1** outProd(T1* vec1, T2* vec2, uint size);
+template<typename T>T matDet(T** A, int M);
 
 int main() {
 
+    srand((unsigned) time(0));
+    uint randRange = 1; // to get random integers in range [-randRange, randRange]
     uint M = 3;
     uint i, j, k, l, lbar, tmp;
     double alpha, beta;
@@ -49,8 +55,6 @@ int main() {
     // Initialize all matrices with zeros
     for (i = 0; i < M; i++) {
         for (j = 0; j < M; j++) {
-            A[i][j] = 0;
-            //Ainv[i][j] = 0;
             A0[i][j] = 0;
             A0inv[i][j] = 0;
             Ar[i][j] = 0;
@@ -67,10 +71,17 @@ int main() {
         }
     }
 
-    // Initialize A
-    A[0][0] = 1;    A[0][1] = 1;    A[0][2] = -1;
-    A[1][0] = 1;    A[1][1] = 1;    A[1][2] = 0;
-    A[2][0] = -1;   A[2][1] = 0;    A[2][2] = -1;
+    // // Initialize A with M=3 and Eq. (17) from paper
+    // A[0][0] = 1;    A[0][1] = 1;    A[0][2] = -1;
+    // A[1][0] = 1;    A[1][1] = 1;    A[1][2] = 0;
+    // A[2][0] = -1;   A[2][1] = 0;    A[2][2] = -1;
+
+    // Fill A with random numbers from [-1,1]
+    for (i = 0; i < M; i++) {
+        for (j = 0; j < M; j++) {
+            A[i][j] = rand()%(2*randRange+1)-randRange;
+        }
+    }
 
     // Define identity matrix, A0, A0inv and p
     p[0] = 0;
@@ -100,7 +111,7 @@ int main() {
         }
         // showVector(ylk[0][k], M+1, "y0k");
     }
-    showMatrixT(ylk[0], M+1, "y0k");
+    // showMatrixT(ylk[0], M+1, "y0k");
 
     // Calculate all the ylk from the y0k
     // showVector(p, M+1, "p");
@@ -131,9 +142,8 @@ int main() {
             // showVector(ylk[l][p[k]], M+1, "ylk");
         }
     }
-    showMatrixT(ylk[1], M+1, "y1k");
-    showMatrixT(ylk[2], M+1, "y2k");
-    // EVERYTHING WORKS UPTO HERE
+    // showMatrixT(ylk[1], M+1, "y1k");
+    // showMatrixT(ylk[2], M+1, "y2k");
 
     // Construct A-inverse from A0-inverse and the ylk
     double** U;
@@ -257,4 +267,40 @@ T1** outProd(T1* vec1, T2* vec2, uint size) {
         }
     }
     return C;
+}
+
+template<typename T>
+T matDet(T** A, int M) {
+    int det = 0, p, h, k, i, j;
+    T** temp = new T*[M];
+    for (int i = 0; i < M; i++) temp[i] = new T[M];
+    if(M==1) {
+        return A[0][0];
+    }
+    else if(M == 2) {
+        det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]);
+        return det;
+    }
+    else {
+        for(p = 0; p < M; p++) {
+            h = 0;
+            k = 0;
+            for(i = 1; i < M; i++) {
+                for( j = 0; j < M; j++) {
+                    if(j == p) {
+                        continue;
+                    }
+                    temp[h][k] = A[i][j];
+                    k++;
+                    if(k == M-1) {
+                        h++;
+                        k = 0;
+                    }
+                }
+            }
+            det = det + A[0][p] * pow(-1, p) * matDet(temp, M-1);
+        }
+        return det;
+    }
+    delete [] temp;
 }
\ No newline at end of file