Merge pull request #16 from pablooliveira/slagel

New implementations for Sherman Morisson

Makefile

@@ -1,11 +1,11 @@
 ## Compilers
 ARCH =
-CXX = icpc
+CXX = clang++
-FC = ifort
+FC = flang-7
 H5CXX = h5c++
 
 ## Compiler flags
-CXXFLAGS = -O0 -debug full -traceback
+CXXFLAGS = -O0 
 FFLAGS = $(CXXFLAGS)
 H5CXXFLAGS = $(CXXFLAGS) -fPIC
 FLIBS = -lstdc++

src/SM_Standard.cpp

@@ -4,15 +4,18 @@
 #include "Helpers.hpp"
 #include <iostream>
 
-void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+
+// Naïve Sherman Morrison
+void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
         double *Updates, unsigned int *Updates_index) {
 
+  std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
   double C[Dim];
   double D[Dim];
 
+  unsigned int l = 0;
   // For each update
-  for (unsigned int l = 0; l < N_updates; l++) {
+  while (l < N_updates) {
-
     // C = A^{-1} x U_l
     for (unsigned int i = 0; i < Dim; i++) {
       C[i] = 0;
@@ -23,8 +26,8 @@ void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
-    if (den < 1e-6) {
+    if (fabs(den) < 1e-6) {
-      std::cerr << "Breakdown condition triggered" << std::endl;
+      std::cerr << "Breakdown condition triggered at " << l << std::endl;
     }
     double iden = 1 / den;
 
@@ -40,5 +43,133 @@ void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
         Slater_inv[i * Dim + j] -= update;
       }
     }
+
+    l += 1;
   }
 }
+
+// Sherman Morrison, with J. Slagel splitting
+// http://hdl.handle.net/10919/52966
+void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+        double *Updates, unsigned int *Updates_index) {
+
+  std::cerr << "Called SM2 with updates " << N_updates << std::endl;
+  double C[Dim];
+  double D[Dim];
+
+  double later_updates[Dim*N_updates];
+  unsigned int later_index[N_updates];
+  unsigned int later = 0;
+
+  unsigned int l = 0;
+  // For each update
+  while (l < N_updates) {
+    // C = A^{-1} x U_l
+    for (unsigned int i = 0; i < Dim; i++) {
+      C[i] = 0;
+      for (unsigned int j = 0; j < Dim; j++) {
+        C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
+      }
+    }
+
+    // Denominator
+    double den = 1 + C[Updates_index[l] - 1];
+    if (fabs(den) < 1e-6) {
+      std::cerr << "Breakdown condition triggered at " << l << std::endl;
+
+      // U_l = U_l / 2 (do the split)
+      for (unsigned int j = 0; j < Dim; j++) {
+        later_updates[later*Dim+j] = Updates[l*Dim+j]/2.0;
+        C[j] /= 2.0;
+      }
+      later_index[later] = Updates_index[l];
+      later++;
+
+      den = 1+C[Updates_index[l]-1];
+    }
+    double iden = 1 / den;
+
+    // D = v^T x A^{-1}
+    for (unsigned int j = 0; j < Dim; j++) {
+      D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
+    }
+
+    // A^{-1} = A^{-1} - C x D / den
+    for (unsigned int i = 0; i < Dim; i++) {
+      for (unsigned int j = 0; j < Dim; j++) {
+        double update = C[i] * D[j] * iden;
+        Slater_inv[i * Dim + j] -= update;
+      }
+    }
+    l += 1;
+  }
+
+  if (later > 0) {
+    SM2(Slater_inv, Dim, later, later_updates, later_index);
+  }
+}
+
+// Sherman Morrison, leaving zero denominators for later
+void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+        double *Updates, unsigned int *Updates_index) {
+
+  std::cerr << "Called SM2 with updates " << N_updates << std::endl;
+  double C[Dim];
+  double D[Dim];
+
+  double later_updates[Dim*N_updates];
+  unsigned int later_index[N_updates];
+  unsigned int later = 0;
+
+  unsigned int l = 0;
+  // For each update
+  while (l < N_updates) {
+    // C = A^{-1} x U_l
+    for (unsigned int i = 0; i < Dim; i++) {
+      C[i] = 0;
+      for (unsigned int j = 0; j < Dim; j++) {
+        C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
+      }
+    }
+
+    // Denominator
+    double den = 1 + C[Updates_index[l] - 1];
+    if (fabs(den) < 1e-6) {
+      std::cerr << "Breakdown condition triggered at " << l << std::endl;
+
+      for (unsigned int j = 0; j < Dim; j++) {
+        later_updates[later*Dim+j] = Updates[l*Dim+j];
+      }
+      later_index[later] = Updates_index[l];
+      later++;
+      l += 1;
+      continue;
+    }
+    double iden = 1 / den;
+
+    // D = v^T x A^{-1}
+    for (unsigned int j = 0; j < Dim; j++) {
+      D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
+    }
+
+    // A^{-1} = A^{-1} - C x D / den
+    for (unsigned int i = 0; i < Dim; i++) {
+      for (unsigned int j = 0; j < Dim; j++) {
+        double update = C[i] * D[j] * iden;
+        Slater_inv[i * Dim + j] -= update;
+      }
+    }
+    l += 1;
+  }
+
+  if (later > 0) {
+    SM3(Slater_inv, Dim, later, later_updates, later_index);
+  }
+}
+
+
+
+void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+        double *Updates, unsigned int *Updates_index) {
+  SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
+}

tests/test_external_h5.cpp

@@ -1,116 +0,0 @@
-#include <iostream>
-#include <string>
-#include "hdf5/serial/hdf5.h"
-#include "hdf5/serial/H5Cpp.h"
-
-#include "SM_MaponiA3.hpp"
-#include "SM_Standard.hpp"
-#include "Helpers.hpp"
-
-using namespace H5;
-//#define DEBUG
-
-const H5std_string FILE_NAME( "datasets.hdf5" );
-
-void read_int(H5File file, std::string key, unsigned int * data) {
-  DataSet ds = file.openDataSet(key);
-  ds.read(data, PredType::STD_U32LE);
-  ds.close();
-}
-
-void read_double(H5File file, std::string key, double * data) {
-  DataSet ds = file.openDataSet(key);
-  ds.read(data, PredType::IEEE_F64LE);
-  ds.close();
-}
-
-int test_cycle(H5File file, int cycle) {
-
-  /* Read the data */
-
-  std::string group = "cycle_" + std::to_string(cycle);
-
-  unsigned int dim, nupdates, col, i, j;
-  read_int(file, group + "/slater_matrix_dim", &dim);
-  read_int(file, group + "/nupdates", &nupdates);
-
-  double * slater_matrix = new double[dim*dim];
-  read_double(file, group + "/slater_matrix", slater_matrix);
-
-  double * slater_inverse = new double[dim*dim];
-  read_double(file, group + "/slater_inverse", slater_inverse);
-  slater_inverse = transpose(slater_inverse, dim);
-
-  unsigned int * col_update_index = new unsigned int[nupdates];
-  read_int(file, group + "/col_update_index", col_update_index);
-
-  double * updates = new double[nupdates*dim];
-  read_double(file, group + "/updates", updates);
-
-  double * u = new double[nupdates*dim];
-
-
-  /* Test */
-#ifdef DEBUG
-  showMatrix(slater_matrix, dim, "OLD Slater");
-#endif
-
-#ifdef DEBUG
-  showMatrix(slater_inverse, dim, "OLD Inverse");
-#endif
-
-  for (j = 0; j < nupdates; j++) {
-    for (i = 0; i < dim; i++) {
-      col = col_update_index[j];
-      u[i + j*dim] = updates[i + j*dim] - slater_matrix[i*dim + (col - 1)];
-      slater_matrix[i*dim + (col - 1)] = updates[i + j*dim];
-    }
-  }
-
-  //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
-  SM(slater_inverse, dim, nupdates, u, col_update_index);
-
-#ifdef DEBUG
-  showMatrix(slater_matrix, dim, "NEW Slater");
-#endif
-
-#ifdef DEBUG
-  showMatrix(slater_inverse, dim, "NEW Inverse");
-#endif
-
-  double * res = new double[dim*dim] {0};
-  matMul(slater_matrix, slater_inverse, res, dim);
-  bool ok = is_identity(res, dim, 0.5e-4);
-
-#ifdef DEBUG
-  showMatrix(res, dim, "Result");
-#endif
-
-  delete [] res, updates, u, col_update_index,
-            slater_matrix, slater_inverse;
-
-  return ok;
-}
-
-int main(int argc, char **argv) {
-  if (argc != 2) {
-    std::cerr << "Execute from within 'datasets/'" << std::endl;
-    std::cerr << "usage: test_external_h5 <cycle>" << std::endl;
-    return 1;
-  }
-  int cycle = std::stoi(argv[1]);
-  H5File file(FILE_NAME, H5F_ACC_RDONLY);
-
-  bool ok = test_cycle(file, cycle);
-
-  if (ok) {
-    std::cerr << "ok -- cycle " << std::to_string(cycle)
-              << std::endl;
-  }
-  else {
-    std::cerr << "failed -- cycle " << std::to_string(cycle)
-              << std::endl;
-  }
-  return ok;
-}
-

tests/test_internal_h5.cpp

@@ -12,6 +12,18 @@ using namespace H5;
 
 const H5std_string FILE_NAME( "datasets.hdf5" );
 
+
+double residual2(double * A, unsigned int Dim) {
+  double res = 0.0;
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      double delta = (A[i * Dim + j] - (i == j));
+      res += delta*delta;
+    }
+  }
+  return res;
+}
+
 void read_int(H5File file, std::string key, unsigned int * data) {
   DataSet ds = file.openDataSet(key);
   ds.read(data, PredType::STD_U32LE);
@@ -66,8 +78,8 @@ int test_cycle(H5File file, int cycle) {
     }
   }
 
-  MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
+  //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
-  //SM(slater_inverse, dim, nupdates, updates, col_update_index);
+  SM(slater_inverse, dim, nupdates, u, col_update_index);
 
 #ifdef DEBUG
   showMatrix(slater_matrix, dim, "NEW Slater");
@@ -79,7 +91,10 @@ int test_cycle(H5File file, int cycle) {
 
   double * res = new double[dim*dim] {0};
   matMul(slater_matrix, slater_inverse, res, dim);
-  bool ok = is_identity(res, dim, 0.5e-4);
+  bool ok = is_identity(res, dim, 1e-4);
+
+  double res2 = residual2(res, dim);
+  std::cerr << "Residual = " << res2 << std::endl;
 
 #ifdef DEBUG
   showMatrix(res, dim, "Result");