First implementation of Woodbury 2x2 and 3x3 kernels.

include/Helpers.hpp

@@ -56,6 +56,45 @@ void showMatrix(T *matrix, unsigned int M, std::string name) {
   std::cout << std::endl;
 }
 
+template <typename T>
+void showMatrix2(T *matrix, unsigned int M, unsigned int N, std::string name) {
+  std::cout.precision(17);
+  std::cout << name << " = [" << std::endl;
+  for (unsigned int i = 0; i < M; i++) {
+    std::cout << "[";
+    for (unsigned int j = 0; j < N; j++) {
+      if (matrix[i * N + j] >= 0) {
+        std::cout << "  " << matrix[i * N + j] << ",";
+      } else {
+        std::cout << " " << matrix[i * N + j] << ",";
+      }
+    }
+    std::cout << " ]," << std::endl;
+  }
+  std::cout << "]" << std::endl;
+  std::cout << std::endl;
+}
+
+
+template <typename T>
+void showMatrixNS(T *matrix, unsigned int M, unsigned int N, std::string name) {
+  std::cout.precision(17);
+  std::cout << name << " = [" << std::endl;
+  for (unsigned int i = 0; i < M; i++) {
+    std::cout << "[";
+    for (unsigned int j = 0; j < N; j++) {
+      if (matrix[i * N + j] >= 0) {
+        std::cout << "  " << matrix[i * N + j] << ",";
+      } else {
+        std::cout << " " << matrix[i * N + j] << ",";
+      }
+    }
+    std::cout << " ]," << std::endl;
+  }
+  std::cout << "]" << std::endl;
+  std::cout << std::endl;
+}
+
 template <typename T> T *transpose(T *A, unsigned int M) {
   T *B = new T[M * M];
   for (unsigned int i = 0; i < M; i++) {
@@ -77,6 +116,18 @@ 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) {
+  for(unsigned int i = 0; i < M; i++) {
+    for(unsigned int j = 0; j < P; j++) {
+      C[i * P + j] = 0;
+      for(unsigned int k = 0; k < N; k++) {
+        C[i * P + j] += A[i * N + k] * B[k * P + j];
+      }
+    }
+  }
+}
+
 template <typename T1, typename T2>
 T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
   T1 *C = new T1[M * M];

include/Woodbury.hpp

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

smvars.sh

@@ -3,6 +3,10 @@ unset THRESHOLD
 unset MKL
 ENV=$1
 
+echo
+echo "Sherman-Morrison-Woodbury parameters"
+echo "------------------------------------"
+
 ## Set Sherman-Morrison root dir
 PWD=$(pwd)
 SRCDIR=$(dirname $BASH_SOURCE)

src/SMWB.cpp

@@ -5,13 +5,16 @@
 
 // Sherman-Morrison-Woodbury kernel 1
 // WB2, WB3, SM2 mixing scheme 1
-void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
-         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);
+           bool ok;
+           ok = WB2(Slater_inv, Dim, Updates, Updates_index);
+           if (!ok) {
+             std::cerr << "Woodbury kernel failed!" << std::endl;
              SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
            }
+         }
 
 extern "C" {
 void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,

src/Woodbury.cpp

@@ -5,65 +5,126 @@
 // (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
+// D := V * S^{-1}, 2 x dim
+//
+// All matrices are stored in row-major order
+
 #include "Woodbury.hpp"
 #include "Helpers.hpp"
 
 // Woodbury 2x2 kernel:
-void WB2(double *Slater_inv, unsigned int Dim,
+bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
-         double *Updates, unsigned int *Updates_index) {
+         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};
+  unsigned int V[2 * Dim]{0}; // 2 x Dim matrix stored in row-major order
-  for (unsigned int i = 0; i < Dim; i++) {
+  V[Updates_index[0] - 1] = 1;
-    if (i == Updates_index[0] - 1) V[i] = 1;
+  V[Dim + Updates_index[1] - 1] = 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++) {
+  double C[2 * Dim]{0};
-      for (unsigned int j = 0; j < 2; j++) {
+  matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
-        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++) {
+  double B[4]{0};
-      for (unsigned int j = 0; j < 2; j++) {
+  matMul2(V, C, B, 2, Dim, 2);
-        for (unsigned int k = 0; k < Dim; k++) {
+  // Compute 1 + B
-          B[i * Dim + j] += (i == j) + V[i * Dim + k] * C[k * Dim + j];
+  B[0] += 1;
-        }
+  B[3] += 1;
-      }
-    }
 
-  // Invert B with explicit formula for 2x2 inversion
+  // Invert 1 + B with explicit formula for 2x2 inversion
   double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]);
+  double Binv[4]{0};
   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
+  // Check if determinant of inverted matrix is not zero
+  double det = B[0] * B[3] - B[1] * B[2];
+  if (std::fabs(det) < threshold()) {
+    std::cerr << "Determinant approached 0!" << std::endl;
+    return false;
+  }
 
+  // Compute (S + U * V)^{-1} with Woobury identity
+  double D[2 * Dim]{0};
+  matMul2(V, Slater_inv, D, 2, Dim, Dim);
+  double tmp[2 * Dim]{0};
+  matMul2(Binv, D, tmp, 2, 2, Dim);
+  double tmp2[Dim * Dim]{0};
+  matMul2(C, tmp, tmp2, Dim, 2, Dim);
+  for (unsigned int i = 0; i < Dim * Dim; i++) {
+    Slater_inv[i] -= tmp2[i];
+  }
+  return true;
 }
 
 // Woodbury 3x3 kernel
-void WB3(double *Slater_inv, unsigned int Dim,
+bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
-         double *Updates, unsigned int *Updates_index) {
+         unsigned int *Updates_index) {
   std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
   
   // Construct V from Updates_index
+  unsigned int V[3 * Dim]{0}; // 2 x Dim matrix stored in row-major order
+  V[Updates_index[0] - 1] = 1;
+  V[Dim + Updates_index[1] - 1] = 1;
+  V[2 * Dim + Updates_index[2] - 1] = 1;
 
-           // Compute B from U, V and Slater_inv
+  // Compute C
+  double C[3 * Dim]{0};
+  matMul2(Slater_inv, Updates, C, Dim, Dim, 3);
+  // Compute B
+  double B[9]{0};
+  matMul2(V, C, B, 3, Dim, 3);
+  // Compute 1 + B
+  B[0] += 1;
+  B[4] += 1;
+  B[8] += 1;
 
-           // Invert B with explicit formula for 3x3 inversion
+  double Binv[9];
+  Binv[0] = B[4] * B[8] - B[5] * B[7];
+  Binv[3] = B[5] * B[6] - B[3] * B[8];
+  Binv[6] = B[3] * B[7] - B[4] * B[6];
+
+  Binv[1] = B[2] * B[7] - B[1] * B[8];
+  Binv[4] = B[0] * B[8] - B[2] * B[6];
+  Binv[7] = B[1] * B[6] - B[0] * B[7];
+
+  Binv[2] = B[1] * B[5] - B[2] * B[4];
+  Binv[5] = B[2] * B[3] - B[0] * B[5];
+  Binv[8] = B[0] * B[4] - B[1] * B[3];
+
+  // Check if determinant of inverted matrix is not zero
+  // If so, exigt and return false.
+  double det;
+  det = B[0] * (B[4] * B[8] - B[5] * B[7]) -
+        B[1] * (B[3] * B[8] - B[5] * B[6]) + B[2] * (B[3] * B[7] - B[4] * B[6]);
+  if (std::fabs(det) < threshold()) {
+    std::cerr << "Determinant approached 0!" << std::endl;
+    return false;
+  }
 
   // Compute (S + U * V)^{-1} with Woobury identity
-           
+  double D[3 * Dim]{0};
+  matMul2(V, Slater_inv, D, 3, Dim, Dim);
+  double tmp[3 * Dim]{0};
+  matMul2(Binv, D, tmp, 3, 3, Dim);
+  double tmp2[Dim * Dim]{0};
+  matMul2(C, tmp, tmp2, Dim, 3, Dim);
+  for (unsigned int i = 0; i < Dim * Dim; i++) {
+    Slater_inv[i] -= tmp2[i];
+  }
+  return true;
+}
+
+extern "C" {
+bool WB2_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
+           unsigned int **Updates_index) {
+  WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
+}
+bool WB3_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
+           unsigned int **Updates_index) {
+  WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
+}
 }

todo.txt

@@ -1,10 +0,0 @@
--=={TODO}==-
-
-* Looking at the eigenvalues of the intermediate matrix after a problematic update is
-* plotting the ratio of the moduli of the maximum and the minimum eigenvalues of
-
-* Use valgrind to find the problem with MaponiA3S
-* Contact Claudia if she would be interested in sharing datasets from Champ to test with SM
-
-* Make a representative selection of update cycles for presentatoin puroses
-