* Modified Maponi algo 3 so that it can be used with a number of updates

that is smaller than the size of the Slater-matrix
* Removed the Slater-matrix as an argument, since it is not used in the
  algo.
* Added a manual 4x4 example to debug MaponiA3 to work with a number of
  updates that is smaller than the size of the Slater-matrix
* Added a new Octave script to quickly check if the computes inverse is
  correct.

.gitignore

@@ -4,4 +4,5 @@
 cMaponiA3_test
 fMaponiA3_test
 QMCChem_dataset_test
-
+Slater_*
+Updates*

Makefile

@@ -38,7 +38,7 @@ clean:
 	@rm -vf *.o *.mod
 
 distclean: clean
-	@rm -vf cMaponiA3_test fMaponiA3_test QMCChem_dataset_test
+	@rm -vf cMaponiA3_test fMaponiA3_test QMCChem_dataset_test Slater_* Updates.dat
 
 ## Linking the C++ example program
 cMaponiA3_test: $(cMaponiA3_testOBJ)

Octave-tests/SMMaponiA3.m

@@ -2,12 +2,12 @@
 ## p. 283, doi:10.1016/j.laa.2006.07.007
 clc ## Clear the screen
 
-## Define the matrix to be inverted. This is example 8 from the paper
+### Define the matrix to be inverted. This is example 8 from the paper
-## In the future this matrix needs to be read from the function call arguments
+### In the future this matrix needs to be read from the function call arguments
-A=[1,1,-1; ...
+#A=[1,1,-1; ...
-   1,1,0; ...
+#   1,1,0; ...
-  -1,0,-1];
+#  -1,0,-1];
-A0=diag(diag(A));  ## The diagonal part of A
+#A0=diag(diag(A));  ## The diagonal part of A
 
 ### The modified example that gives all singular updates at some point  
 #A=[1,1,1; ...
@@ -17,10 +17,18 @@ A0=diag(diag(A));  ## The diagonal part of A
 
 ### A square uniform distributed random integer matrix with entries in [-1,1]  
 #do
-#  A=randi([-1,1],3,3);
+#  A=randi([-1,1],4,4);
 #  A0=diag(diag(A));  ## The diagonal part of A
 #until (det(A)!=0 && det(A0)!=0) ## We need both matrices to be simultaniously non-singular
 
+## Define the matrix to be inverted. This is example 8 from the paper
+## In the future this matrix needs to be read from the function call arguments
+A=[1,0, 1,-1; ...
+   0,-1, 1, -1; ...
+  -1,0,-1, 0; ...
+   1,1, 1, 1];
+A0=diag(diag(A));  ## The diagonal part of A
+
 ### A square uniform distributed random float matrix with entries in (0,1)  
 #do
 #  A=rand(5);
@@ -35,6 +43,7 @@ Ainv=zeros(nCols,nCols);
 ylk=zeros(nCols,nCols,nCols);
 p=zeros(nCols,1);
 breakdown=zeros(nCols,1);
+cutOff=10e-6
 
 A,A0
 printf("Determinant of A  is: %d\n",det(A))
@@ -109,4 +118,4 @@ else
     printf("Still not found. Giving up!\n");
     IdTest
   endif
-endif
+endif

QMCChem_dataset_test.f90

@@ -20,7 +20,7 @@ program QMCChem_dataset_test
 
     !! Write current S and S_inv to file for check in Octave
     open(unit = 2000, file = "Slater_old.dat")
-    open(unit = 3000, file = "Slater_inv_old.dat")
+    open(unit = 3000, file = "Slater_old_inv.dat")
     do i=1,dim
         do j=1,dim
             write(2000,"(E23.15, 1X)", advance="no") S(i,j)
@@ -29,38 +29,44 @@ program QMCChem_dataset_test
         write(2000,*)
         write(3000,*)
     end do
-    call flush(2000)
-    call flush(3000)
     close(2000)
     close(3000)
 
-    !! Send S, S_inv and Updates to MaponiA3 algo
+    !! Write Updates to file to check
-    call MaponiA3(S, S_inv, dim, n_updates, Updates, Updates_index)
+    open(unit = 2000, file = "Updates.dat")
+    do i=1,dim
+        do j=1,n_updates
+            write(2000,"(E23.15, 1X)", advance="no") Updates(i,j)
+        end do
+        write(2000,*)
+    end do
+    close(2000)
 
-    !! Update S itself
+    !! Update S
     do j=1,n_updates
         do i=1,dim
             col = Updates_index(j)
-            S(i,col) = S(i,col) + Updates(i,col)
+            S(i,col) = S(i,col) + Updates(i,j)
         end do
     end do
 
+    !! Send S_inv and Updates to MaponiA3 algo
+    call MaponiA3(S_inv, dim, n_updates, Updates, Updates_index)
+
     !! Write new S and S_inv to file for check in Octave
-    open(unit = 2000, file = "Slater_new.dat")
+    open(unit = 4000, file = "Slater.dat")
-    open(unit = 3000, file = "Slater_inv_new.dat")
+    open(unit = 5000, file = "Slater_inv.dat")
     do i=1,dim
         do j=1,dim
-            write(2000,"(E23.15, 1X)", advance="no") S(i,j)
+            write(4000,"(E23.15, 1X)", advance="no") S(i,j)
-            write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j)
+            write(5000,"(E23.15, 1X)", advance="no") S_inv(i,j)
         end do
-        write(2000,*)
+        write(4000,*)
-        write(3000,*)
+        write(5000,*)
     end do
-    call flush(2000)
-    call flush(3000)
 
-    close(2000)
+    close(4000)
-    close(3000)
+    close(5000)
 
     deallocate(S, S_inv, Updates, Updates_index)
 end program

QMCChem_dataset_test.m

@@ -0,0 +1,12 @@
+Sold = dlmread('Slater_old.dat');
+Sold_inv = dlmread('Slater_old_inv.dat');
+S = dlmread('Slater.dat');
+S_inv = dlmread('Slater_inv.dat');
+det(Sold*transpose(Sold_inv))
+trace(Sold*transpose(Sold_inv))
+norm(Sold*transpose(Sold_inv))
+
+det(S*transpose(S_inv))
+trace(S*transpose(S_inv))
+norm(S*transpose(S_inv))
+

SM_MaponiA3.cpp

@@ -4,54 +4,75 @@
 #include "SM_MaponiA3.hpp"
 #include "Helpers.hpp"
 
-void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M,
+void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-              unsigned int N_updates, double *Updates,
+              double *Updates, unsigned int *Updates_index) {
-              unsigned int *Updates_index) {
 
-  unsigned int k, l, lbar, i, j, tmp = M;
+  unsigned int k, l, lbar, i, j, tmp, component;
-  unsigned int *p = new unsigned int[M + 1];
+  unsigned int *p = new unsigned int[N_updates + 1];
   double alpha, beta;
-  double *breakdown = new double[M + 1];
+  double *breakdown = new double[N_updates + 1];
-  double *Al = new double[M * M];
+  double *Al = new double[Dim * Dim];
   p[0] = 0;
-  for (i = 0; i < M; i++) {
+  for (i = 0; i < N_updates; i++) {
     p[i + 1] = i + 1;
   }
 
   // Declare auxiliary solution matrix ylk
-  double ***ylk = new double **[M];
+  double ***ylk = new double **[N_updates];
-  for (l = 0; l < M; l++) {
+  for (l = 0; l < N_updates; l++) {
-    ylk[l] = new double *[M + 1];
+    ylk[l] = new double *[N_updates + 1];
-    for (k = 0; k < M + 1; k++) {
+    for (k = 0; k < N_updates + 1; k++) {
-      ylk[l][k] = new double[M + 1] {0};
+      ylk[l][k] = new double[Dim + 1] {0};
     }
   }
-
+  
-  // Calculate all the y0k in M^2 multiplications instead of M^3
+  cout << endl;
+  // Calculate all the y0k in Dim^2 multiplications instead of Dim^3
   for (k = 1; k < N_updates + 1; k++) {
-    for (i = 1; i < M + 1; i++) {
+    for (i = 1; i < Dim + 1; i++) {
+      cout << k << "," << i << endl;
       ylk[0][k][i] =
-          Slater_inv[(i - 1) * M + (i - 1)] * Updates[(i - 1) * M + (k - 1)];
+          Slater_inv[(i - 1) * Dim + (i - 1)] * Updates[(i - 1) * Dim + (k - 1)];
     }
   }
 
+  cout <<  "y0k = " << endl;
+  for (k = 1; k < N_updates + 1; k++) {
+      cout << "[ ";
+      for (i = 1; i < Dim + 1; i++) {
+          cout << ylk[0][k][i] << " ";
+      }
+      cout << " ]" << endl;
+  }
+  cout << endl;
+
+  showVector(p, N_updates + 1, "p");
+  showVector(Updates_index, N_updates, "Updates_index");
+ 
   // Calculate all the ylk from the y0k
-  for (l = 1; l < M; l++) {
+  for (l = 1; l < N_updates; l++) {
-    for (j = l; j < M + 1; j++) {
+    for (j = l; j < N_updates + 1; j++) {
-      breakdown[j] = abs(1 + ylk[l - 1][p[j]][p[j]]);
+      component = Updates_index[p[j] - 1];
+      cout << component << endl;
+      breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]);
+    }
+    lbar = getMaxIndex(breakdown, N_updates + 1);
+    // Reset breakdown back to 0 for next round to avoid case where its first element is always the largest
+    for (i = 0; i < N_updates + 1; i++) {
+      breakdown[i] = 0;
     }
-    lbar = getMaxIndex(breakdown, M + 1);
     tmp = p[l];
     p[l] = p[lbar];
     p[lbar] = tmp;
-    for (k = l + 1; k < M + 1; k++) {
+    component = Updates_index[p[l] - 1];
-      beta = 1 + ylk[l - 1][p[l]][p[l]];
+    beta = 1 + ylk[l - 1][p[l]][component];
-      if (beta == 0) {
+    if (beta == 0) {
-        cout << "Break-down condition occured. Exiting..." << endl;
+      cout << "Break-down occured. Exiting..." << endl;
-        exit;
+      exit;
-      }
+    }
-      for (i = 1; i < M + 1; i++) {
+    for (k = l + 1; k < N_updates + 1; k++) {
-        alpha = ylk[l - 1][p[k]][p[l]] / beta;
+      alpha = ylk[l - 1][p[k]][component] / beta;
+      for (i = 1; i < Dim + 1; i++) {
         ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - alpha * ylk[l - 1][p[l]][i];
       }
     }
@@ -63,26 +84,27 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M,
   double *copy = Slater_inv;
 
   // Construct A-inverse from A0-inverse and the ylk
-  for (l = 0; l < M; l++) {
+  for (l = 0; l < N_updates; l++) {
     k = l + 1;
-    beta = 1 + ylk[l][p[k]][p[k]];
+    component = Updates_index[p[k] - 1];
-    for (i = 0; i < M; i++) {
+    beta = 1 + ylk[l][p[k]][component];
-      for (j = 0; j < M; j++) {
+    for (i = 0; i < Dim; i++) {
-        Al[i * M + j] = (i == j) - (j == p[k]-1) * ylk[l][p[k]][i + 1] / beta;
+      for (j = 0; j < Dim; j++) {
+        Al[i * Dim + j] = (i == j) - (j == p[k]-1) * ylk[l][p[k]][i + 1] / beta;
       }
     }
-    Slater_inv = matMul(Al, Slater_inv, M);
+    Slater_inv = matMul(Al, Slater_inv, Dim);
   }
 
   // Assign the new values of 'Slater_inv' to the old values in 'copy[][]'
-  for (i = 0; i < M; i++) {
+  for (i = 0; i < Dim; i++) {
-    for (j = 0; j < M; j++) {
+    for (j = 0; j < Dim; j++) {
-      copy[i * M + j] = Slater_inv[i * M + j];
+      copy[i * Dim + j] = Slater_inv[i * Dim + j];
     }
   }
 
-  for (l = 0; l < M; l++) {
+  for (l = 0; l < N_updates; l++) {
-    for (k = 0; k < M + 1; k++) {
+    for (k = 0; k < N_updates + 1; k++) {
       delete[] ylk[l][k];
     }
     delete[] ylk[l];
@@ -92,11 +114,11 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M,
 }
 
 extern "C" {
-  void MaponiA3_f(double **linSlater0, double **linSlater_inv, unsigned int *Dim,
+  void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
                   unsigned int *N_updates, double **linUpdates,
                   unsigned int **Updates_index) {
-                    MaponiA3(*linSlater0, *linSlater_inv, *Dim, 
+    MaponiA3(*linSlater_inv, *Dim, 
-                             *N_updates, *linUpdates,
+             *N_updates, *linUpdates,
-                             *Updates_index);
+             *Updates_index);
   }
 }

SM_MaponiA3.hpp

@@ -1,3 +1,3 @@
-void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M,
+void MaponiA3(double *Slater_inv, unsigned int Dim,
               unsigned int N_updates, double *Updates,
               unsigned int *Updates_index);

SM_MaponiA3_mod.f90

@@ -1,10 +1,10 @@
 module Sherman_Morrison
     interface
-        subroutine MaponiA3(Slater0, Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_f")
+        subroutine MaponiA3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_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) :: Slater0, Updates
+            real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
             real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
         end subroutine MaponiA3
     end interface

cMaponiA3_test.cpp

@@ -48,7 +48,7 @@ int main() {
     // Define pointers dim and n_updates to use in Sherman-Morrison(...) function call
     unsigned int dim = M;
     unsigned int n_updates = M;
-    MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index);
+    MaponiA3(A0_inv, dim, n_updates, Ar, Ar_index);
     showMatrix(A0_inv, M, "A0_inv");
 
     // Deallocate all vectors and matrices

fMaponiA3_test.f90

@@ -3,50 +3,181 @@ program Interface_test
     use, intrinsic :: iso_c_binding, only : c_int, c_double
     implicit none
 
-    integer i, j !! Iterators
+    integer i, j, col !! Iterators
-    integer(c_int) :: dim, N_updates
+    integer(c_int) :: Dim, N_updates
-    integer(c_int), dimension(:), allocatable :: Ar_index
+    integer(c_int), dimension(:), allocatable :: Updates_index
-    real(c_double), dimension(:,:), allocatable :: A, A0, Ar
+    real(c_double), dimension(:,:), allocatable :: S, A, Updates
-    real(c_double), dimension(:,:), allocatable :: A0_inv
+    real(c_double), dimension(:,:), allocatable :: S_inv, A_inv
 
-    dim = 3
+    Dim = 4
-    N_updates = 3
+    N_updates = 2
-    allocate(Ar_index(dim), A(dim,dim), A0(dim,dim), Ar(dim,dim), A0_inv(dim,dim))
+    allocate( S(Dim,Dim), S_inv(Dim,Dim), A(Dim,Dim), A_inv(Dim,Dim), &
+              Updates(Dim,N_updates), Updates_index(N_updates))
+
+    !! Initialize S, S_inv, A and A_inv
+    S(1,1) = 1.0d0
+    S(1,2) = 0.0d0
+    S(1,3) = 1.0d0
+    S(1,4) = -1.0d0
+    S(2,1) = 0.0d0
+    S(2,2) = 1.0d0
+    S(2,3) = 1.0d0
+    S(2,4) = 0.0d0
+    S(3,1) = -1.0d0
+    S(3,2) = 0.0d0
+    S(3,3) = -1.0d0
+    S(3,4) = 0.0d0
+    S(4,1) = 1.0d0
+    S(4,2) = 1.0d0
+    S(4,3) = 1.0d0
+    S(4,4) = 1.0d0
+
+    S_inv(1,1) = 1.0d0
+    S_inv(1,2) = -1.0d0
+    S_inv(1,3) = 1.0d0
+    S_inv(1,4) = 1.0d0
+    S_inv(2,1) = 1.0d0
+    S_inv(2,2) = 0.0d0
+    S_inv(2,3) = 2.0d0
+    S_inv(2,4) = 1.0d0
+    S_inv(3,1) = -1.0d0
+    S_inv(3,2) = 1.0d0
+    S_inv(3,3) = -2.0d0
+    S_inv(3,4) = -1.0d0
+    S_inv(4,1) = -1.0d0
+    S_inv(4,2) = 0.0d0
+    S_inv(4,3) = -1.0d0
+    S_inv(4,4) = 0.0d0
 
-    !! Initialize A with M=3 and fill acc. to Eq. (17) from paper
     A(1,1) = 1.0d0
-    A(1,2) = 1.0d0
+    A(1,2) = 0.0d0
-    A(1,3) = -1.0d0
+    A(1,3) = 1.0d0
-    A(2,1) = 1.0d0
+    A(1,4) = -1.0d0
-    A(2,2) = 1.0d0
+    A(2,1) = 0.0d0
-    A(2,3) = 0.0d0
+    A(2,2) = -1.0d0
+    A(2,3) = 1.0d0
+    A(2,4) = -1.0d0
     A(3,1) = -1.0d0
     A(3,2) = 0.0d0
     A(3,3) = -1.0d0
+    A(3,4) = 0.0d0
+    A(4,1) = 1.0d0
+    A(4,2) = 1.0d0
+    A(4,3) = 1.0d0
+    A(4,4) = 1.0d0
 
-    !! Prepare the diagonal matrix A0 and the update matrix Ar
+    A_inv(1,1) = 0.0d0
-    do i=1,dim
+    A_inv(1,2) = -1.0d0
-        Ar_index(i) = i
+    A_inv(1,3) = -2.0d0
-        do j=1,dim
+    A_inv(1,4) = -1.0d0
-            if (i == j) then
+    A_inv(2,1) = 1.0d0
-                A0(i,j) = A(i,j)
+    A_inv(2,2) = 0.0d0
-                A0_inv(i,j) = 1.0d0 / A0(i,j)
+    A_inv(2,3) = 2.0d0
-            else
+    A_inv(2,4) = 1.0d0
-                A0(i,j) = 0.0d0
+    A_inv(3,1) = 0.0d0
-                A0_inv(i,j) = 0.0d0
+    A_inv(3,2) = 1.0d0
-            end if
+    A_inv(3,3) = 1.0d0
-            Ar(i,j) = A(i,j) - A0(i,j)
+    A_inv(3,4) = 1.0d0
-        end do
+    A_inv(4,1) = -1.0d0
-    end do
+    A_inv(4,2) = 0.0d0
+    A_inv(4,3) = -1.0d0
+    A_inv(4,4) = 0.0d0
 
-    call MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index)
+    !! Prepare Updates and Updates_index
+    Updates(1,1) = 0
+    Updates(1,2) = 0
+    Updates(2,1) = -2
+    Updates(2,2) = -1
+    Updates(3,1) = 0
+    Updates(3,2) = 0
+    Updates(4,1) = 0
+    Updates(4,2) = 0
 
-    do i=1,dim
+    Updates_index(1) = 2
-        do j=1,dim
+    Updates_index(2) = 4
-            write(*,"(F3.0,3X)", advance="no") A0_inv(i,j)
+
+    write(*,*)
+    write(*,*) "Old S = "
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") S(i,j)
         end do
         write(*,*)
     end do
 
-    deallocate(Ar_index, A, A0, Ar, A0_inv)
+    write(*,*)
+    write(*,*) "Old S_inv = "
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") S_inv(i,j)
+        end do
+        write(*,*)
+    end do
+
+    write(*,*)
+    write(*,*) "Updates = "
+    do i=1,Dim
+        do j=1,N_updates
+            write(*,"(F3.0,3X)", advance="no") Updates(i,j)
+        end do
+        write(*,*)
+    end do
+
+    write(*,*)
+    write(*,*) "Updates_index = "
+    do i=1,N_updates
+            write(*,"(I1,3X)", advance="no") Updates_index(i)
+    end do
+    write(*,*)
+
+    !! Update S
+    do i=1,N_updates
+        do j=1,Dim
+            col = Updates_index(i)
+            S(j,col) = S(j,col) + Updates(j,i)
+        end do
+    end do
+            
+    !! Update S_inv
+    call MaponiA3(S_inv, Dim, N_updates, Updates, Updates_index)
+
+    write(*,*)
+    write(*,*)
+    write(*,*) "New computed S = "
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") S(i,j)
+        end do
+        write(*,*)
+    end do
+
+    write(*,*)
+    write(*,*) "New actual S = "
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") A(i,j)
+        end do
+        write(*,*)
+    end do
+
+    write(*,*)
+    write(*,*)
+    write(*,*) "New computed S_inv = "
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") S_inv(i,j)
+        end do
+        write(*,*)
+    end do
+    
+    write(*,*)
+    write(*,*) "New actual S_inv = "
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") A_inv(i,j)
+        end do
+        write(*,*)
+    end do
+
+    deallocate(Updates_index, A, A_inv, S, Updates, S_inv)
 end program

fMaponiA3_test.f90.org

@@ -0,0 +1,52 @@
+program Interface_test
+    use Sherman_Morrison, only : MaponiA3
+    use, intrinsic :: iso_c_binding, only : c_int, c_double
+    implicit none
+
+    integer i, j !! Iterators
+    integer(c_int) :: Dim, N_updates
+    integer(c_int), dimension(:), allocatable :: Updates_index
+    real(c_double), dimension(:,:), allocatable :: A, S, Updates
+    real(c_double), dimension(:,:), allocatable :: S_inv
+
+    Dim = 3
+    N_updates = 3
+    allocate(Updates_index(Dim), A(Dim,Dim), S(Dim,Dim), Updates(Dim,Dim), S_inv(Dim,Dim))
+
+    !! Initialize A with M=3 and fill acc. to Eq. (17) from paper
+    A(1,1) = 1.0d0
+    A(1,2) = 1.0d0
+    A(1,3) = -1.0d0
+    A(2,1) = 1.0d0
+    A(2,2) = 1.0d0
+    A(2,3) = 0.0d0
+    A(3,1) = -1.0d0
+    A(3,2) = 0.0d0
+    A(3,3) = -1.0d0
+
+    !! Prepare the diagonal matrix S and the update matrix Updates
+    do i=1,Dim
+        Updates_index(i) = i
+        do j=1,Dim
+            if (i == j) then
+                S(i,j) = A(i,j)
+                S_inv(i,j) = 1.0d0 / S(i,j)
+            else
+                S(i,j) = 0.0d0
+                S_inv(i,j) = 0.0d0
+            end if
+            Updates(i,j) = A(i,j) - S(i,j)
+        end do
+    end do
+
+    call MaponiA3(S_inv, Dim, N_updates, Updates, Updates_index)
+
+    do i=1,Dim
+        do j=1,Dim
+            write(*,"(F3.0,3X)", advance="no") S_inv(i,j)
+        end do
+        write(*,*)
+    end do
+
+    deallocate(Updates_index, A, S, Updates, S_inv)
+end program