Merge pull request #12 from pablooliveira/SMStandard

This commit is contained in:
François Coppens 2021-03-07 17:17:26 +01:00 committed by GitHub
commit f352a09815
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 72 additions and 18 deletions

View File

@ -1,16 +1,16 @@
## Compilers
ARCH = -xCORE-AVX2
ARCH =
H5CXX = h5c++
CXX = icpc
FC = ifort
CXX = clang++
FC = flang
## Compiler flags
H5CXXFLAGS = -O0 -g
CXXFLAGS = -O0 -g -traceback
FFLAGS = -O0 -g -traceback
CXXFLAGS = -O0 -g
FFLAGS = -O0 -g
INCLUDE = -I $(INC_DIR)/
DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o
DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o
DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
FLIBS = -lstdc++
@ -66,6 +66,8 @@ $(OBJ_DIR)/%.o: $(TST_DIR)/%.f90 | $(OBJ_DIR)
$(OBJ_DIR)/SM_MaponiA3.o: $(SRC_DIR)/SM_MaponiA3.cpp $(INC_DIR)/* | $(OBJ_DIR)
$(CXX) $(CXXFLAGS) -fPIC $(ARCH) $(INCLUDE) -c -o $@ $<
$(OBJ_DIR)/SM_Standard.o: $(SRC_DIR)/SM_Standard.cpp $(INC_DIR)/* | $(OBJ_DIR)
$(CXX) $(CXXFLAGS) -fPIC $(ARCH) $(INCLUDE) -c -o $@ $<
#### LINKING

View File

@ -35,20 +35,21 @@ void showVector(T *vector, unsigned int size, string name) {
template<typename T>
void showMatrix(T *matrix, unsigned int M, string name) {
cout << name << " = " << endl;
cout.precision(17);
cout << name << " = [" << endl;
for (unsigned int i = 0; i < M; i++) {
cout << "[";
for (unsigned int j = 0; j < M; j++) {
if (matrix[i*M + j] >= 0) {
cout << " " << matrix[i*M + j];
cout << " " << matrix[i*M + j] << ",";
}
else {
cout << " " << matrix[i*M + j];
cout << " " << matrix[i*M + j] << "," ;
}
}
cout << " ]" << endl;
cout << " ]," << endl;
}
cout << endl;
cout << "]" << endl;
}
template<typename T>

3
include/SM_Standard.hpp Normal file
View File

@ -0,0 +1,3 @@
void SM(double *Slater_inv, unsigned int Dim,
unsigned int N_updates, double *Updates,
unsigned int *Updates_index);

View File

@ -5,7 +5,7 @@
#include "Helpers.hpp"
void MaponiA3(double *Slater_inv, unsigned int Dim,
unsigned int N_updates, double *Updates,
unsigned int N_updates, double *Updates,
unsigned int *Updates_index) {
unsigned int k, l, lbar, i, j, tmp, component;
@ -27,7 +27,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
ylk[l][k] = new double[Dim + 1] {0};
}
}
// Calculate the y0k
for (k = 1; k < N_updates + 1; k++) {
for (i = 1; i < Dim + 1; i++) {
@ -55,9 +55,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
p[lbar] = tmp;
component = Updates_index[p[l] - 1];
beta = 1 + ylk[l - 1][p[l]][component];
if (beta == 0) {
if (fabs(beta) < 1e-6) {
cout << "Break-down occured. Exiting..." << endl;
exit;
exit(1);
}
for (k = l + 1; k < N_updates + 1; k++) {
alpha = ylk[l - 1][p[k]][component] / beta;
@ -110,7 +110,7 @@ extern "C" {
void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
unsigned int *N_updates, double **linUpdates,
unsigned int **Updates_index) {
MaponiA3(*linSlater_inv, *Dim,
MaponiA3(*linSlater_inv, *Dim,
*N_updates, *linUpdates,
*Updates_index);
}

44
src/SM_Standard.cpp Normal file
View File

@ -0,0 +1,44 @@
// SM-Standard.cpp
// Standard Sherman Morrison with multiple updates
#include "SM_Standard.hpp"
#include "Helpers.hpp"
#include <iostream>
void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index) {
double C[Dim];
double D[Dim];
// For each update
for (unsigned int l = 0; l < N_updates; l++) {
// 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 (den < 1e-6) {
std::cerr << "Breakdown condition triggered" << std::endl;
}
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;
}
}
}
}

View File

@ -4,10 +4,11 @@
#include "hdf5/serial/H5Cpp.h"
#include "SM_MaponiA3.hpp"
#include "SM_Standard.hpp"
#include "Helpers.hpp"
using namespace H5;
//#define DEBUG 0
//#define DEBUG
const H5std_string FILE_NAME( "datasets.hdf5" );
@ -62,7 +63,8 @@ int test_cycle(H5File file, int cycle) {
}
}
MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
//MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
SM(slater_inverse, dim, nupdates, updates, col_update_index);
#ifdef DEBUG
showMatrix(slater_matrix, dim, "NEW Slater");

View File

@ -4,6 +4,7 @@
#include "hdf5/serial/H5Cpp.h"
#include "SM_MaponiA3.hpp"
#include "SM_Standard.hpp"
#include "Helpers.hpp"
using namespace H5;
@ -63,6 +64,7 @@ int test_cycle(H5File file, int cycle) {
}
MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
//SM(slater_inverse, dim, nupdates, updates, col_update_index);
#ifdef DEBUG
showMatrix(slater_matrix, dim, "NEW Slater");