TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2024-11-03 20:54:08 +01:00
Code Issues Projects Releases Wiki Activity

The algorithm now works for the following 4x4 example with 2 updates:

Browse Source 
S = [1,0,1,-1; 0,1,1,0; -1,0,-1,0; 1,1,1,1]
S_inv = [1,-1,1,1; 1,0,2,1; -1,1,-2,-1; -1,0,-1,0]
u1 = [0,-2,0,0]
u2 = [0,-1,0,0]
upd_idx = [2,4]

To go from Maponi's examples where the number of updates is always equal
to the the dimension of the matrix, and the decomposition is always
diagonal, to cases with a non-diagonal decomposition and a number of
updates unequal to its size, the following changed needed to be made:

* in the calculation of the {y0k} an extra inner for-loop needs to be
  added to make it a full matrix-vector multiplication due to the fact
  that A0 is not a diagonal matrix

* in some places the use of the update-order vector p needs
  the be replaced with that of upd_idx to make sure the correct
  component of the ylk is selected and the proper rank-1 matrices are
  constructed

* when a matrix is passed from Fortran to C++ with 2D adressing, it is
  passed in colum-major order. The passed matrix needs to be transposed
  before passing to C++. Doing this inside the algorithm will break
  compatibility with called from C/C++.
This commit is contained in:
François Coppens 2021-02-21 18:28:08 +01:00
parent 5b56a8e8a7
commit 8d63dd1701
10 changed files with 107 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

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

21
Helpers.hpp
View File

@ -37,9 +37,14 @@ template<typename T>
void showMatrix(T *matrix, unsigned int M, string name) { void showMatrix(T *matrix, unsigned int M, string name) {
cout << name << " = " << endl; cout << name << " = " << endl;
for (unsigned int i = 0; i < M; i++) { for (unsigned int i = 0; i < M; i++) {
cout << "[ "; cout << "[";
for (unsigned int j = 0; j < M; j++) { for (unsigned int j = 0; j < M; j++) {
cout << matrix[i*M+j] << " "; if (matrix[i*M+j] >= 0) {
cout << " " << matrix[i*M+j];
}
else {
cout << " " << matrix[i*M+j];
}
} }
cout << " ]" << endl; cout << " ]" << endl;
} }
@ -59,6 +64,17 @@ void showMatrixT(T **matrix, unsigned int size, string name) {
cout << endl; cout << 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++) {
for (unsigned int j = 0; j < M; j++) {
B[i*M + j] = A[i + j*M];
}
}
return B;
}
template<typename T> template<typename T>
T *matMul(T *A, T *B, unsigned int M) { T *matMul(T *A, T *B, unsigned int M) {
T *C = new T[M*M]; T *C = new T[M*M];
@ -72,7 +88,6 @@ T *matMul(T *A, T *B, unsigned int M) {
return C; return C;
} }
template<typename T1, typename T2> template<typename T1, typename T2>
T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) { T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
T1 *C = new T1[M*M]; T1 *C = new T1[M*M];

68
Makefile
View File

@ -1,56 +1,64 @@
# ARCH = -xCORE-AVX2
## Used compilers ## Used compilers
H5CXX = h5c++ H5CXX = h5c++
CXX = clang++ CXX = clang++
FC = flang FC = flang
## Compiler flags ## Compiler flags & common obs & libs
CXXFLAGS = -O0 CXXFLAGS = -O0 -debug full -traceback
FFLAGS = -O0 FFLAGS = -O0 -debug full -traceback
FLIBS = -lstdc++
OBJS = SM_MaponiA3.o
## Deps & objs for C++ cMaponiA3_test
cMaponiA3_testDEP = cMaponiA3_test.cpp SM_MaponiA3.cpp SM_MaponiA3.hpp Helpers.hpp
cMaponiA3_testOBJ = cMaponiA3_test.o SM_MaponiA3.o
## Deps & objs for Fortran fMaponiA3_test
fMaponiA3_testDEP = fMaponiA3_test.f90 SM_MaponiA3_mod.f90
fMaponiA3_testOBJ = SM_MaponiA3.o SM_MaponiA3_mod.o fMaponiA3_test.o
fMaponiA3_testLIB = -lstdc++
## Deps & objs for C++ cMaponiA3_test_3x3_3
cMaponiA3_test_3x3_3OBJ = cMaponiA3_test_3x3_3.o
## Deps & objs for Fortran fMaponiA3_test_3x3_3
fMaponiA3_test_3x3_3OBJ = SM_MaponiA3_mod.o fMaponiA3_test_3x3_3.o
## Deps & objs for Fortran fMaponiA3_test_4x4_2
fMaponiA3_test_4x4_2OBJ = SM_MaponiA3_mod.o fMaponiA3_test_4x4_2.o
## Deps & objs for Fortran QMCChem_dataset_test ## Deps & objs for Fortran QMCChem_dataset_test
QMCChem_dataset_testDEP = QMCChem_dataset_test.f90 SM_MaponiA3_mod.f90 Utils_mod.f90 QMCChem_dataset_testOBJ = Utils_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o
QMCChem_dataset_testOBJ = SM_MaponiA3.o Utils_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o
QMCChem_dataset_testLIB = -lstdc++
## Compile recipes for C++ cMaponiA3_test
%.o: %.cpp $(cMaponiA3_testDEP) ## Default build target: build everything
all: cMaponiA3_test_3x3_3 fMaponiA3_test_3x3_3 fMaponiA3_test_4x4_2 QMCChem_dataset_test
## Compile recipes for C++
%.o: %.cpp
$(CXX) $(ARCH) $(CXXFLAGS) -c -o $@ $< $(CXX) $(ARCH) $(CXXFLAGS) -c -o $@ $<
## Compile recepies for Fortran fMaponiA3_test ## Compile recepies for Fortran
%.o: %.f90 $(fMaponiA3_testDEP) %.o: %.f90
$(FC) $(ARCH) $(FFLAGS) -c -o $@ $< $(FC) $(ARCH) $(FFLAGS) -c -o $@ $<
## Explicit recipe to trigger rebuild and relinking when headerfile is changed
SM_MaponiA3.o: SM_MaponiA3.cpp Helpers.hpp
$(CXX) $(ARCH) $(CXXFLAGS) -c -o $@ $<
## Build tagets ## Build tagets
.PHONY: all clean distclean .PHONY: all clean distclean
all: cMaponiA3_test fMaponiA3_test QMCChem_dataset_test tests/test
clean: clean:
@rm -vf *.o *.mod @rm -vf *.o *.mod
distclean: clean distclean: clean
@rm -vf cMaponiA3_test fMaponiA3_test QMCChem_dataset_test Slater_* Updates.dat @rm -vf cMaponiA3_test_3x3_3 fMaponiA3_test_3x3_3 QMCChem_dataset_test Slater_* Updates.dat
## Linking the C++ example program ## Linking the C++ example program
cMaponiA3_test: $(cMaponiA3_testOBJ) cMaponiA3_test_3x3_3: $(cMaponiA3_test_3x3_3OBJ) $(OBJS)
$(CXX) $(ARCH) $(CXXFLAGS) -o $@ $^ $(CXX) $(ARCH) $(CXXFLAGS) -o $@ $^
## Linking Fortran example program calling the C++ function 'Sherman_Morrison()' ## Linking Fortran example program calling the C++ function
fMaponiA3_test: $(fMaponiA3_testOBJ) fMaponiA3_test_3x3_3: $(fMaponiA3_test_3x3_3OBJ) $(OBJS)
$(FC) $(ARCH) $(FFLAGS) $(fMaponiA3_testLIB) -o $@ $^ $(FC) $(ARCH) $(FFLAGS) $(FLIBS) -o $@ $^
## Linking Fortran example program calling the C++ function 'Sherman_Morrison()' fMaponiA3_test_4x4_2: $(fMaponiA3_test_4x4_2OBJ) $(OBJS)
QMCChem_dataset_test: $(QMCChem_dataset_testOBJ) $(FC) $(ARCH) $(FFLAGS) $(FLIBS) -o $@ $^
$(FC) $(ARCH) $(FFLAGS) $(QMCChem_dataset_testLIB) -o $@ $^
tests/test: tests/test.cpp SM_MaponiA3.o QMCChem_dataset_test: $(QMCChem_dataset_testOBJ) $(OBJS)
$(H5CXX) $(ARCH) $(CXXFLAGS) -o $@ $^ $(FC) $(ARCH) $(FFLAGS) $(FLIBS) -o $@ $^

14
QMCChem_dataset_test.f90
View File

@ -7,7 +7,7 @@ program QMCChem_dataset_test
integer :: i, j, col !! Iterators integer :: i, j, col !! Iterators
integer :: cycle_id, dim, n_updates integer :: cycle_id, dim, n_updates
integer(c_int), dimension(:), allocatable :: Updates_index integer(c_int), dimension(:), allocatable :: Updates_index
real(c_double), dimension(:,:), allocatable :: S, S_inv, Updates real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_trans, Updates
call Read_dataset("test.dataset.dat", & call Read_dataset("test.dataset.dat", &
cycle_id, & cycle_id, &
@ -17,6 +17,13 @@ program QMCChem_dataset_test
S_inv, & S_inv, &
Updates_index, & Updates_index, &
Updates) Updates)
allocate(S_inv_trans(dim,dim))
do i=1,dim
do j=1,dim
S_inv_trans(i,j) = S_inv(j,i)
end do
end do
!! Write current S and S_inv to file for check in Octave !! Write current S and S_inv to file for check in Octave
open(unit = 2000, file = "Slater_old.dat") open(unit = 2000, file = "Slater_old.dat")
@ -24,7 +31,7 @@ program QMCChem_dataset_test
do i=1,dim do i=1,dim
do j=1,dim do j=1,dim
write(2000,"(E23.15, 1X)", advance="no") S(i,j) write(2000,"(E23.15, 1X)", advance="no") S(i,j)
write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) write(3000,"(E23.15, 1X)", advance="no") S_inv_trans(i,j)
end do end do
write(2000,*) write(2000,*)
write(3000,*) write(3000,*)
@ -51,7 +58,7 @@ program QMCChem_dataset_test
end do end do
!! Send S_inv and Updates to MaponiA3 algo !! Send S_inv and Updates to MaponiA3 algo
call MaponiA3(S_inv, dim, n_updates, Updates, Updates_index) call MaponiA3(S_inv_trans, dim, n_updates, Updates, Updates_index)
!! Write new S and S_inv to file for check in Octave !! Write new S and S_inv to file for check in Octave
open(unit = 4000, file = "Slater.dat") open(unit = 4000, file = "Slater.dat")
@ -64,7 +71,6 @@ program QMCChem_dataset_test
write(4000,*) write(4000,*)
write(5000,*) write(5000,*)
end do end do
close(4000) close(4000)
close(5000) close(5000)

12
QMCChem_dataset_test.m
View File

@ -2,11 +2,11 @@ Sold = dlmread('Slater_old.dat');
Sold_inv = dlmread('Slater_old_inv.dat'); Sold_inv = dlmread('Slater_old_inv.dat');
S = dlmread('Slater.dat'); S = dlmread('Slater.dat');
S_inv = dlmread('Slater_inv.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)) det(Sold*Sold_inv)
trace(S*transpose(S_inv)) trace(Sold*Sold_inv)
norm(S*transpose(S_inv)) norm(Sold*Sold_inv)
det(S*S_inv)
trace(S*S_inv)
norm(S*S_inv)

45
SM_MaponiA3.cpp
View File

@ -30,52 +30,28 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
for (k = 1; k < N_updates + 1; k++) { for (k = 1; k < N_updates + 1; k++) {
for (i = 1; i < Dim + 1; i++) { for (i = 1; i < Dim + 1; i++) {
for (j = 1; j < Dim + 1; j++) { for (j = 1; j < Dim + 1; j++) {
ylk[0][k][i] += Slater_inv[(i-1)+(j-1)*Dim] ylk[0][k][i] += Slater_inv[(i-1) + (j-1)*Dim]
* Updates[(k-1)*Dim+(j-1)]; * Updates[(k-1)*Dim + (j-1)];
} }
} }
} }
cout << endl;
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;
cout << endl;
cout << "EVERYTHING LOOKS CORRECT UPTO THIS POINT" << endl;
cout << "NEXT STEP: COMP. YLKs FROM Y0Ks" << endl << endl;
showVector(p, N_updates + 1, "p");
showVector(Updates_index, N_updates, "Updates_index");
// Calculate all the ylk from the y0k // Calculate all the ylk from the y0k
for (l = 1; l < N_updates; l++) { for (l = 1; l < N_updates; l++) {
for (j = l; j < N_updates + 1; j++) { for (j = l; j < N_updates + 1; j++) {
component = Updates_index[p[j] - 1]; component = Updates_index[p[j] - 1];
cout << component << endl;
breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]); breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]);
} }
showVector(breakdown, N_updates+1, "breakdown");
lbar = getMaxIndex(breakdown, N_updates + 1); lbar = getMaxIndex(breakdown, N_updates + 1);
cout << "lbar = " << lbar << endl; // Reset breakdown back to 0 for next round to avoid case where its
// Reset breakdown back to 0 for next round to avoid case where its first element is always the largest // first element is always the largest
for (i = 0; i < N_updates + 1; i++) { for (i = 0; i < N_updates + 1; i++) {
breakdown[i] = 0; breakdown[i] = 0;
} }
cout << "l = " << l << endl;
cout << "p[l] = " << p[l] << endl;
cout << "p[lbar] = " << p[lbar] << endl;
tmp = p[l]; tmp = p[l];
p[l] = p[lbar]; p[l] = p[lbar];
p[lbar] = tmp; p[lbar] = tmp;
component = Updates_index[p[l] - 1]; component = Updates_index[p[l] - 1];
cout << "component = " << component << endl;
beta = 1 + ylk[l - 1][p[l]][component]; beta = 1 + ylk[l - 1][p[l]][component];
if (beta == 0) { if (beta == 0) {
cout << "Break-down occured. Exiting..." << endl; cout << "Break-down occured. Exiting..." << endl;
@ -83,7 +59,6 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
} }
for (k = l + 1; k < N_updates + 1; k++) { for (k = l + 1; k < N_updates + 1; k++) {
alpha = ylk[l - 1][p[k]][component] / beta; alpha = ylk[l - 1][p[k]][component] / beta;
cout << "alpha = " << alpha << endl;
for (i = 1; i < Dim + 1; i++) { 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]; ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - alpha * ylk[l - 1][p[l]][i];
} }
@ -95,19 +70,23 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
// a new pointer 'copy' that points to whereever 'Slater_inv' points to now. // a new pointer 'copy' that points to whereever 'Slater_inv' points to now.
double *copy = Slater_inv; double *copy = Slater_inv;
Slater_inv = transpose(Slater_inv, Dim);
// Construct A-inverse from A0-inverse and the ylk // Construct A-inverse from A0-inverse and the ylk
for (l = 0; l < N_updates; l++) { for (l = 0; l < N_updates; l++) { // l = 0, 1
k = l + 1; k = l + 1; // k = 1, 2
component = Updates_index[p[k] - 1]; component = Updates_index[p[k] - 1]; // comp = 2, 4
beta = 1 + ylk[l][p[k]][component]; beta = 1 + ylk[l][p[k]][component];
for (i = 0; i < Dim; i++) { for (i = 0; i < Dim; i++) {
for (j = 0; j < Dim; j++) { for (j = 0; j < Dim; j++) {
Al[i * Dim + j] = (i == j) - (j == p[k]-1) * ylk[l][p[k]][i + 1] / beta; Al[i * Dim + j] = (i == j) - (j == component - 1) * ylk[l][p[k]][i + 1] / beta;
} }
} }
Slater_inv = matMul(Al, Slater_inv, Dim); Slater_inv = matMul(Al, Slater_inv, Dim);
} }
Slater_inv = transpose(Slater_inv, Dim);
// Assign the new values of 'Slater_inv' to the old values in 'copy[][]' // Assign the new values of 'Slater_inv' to the old values in 'copy[][]'
for (i = 0; i < Dim; i++) { for (i = 0; i < Dim; i++) {
for (j = 0; j < Dim; j++) { for (j = 0; j < Dim; j++) {

0
cMaponiA3_test.cpp → cMaponiA3_test_3x3_3.cpp
View File

0
fMaponiA3_test.f90.org → fMaponiA3_test_3x3_3.f90
View File

BIN
fMaponiA3_test_4x4_2 Executable file
View File

Binary file not shown.

88
fMaponiA3_test.f90 → fMaponiA3_test_4x4_2.f90
View File

@ -96,39 +96,19 @@ program Interface_test
Updates_index(1) = 2 Updates_index(1) = 2
Updates_index(2) = 4 Updates_index(2) = 4
write(*,*) !! Write current S and S_inv to file for check in Octave
write(*,*) "Old S = " open(unit = 2000, file = "Slater_old.dat")
do i=1,Dim open(unit = 3000, file = "Slater_old_inv.dat")
do j=1,Dim do i=1,dim
write(*,"(F3.0,3X)", advance="no") S(i,j) do j=1,dim
write(2000,"(E23.15, 1X)", advance="no") S(i,j)
write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j)
end do end do
write(*,*) write(2000,*)
write(3000,*)
end do end do
close(2000)
write(*,*) close(3000)
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 !! Update S
do i=1,N_updates do i=1,N_updates
@ -141,43 +121,19 @@ program Interface_test
!! Update S_inv !! Update S_inv
call MaponiA3(S_inv, Dim, N_updates, Updates, Updates_index) call MaponiA3(S_inv, Dim, N_updates, Updates, Updates_index)
write(*,*) !! Write new S and S_inv to file for check in Octave
write(*,*) open(unit = 4000, file = "Slater.dat")
write(*,*) "New computed S = " open(unit = 5000, file = "Slater_inv.dat")
do i=1,Dim do i=1,dim
do j=1,Dim do j=1,dim
write(*,"(F3.0,3X)", advance="no") S(i,j) write(4000,"(E23.15, 1X)", advance="no") S(i,j)
write(5000,"(E23.15, 1X)", advance="no") S_inv(i,j)
end do end do
write(*,*) write(4000,*)
end do write(5000,*)
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 end do
close(4000)
close(5000)
deallocate(Updates_index, A, A_inv, S, Updates, S_inv) deallocate(Updates_index, A, A_inv, S, Updates, S_inv)
end program end program