Created Woodbery kernel infrastructure and started wrinting Woodbery 2x2 kernel.

Makefile

@@ -42,7 +42,9 @@ FFLAGS = $(CXXFLAGS)
 INCLUDE = -I $(INC_DIR)/
 DEPS_CXX = $(OBJ_DIR)/SM_Maponi.o \
 		   $(OBJ_DIR)/SM_Standard.o \
-		   $(OBJ_DIR)/SM_Helpers.o
+		   $(OBJ_DIR)/Woodbury.o \
+		   $(OBJ_DIR)/SMWB.o \
+		   $(OBJ_DIR)/Helpers.o
 DEPS_F = $(DEPS_CXX) \
 		 $(OBJ_DIR)/finterface_mod.o \
 		 $(OBJ_DIR)/helpers_mod.o

include/Helpers.hpp

@@ -1,4 +1,4 @@
-// SM_Helpers.hpp
+// Helpers.hpp
 // Some usefull helper functions to support the Maponi algorithm.
 #include <cmath>
 #include <cstring>

include/SMWB.hpp

@@ -0,0 +1,4 @@
+// Sherman-Morrison-Woodbury kernel 1
+// WB2, WB3, SM2 mixing scheme 1
+void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index);

include/Woodbury.hpp

@@ -0,0 +1,7 @@
+// Woodbury 2x2 kernel
+void WB2(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index);
+
+// Woodbury 3x3 kernel
+void WB3(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index);

src/Helpers.cpp

@@ -1,4 +1,4 @@
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 
 // Set common break-down threshold
 double threshold() {

src/SMWB.cpp

@@ -0,0 +1,21 @@
+#include "SMWB.hpp"
+#include "SM_Standard.hpp"
+#include "Woodbury.hpp"
+#include "Helpers.hpp"
+
+// Sherman-Morrison-Woodbury kernel 1
+// WB2, WB3, SM2 mixing scheme 1
+void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index) {
+           std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
+           WB3(Slater_inv, Dim, Updates, Updates_index);
+           WB2(Slater_inv, Dim, Updates, Updates_index);
+           SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
+         }
+
+extern "C" {
+void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
+           double **linUpdates, unsigned int **Updates_index) {
+  SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+}
+}

src/SM_Maponi.cpp

@@ -2,7 +2,7 @@
 // Algorithm 3 from P. Maponi,
 // p. 283, doi:10.1016/j.laa.2006.07.007
 #include "SM_Maponi.hpp"
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 
 // #define DEBUG
 

src/SM_Standard.cpp

@@ -1,7 +1,7 @@
 // SM-Standard.cpp
 // Standard Sherman Morrison with multiple updates
 #include "SM_Standard.hpp"
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 
 // Naïve Sherman Morrison
 void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,

src/Woodbury.cpp

@@ -0,0 +1,69 @@
+// Woodbury.cpp
+// Woodbury 2x2 and 3x3 kernels
+//
+// Woodbury matrix identity:
+// (S + U * V)^{-1} = S^{-1} - S^{-1} * U * B^{-1} * V * S^{-1}
+// B := 1 + V * C,  2 x 2
+// C := S^{-1} * U, dim x 2
+#include "Woodbury.hpp"
+#include "Helpers.hpp"
+
+// Woodbury 2x2 kernel: 
+void WB2(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index) {
+  std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
+
+  // Construct V from Updates_index
+  unsigned int *V = new unsigned int[2 * Dim]{0};
+  for (unsigned int i = 0; i < Dim; i++) {
+    if (i == Updates_index[0] - 1) V[i] = 1;
+    if (Dim + i == Updates_index[1] - 1) V[Dim + i] = 1;
+  }
+
+  // Compute B from U, V and Slater_inv
+  double B[4]{0};
+  double Binv[4]{0};
+  double C[4]{0};
+  double D[2*Dim]{0};
+    // Compute C
+    for (unsigned int i = 0; i < Dim; i++) {
+      for (unsigned int j = 0; j < 2; j++) {
+        for (unsigned int k = 0; k < Dim; k++) {
+          C[i * Dim + j] += Slater_inv[i * Dim + k] * Updates[k * Dim + j];
+        }
+      }
+    }
+    // Compute B
+    for (unsigned int i = 0; i < 2; i++) {
+      for (unsigned int j = 0; j < 2; j++) {
+        for (unsigned int k = 0; k < Dim; k++) {
+          B[i * Dim + j] += (i == j) + V[i * Dim + k] * C[k * Dim + j];
+        }
+      }
+    }
+
+  // Invert B with explicit formula for 2x2 inversion
+  double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]);
+  Binv[0] = idet * B[3];
+  Binv[1] = -1.0 * idet * B[1];
+  Binv[2] = -1.0 * idet * B[2];
+  Binv[3] = idet * B[0];
+
+  // Compute (S + U * V)^{-1} with Woobury identity
+  
+}
+
+// Woodbury 3x3 kernel
+void WB3(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index) {
+           std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
+
+           // Construct V from Updates_index
+
+           // Compute B from U, V and Slater_inv
+
+           // Invert B with explicit formula for 3x3 inversion
+
+           // Compute (S + U * V)^{-1} with Woobury identity
+           
+         }

src/finterface_mod.f90

@@ -43,5 +43,12 @@ module Sherman_Morrison
             real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
             real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
         end subroutine SM4
-    end interface
+        subroutine SMWB1(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SMWB1_f")
+            use, intrinsic :: iso_c_binding, only : c_int, c_double
+            integer(c_int), intent(in) :: dim, n_updates
+            integer(c_int), dimension(:), allocatable, intent(in) :: Updates_index
+            real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
+            real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
+        end subroutine SMWB1
+        end interface
 end module Sherman_Morrison

tests/cMaponiA3_test_3x3_3.cpp

@@ -1,5 +1,5 @@
 // main.cpp
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 #include "SM_Maponi.hpp"
 
 int main() {

tests/test_h5.cpp

@@ -1,9 +1,10 @@
 #include "hdf5/serial/H5Cpp.h"
 #include "hdf5/serial/hdf5.h"
 
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 #include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
+#include "SMWB.hpp"
 
 using namespace H5;
 // #define DEBUG
@@ -81,6 +82,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
     SM3(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "sm4") {
     SM4(slater_inverse, dim, nupdates, u, col_update_index);
+  } else if (version == "smwb1") {
+    SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
 #ifdef MKL
   } else if (version == "lapack") {
     memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double));

tests/vfc_test_h5.cpp

@@ -10,7 +10,7 @@
 #include <sstream>
 #include <vector>
 
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 #include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
 #include "vfc_probe.h"

tools/restyle.sh

@@ -9,7 +9,7 @@ then
 	exit 1
 fi
 
-for ext in c cpp h hpp
+for ext in c cc cpp h hpp
 do
     find $SMROOT -type f -iname "*.${ext}" -exec echo "$FORMATER $STYLE" {} \;
-done
+done