From 573947fe2df601d84e983b9a38b362864edc93d7 Mon Sep 17 00:00:00 2001
From: Francois Coppens <francois.coppens@uvsq.fr>
Date: Fri, 11 Jun 2021 08:46:39 +0200
Subject: [PATCH] Woodbury debugging...

---
 include/Helpers.hpp |  2 +-
 src/Woodbury.cpp    | 14 ++++++++++++++
 tests/test_h5.cpp   |  2 ++
 3 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/include/Helpers.hpp b/include/Helpers.hpp
index 2384db1..cd6fdbb 100644
--- a/include/Helpers.hpp
+++ b/include/Helpers.hpp
@@ -117,7 +117,7 @@ template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
 }
 
 template <typename T1, typename T2, typename T3>
-static inline void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) {
+void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) {
   for(unsigned int i = 0; i < M; i++) {
     for(unsigned int j = 0; j < P; j++) {
       C[i * P + j] = 0;
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 1c92a68..3b8e26a 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -19,12 +19,26 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
 
   // Construct V from Updates_index
   unsigned int V[2 * Dim]; // 2 x Dim matrix stored in row-major order
+  std::memset(V, 0, 2 * Dim * sizeof(unsigned int));
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 1] = 1;
 
+  showMatrix2(Slater_inv, Dim, Dim, "Slater_inv");
+  showMatrix2(Updates, 2, Dim, "Updates");
+  showMatrix2(Updates_index, 1, 2, "Updates_index");
+  showMatrix2(V, 2, Dim, "V");
+
   // Compute C
   double C[2 * Dim];
   matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
+
+  int A[6] = {1,2,3,4,5,6};
+  int Y[12] = {7,8,9,10,11,12,13,14,15,16,17,18};
+  int Z[8] = {0,0,0,0,0,0,0,0};
+  matMul2(A,Y,Z,2,3,4);
+  showMatrix2(Z, 2,4,"Z");
+
+  showMatrix2(C, 2, Dim, "C = S_inv * Updates");
   // Compute B
   double B[4];
   matMul2(V, C, B, 2, Dim, 2);
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index f750f10..9b49a1c 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -52,6 +52,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   showMatrix(slater_inverse, dim, "OLD Inverse");
 #endif
 
+  showMatrix2(slater_matrix, dim, dim, "Slater");
+  
   // Transform replacement updates in 'updates[]' into additive updates in 'u[]'
   for (j = 0; j < nupdates; j++) {
     for (i = 0; i < dim; i++) {