diff --git a/.gitignore b/.gitignore index 6d31edd..4c5238b 100644 --- a/.gitignore +++ b/.gitignore @@ -4,5 +4,5 @@ cMaponiA3_test fMaponiA3_test QMCChem_dataset_test -Slater_* +Slater* Updates* diff --git a/Helpers.hpp b/Helpers.hpp index d0e6ac5..3c361de 100644 --- a/Helpers.hpp +++ b/Helpers.hpp @@ -37,9 +37,14 @@ template void showMatrix(T *matrix, unsigned int M, string name) { cout << name << " = " << endl; for (unsigned int i = 0; i < M; i++) { - cout << "[ "; + cout << "["; 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; } @@ -59,6 +64,17 @@ void showMatrixT(T **matrix, unsigned int size, string name) { cout << endl; } +template +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 T *matMul(T *A, T *B, unsigned int M) { T *C = new T[M*M]; @@ -72,7 +88,6 @@ T *matMul(T *A, T *B, unsigned int M) { return C; } - template T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) { T1 *C = new T1[M*M]; diff --git a/Makefile b/Makefile index 5c8d229..168bf18 100644 --- a/Makefile +++ b/Makefile @@ -1,56 +1,64 @@ +# ARCH = -xCORE-AVX2 + ## Used compilers H5CXX = h5c++ CXX = clang++ FC = flang -## Compiler flags -CXXFLAGS = -O0 -FFLAGS = -O0 +## Compiler flags & common obs & libs +CXXFLAGS = -O0 -debug full -traceback +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 -QMCChem_dataset_testDEP = QMCChem_dataset_test.f90 SM_MaponiA3_mod.f90 Utils_mod.f90 -QMCChem_dataset_testOBJ = SM_MaponiA3.o Utils_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o -QMCChem_dataset_testLIB = -lstdc++ +QMCChem_dataset_testOBJ = Utils_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o -## 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 $@ $< -## Compile recepies for Fortran fMaponiA3_test -%.o: %.f90 $(fMaponiA3_testDEP) +## Compile recepies for Fortran +%.o: %.f90 $(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 .PHONY: all clean distclean -all: cMaponiA3_test fMaponiA3_test QMCChem_dataset_test tests/test - clean: @rm -vf *.o *.mod 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 -cMaponiA3_test: $(cMaponiA3_testOBJ) +cMaponiA3_test_3x3_3: $(cMaponiA3_test_3x3_3OBJ) $(OBJS) $(CXX) $(ARCH) $(CXXFLAGS) -o $@ $^ -## Linking Fortran example program calling the C++ function 'Sherman_Morrison()' -fMaponiA3_test: $(fMaponiA3_testOBJ) - $(FC) $(ARCH) $(FFLAGS) $(fMaponiA3_testLIB) -o $@ $^ +## Linking Fortran example program calling the C++ function +fMaponiA3_test_3x3_3: $(fMaponiA3_test_3x3_3OBJ) $(OBJS) + $(FC) $(ARCH) $(FFLAGS) $(FLIBS) -o $@ $^ -## Linking Fortran example program calling the C++ function 'Sherman_Morrison()' -QMCChem_dataset_test: $(QMCChem_dataset_testOBJ) - $(FC) $(ARCH) $(FFLAGS) $(QMCChem_dataset_testLIB) -o $@ $^ +fMaponiA3_test_4x4_2: $(fMaponiA3_test_4x4_2OBJ) $(OBJS) + $(FC) $(ARCH) $(FFLAGS) $(FLIBS) -o $@ $^ -tests/test: tests/test.cpp SM_MaponiA3.o - $(H5CXX) $(ARCH) $(CXXFLAGS) -o $@ $^ +QMCChem_dataset_test: $(QMCChem_dataset_testOBJ) $(OBJS) + $(FC) $(ARCH) $(FFLAGS) $(FLIBS) -o $@ $^ diff --git a/QMCChem_dataset_test.f90 b/QMCChem_dataset_test.f90 index 4bebd0d..cc16320 100644 --- a/QMCChem_dataset_test.f90 +++ b/QMCChem_dataset_test.f90 @@ -7,7 +7,7 @@ program QMCChem_dataset_test integer :: i, j, col !! Iterators integer :: cycle_id, dim, n_updates 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", & cycle_id, & @@ -17,6 +17,13 @@ program QMCChem_dataset_test S_inv, & Updates_index, & 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 open(unit = 2000, file = "Slater_old.dat") @@ -24,7 +31,7 @@ program QMCChem_dataset_test do i=1,dim 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) + write(3000,"(E23.15, 1X)", advance="no") S_inv_trans(i,j) end do write(2000,*) write(3000,*) @@ -51,7 +58,7 @@ program QMCChem_dataset_test end do !! 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 open(unit = 4000, file = "Slater.dat") @@ -64,7 +71,6 @@ program QMCChem_dataset_test write(4000,*) write(5000,*) end do - close(4000) close(5000) diff --git a/QMCChem_dataset_test.m b/QMCChem_dataset_test.m index cc075f3..1d536a3 100644 --- a/QMCChem_dataset_test.m +++ b/QMCChem_dataset_test.m @@ -2,11 +2,11 @@ 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)) +det(Sold*Sold_inv) +trace(Sold*Sold_inv) +norm(Sold*Sold_inv) +det(S*S_inv) +trace(S*S_inv) +norm(S*S_inv) diff --git a/SM_MaponiA3.cpp b/SM_MaponiA3.cpp index 313de7d..73f4d6b 100644 --- a/SM_MaponiA3.cpp +++ b/SM_MaponiA3.cpp @@ -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 (i = 1; i < Dim + 1; i++) { for (j = 1; j < Dim + 1; j++) { - ylk[0][k][i] += Slater_inv[(i-1)+(j-1)*Dim] - * Updates[(k-1)*Dim+(j-1)]; + ylk[0][k][i] += Slater_inv[(i-1) + (j-1)*Dim] + * 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 for (l = 1; l < N_updates; l++) { for (j = l; j < N_updates + 1; j++) { component = Updates_index[p[j] - 1]; - cout << component << endl; breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]); } - showVector(breakdown, N_updates+1, "breakdown"); lbar = getMaxIndex(breakdown, N_updates + 1); - cout << "lbar = " << lbar << endl; - // Reset breakdown back to 0 for next round to avoid case where its first element is always the largest + // 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; } - cout << "l = " << l << endl; - cout << "p[l] = " << p[l] << endl; - cout << "p[lbar] = " << p[lbar] << endl; tmp = p[l]; p[l] = p[lbar]; p[lbar] = tmp; component = Updates_index[p[l] - 1]; - cout << "component = " << component << endl; beta = 1 + ylk[l - 1][p[l]][component]; if (beta == 0) { 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++) { alpha = ylk[l - 1][p[k]][component] / beta; - cout << "alpha = " << alpha << endl; 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]; } @@ -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. double *copy = Slater_inv; + Slater_inv = transpose(Slater_inv, Dim); + // Construct A-inverse from A0-inverse and the ylk - for (l = 0; l < N_updates; l++) { - k = l + 1; - component = Updates_index[p[k] - 1]; + for (l = 0; l < N_updates; l++) { // l = 0, 1 + k = l + 1; // k = 1, 2 + component = Updates_index[p[k] - 1]; // comp = 2, 4 beta = 1 + ylk[l][p[k]][component]; for (i = 0; i < Dim; i++) { 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 = transpose(Slater_inv, Dim); + // Assign the new values of 'Slater_inv' to the old values in 'copy[][]' for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { diff --git a/cMaponiA3_test.cpp b/cMaponiA3_test_3x3_3.cpp similarity index 100% rename from cMaponiA3_test.cpp rename to cMaponiA3_test_3x3_3.cpp diff --git a/fMaponiA3_test.f90.org b/fMaponiA3_test_3x3_3.f90 similarity index 100% rename from fMaponiA3_test.f90.org rename to fMaponiA3_test_3x3_3.f90 diff --git a/fMaponiA3_test_4x4_2 b/fMaponiA3_test_4x4_2 new file mode 100755 index 0000000..de1bd47 Binary files /dev/null and b/fMaponiA3_test_4x4_2 differ diff --git a/fMaponiA3_test.f90 b/fMaponiA3_test_4x4_2.f90 similarity index 64% rename from fMaponiA3_test.f90 rename to fMaponiA3_test_4x4_2.f90 index 16c425c..2ca2350 100644 --- a/fMaponiA3_test.f90 +++ b/fMaponiA3_test_4x4_2.f90 @@ -96,39 +96,19 @@ program Interface_test Updates_index(1) = 2 Updates_index(2) = 4 - write(*,*) - write(*,*) "Old S = " - do i=1,Dim - do j=1,Dim - write(*,"(F3.0,3X)", advance="no") S(i,j) + !! Write current S and S_inv to file for check in Octave + open(unit = 2000, file = "Slater_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) + write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) end do - write(*,*) + write(2000,*) + write(3000,*) end do - - 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(*,*) + close(2000) + close(3000) !! Update S do i=1,N_updates @@ -141,43 +121,19 @@ program Interface_test !! 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) + !! Write new S and S_inv to file for check in Octave + open(unit = 4000, file = "Slater.dat") + open(unit = 5000, file = "Slater_inv.dat") + do i=1,dim + do j=1,dim + write(4000,"(E23.15, 1X)", advance="no") S(i,j) + write(5000,"(E23.15, 1X)", advance="no") S_inv(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(*,*) + write(4000,*) + write(5000,*) end do + close(4000) + close(5000) deallocate(Updates_index, A, A_inv, S, Updates, S_inv) end program