mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-12-25 05:43:54 +01:00
Merge pull request #29 from fmgjcoppens/prep
Various small unrelated modifications
This commit is contained in:
François Coppens
GitHub
commit
bd9bd62fe9
26
Makefile
26
Makefile
@ -18,18 +18,24 @@ else ifeq ($(ENV),GNU)
|
||||
OPT = -O0
|
||||
DEBUG = -g
|
||||
else
|
||||
$(error No valid compiler environment set in $$ENV. First run: $$ source smvars.sh {intel | llvm | gnu})
|
||||
$(error No valid compiler environment set in $$ENV. \
|
||||
First run: $$ source smvars.sh {intel | llvm | gnu})
|
||||
endif
|
||||
CXXFLAGS = $(OPT) $(ARCH) $(DEBUG)
|
||||
CXXFLAGS = $(OPT) $(ARCH) $(DEBUG) -fPIC $(THRESHOLD)
|
||||
FFLAGS = $(CXXFLAGS)
|
||||
H5CXX = h5c++
|
||||
H5CXXFLAGS = $(CXXFLAGS) -fPIC
|
||||
H5CXXFLAGS = $(CXXFLAGS)
|
||||
FLIBS = -lstdc++
|
||||
|
||||
## Includes and dependencies
|
||||
INCLUDE = -I $(INC_DIR)/
|
||||
DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o $(OBJ_DIR)/SM_Helpers.o
|
||||
DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
|
||||
DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
|
||||
$(OBJ_DIR)/SM_MaponiA3S.o \
|
||||
$(OBJ_DIR)/SM_Standard.o \
|
||||
$(OBJ_DIR)/SM_Helpers.o
|
||||
DEPS_F = $(DEPS_CXX) \
|
||||
$(OBJ_DIR)/SM_MaponiA3_mod.o \
|
||||
$(OBJ_DIR)/Helpers_mod.o
|
||||
|
||||
## Directory structure
|
||||
SRC_DIR := src
|
||||
@ -67,7 +73,7 @@ $(OBJ_DIR)/%.o: $(TST_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
$(CXX) $(CXXFLAGS) $(INCLUDE) -c -o $@ $<
|
||||
|
||||
$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
$(CXX) $(CXXFLAGS) -fPIE $(INCLUDE) -c -o $@ $<
|
||||
$(CXX) $(CXXFLAGS) $(INCLUDE) -c -o $@ $<
|
||||
|
||||
## HDF5/C++ objects
|
||||
$(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
@ -85,14 +91,6 @@ endif
|
||||
$(OBJ_DIR)/%.o: $(TST_DIR)/%.f90 | $(OBJ_DIR)
|
||||
$(FC) $(FFLAGS) -I $(OBJ_DIR)/ -c -o $@ $<
|
||||
|
||||
### EXPLICIT BUILD RULES
|
||||
## special compiler flag -fPIC otherwise h5c++ builds fail
|
||||
$(OBJ_DIR)/SM_MaponiA3.o: $(SRC_DIR)/SM_MaponiA3.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
$(CXX) $(CXXFLAGS) -fPIC $(INCLUDE) -c -o $@ $<
|
||||
|
||||
$(OBJ_DIR)/SM_Standard.o: $(SRC_DIR)/SM_Standard.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
$(CXX) $(CXXFLAGS) -fPIC $(INCLUDE) -c -o $@ $<
|
||||
|
||||
|
||||
#### LINKING
|
||||
$(BIN_DIR)/cMaponiA3_test_3x3_3: $(OBJ_DIR)/cMaponiA3_test_3x3_3.o $(DEPS_CXX) | $(BIN_DIR)
|
||||
|
||||
@ -4,13 +4,25 @@ CXX = verificarlo-c++
|
||||
FC = verificarlo-f
|
||||
|
||||
## Compiler flags
|
||||
H5FLAGS = "-I/usr/include/hdf5/serial -L/usr/lib/x86_64-linux-gnu/hdf5/serial /usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5_hl_cpp.a /usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5_cpp.a /usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5_hl.a /usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5.a -lpthread -lsz -lz -ldl -lm -Wl,-rpath -Wl,/usr/lib/x86_64-linux-gnu/hdf5/serial"
|
||||
H5FLAGS = "-I/usr/include/hdf5/serial \
|
||||
-L/usr/lib/x86_64-linux-gnu/hdf5/serial \
|
||||
/usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5_hl_cpp.a \
|
||||
/usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5_cpp.a \
|
||||
/usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5_hl.a \
|
||||
/usr/lib/x86_64-linux-gnu/hdf5/serial/libhdf5.a \
|
||||
-lpthread -lsz -lz -ldl -lm \
|
||||
-Wl,-rpath -Wl,/usr/lib/x86_64-linux-gnu/hdf5/serial"
|
||||
CXXFLAGS = -O0 -g $(H5FLAGS)
|
||||
FFLAGS = -O0 -g $(H5FLAGS)
|
||||
FFLAGS = $(CXXFLAGS)
|
||||
|
||||
INCLUDE = -I $(INC_DIR)/
|
||||
DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o $(OBJ_DIR)/SM_Helpers.o
|
||||
DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
|
||||
DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
|
||||
$(OBJ_DIR)/SM_MaponiA3S.o \
|
||||
$(OBJ_DIR)/SM_Standard.o \
|
||||
$(OBJ_DIR)/SM_Helpers.o
|
||||
DEPS_F = $(DEPS_CXX) \
|
||||
$(OBJ_DIR)/SM_MaponiA3_mod.o \
|
||||
$(OBJ_DIR)/Helpers_mod.o
|
||||
FLIBS = -lstdc++
|
||||
|
||||
SRC_DIR := src
|
||||
@ -33,7 +45,7 @@ EXEC := $(BIN_DIR)/cMaponiA3_test_3x3_3 \
|
||||
all: $(EXEC)
|
||||
|
||||
clean:
|
||||
@rm -vrf $(OBJ_DIR)
|
||||
@rm -vrf $(OBJ_DIR) .*.o *.dbg *.cmdx *.cmod *.ilm *.stb
|
||||
|
||||
distclean: clean
|
||||
@rm -vrf $(BIN_DIR) \
|
||||
@ -50,7 +62,7 @@ $(OBJ_DIR)/%.o: $(TST_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
$(CXX) $(CXXFLAGS) $(ARCH) $(INCLUDE) -c -o $@ $<
|
||||
|
||||
$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
$(CXX) $(CXXFLAGS) -fPIE $(INCLUDE) -c -o $@ $<
|
||||
$(CXX) $(CXXFLAGS) $(INCLUDE) -c -o $@ $<
|
||||
|
||||
## HDF5/C++ objects
|
||||
$(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
@ -58,32 +70,19 @@ $(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR)
|
||||
|
||||
## Fortran modules
|
||||
$(OBJ_DIR)/%_mod.o: $(SRC_DIR)/%_mod.f90 | $(OBJ_DIR)
|
||||
$(FC) $(FFLAGS) $(ARCH) -module $(OBJ_DIR)/ -c -o $@ $<
|
||||
$(FC) $(FFLAGS) $(ARCH) -J $(OBJ_DIR)/ -c -o $@ $<
|
||||
|
||||
## Fortran objects
|
||||
$(OBJ_DIR)/%.o: $(TST_DIR)/%.f90 | $(OBJ_DIR)
|
||||
$(FC) $(FFLAGS) $(ARCH) -I $(OBJ_DIR)/ -c -o $@ $<
|
||||
|
||||
### EXPLICIT BUILD RULES
|
||||
## special compiler flag -fPIC otherwise h5c++ builds fail
|
||||
$(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
|
||||
$(BIN_DIR)/cMaponiA3_test_3x3_3: $(OBJ_DIR)/cMaponiA3_test_3x3_3.o $(DEPS_CXX) | $(BIN_DIR)
|
||||
$(CXX) -o $@ $^
|
||||
|
||||
#begin H5
|
||||
$(BIN_DIR)/test_h5: $(OBJ_DIR)/test_h5.o $(DEPS_CXX) | $(BIN_DIR)
|
||||
$(CXX) -o $@ $^ $(H5FLAGS)
|
||||
#$(BIN_DIR)/test_external_h5: $(OBJ_DIR)/test_external_h5.o $(DEPS_CXX) | $(BIN_DIR)
|
||||
# $(CXX) -o $@ $^ $(H5FLAGS)
|
||||
#end H5
|
||||
|
||||
|
||||
$(BIN_DIR)/fMaponiA3_test_3x3_3: $(DEPS_F) $(OBJ_DIR)/fMaponiA3_test_3x3_3.o | $(BIN_DIR)
|
||||
$(FC) $(FLIBS) -o $@ $^
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
Binary file not shown.
@ -535,12 +535,12 @@
|
||||
(04,04) 0.100000000000000E+01 0.000000000000000E+00
|
||||
#COL_UPDATE_INDEX: 2
|
||||
#COL_UPDATE_COMP_(01): 0.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(02): -0.200000000000000E+01
|
||||
#COL_UPDATE_COMP_(02): -1.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(03): 0.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(04): 0.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(04): 1.000000000000000E+00
|
||||
#COL_UPDATE_INDEX: 4
|
||||
#COL_UPDATE_COMP_(01): 0.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(02): -0.100000000000000E+01
|
||||
#COL_UPDATE_COMP_(01): -1.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(02): -1.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(03): 0.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(04): 0.000000000000000E+00
|
||||
#COL_UPDATE_COMP_(04): -1.000000000000000E+00
|
||||
#END_PACKET
|
||||
Binary file not shown.
@ -1,127 +1,129 @@
|
||||
// SM_Helpers.hpp
|
||||
// Some usefull helper functions to support the Maponi algorithm.
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <cstring>
|
||||
#include <cmath>
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
#ifndef THRESHOLD
|
||||
#define THRESHOLD 1e-3
|
||||
#endif
|
||||
double threshold();
|
||||
|
||||
void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
|
||||
|
||||
void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
||||
unsigned int *Updates_index, unsigned int *p,
|
||||
double ***ylk);
|
||||
void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
||||
unsigned int *Updates_index, unsigned int *p,
|
||||
double ***ylk);
|
||||
|
||||
template<typename T>
|
||||
void showScalar(T scalar, std::string name) {
|
||||
std::cout << name << " = " << scalar << std::endl << std::endl;
|
||||
template <typename T> void showScalar(T scalar, std::string name) {
|
||||
std::cout << name << " = " << scalar << std::endl << std::endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
template <typename T>
|
||||
void showVector(T *vector, unsigned int size, std::string name) {
|
||||
std::cout << name << " = " << std::endl;
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
std::cout << "[ " << vector[i] << " ]" << std::endl;
|
||||
}
|
||||
std::cout << std::endl;
|
||||
std::cout << name << " = " << std::endl;
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
std::cout << "[ " << vector[i] << " ]" << std::endl;
|
||||
}
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
template <typename T>
|
||||
void showMatrix(T *matrix, unsigned int M, std::string name) {
|
||||
std::cout.precision(17);
|
||||
std::cout << name << " = [" << std::endl;
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
std::cout << "[";
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
if (matrix[i*M + j] >= 0) {
|
||||
std::cout << " " << matrix[i*M + j] << ",";
|
||||
}
|
||||
else {
|
||||
std::cout << " " << matrix[i*M + j] << "," ;
|
||||
}
|
||||
}
|
||||
std::cout << " ]," << std::endl;
|
||||
std::cout.precision(17);
|
||||
std::cout << name << " = [" << std::endl;
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
std::cout << "[";
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
if (matrix[i * M + j] >= 0) {
|
||||
std::cout << " " << matrix[i * M + j] << ",";
|
||||
} else {
|
||||
std::cout << " " << matrix[i * M + j] << ",";
|
||||
}
|
||||
}
|
||||
std::cout << "]" << std::endl;
|
||||
std::cout << std::endl;
|
||||
std::cout << " ]," << std::endl;
|
||||
}
|
||||
std::cout << "]" << std::endl;
|
||||
std::cout << std::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];
|
||||
}
|
||||
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;
|
||||
}
|
||||
return B;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void matMul(T *A, T *B, T *C, unsigned int M) {
|
||||
memset(C, 0, M*M*sizeof(T));
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
for (unsigned int k = 0; k < M; k++) {
|
||||
C[i*M+j] += A[i*M+k] * B[k*M+j];
|
||||
}
|
||||
}
|
||||
template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
|
||||
memset(C, 0, M * M * sizeof(T));
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
for (unsigned int k = 0; k < M; k++) {
|
||||
C[i * M + j] += A[i * M + k] * B[k * M + j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename T1, typename T2>
|
||||
template <typename T1, typename T2>
|
||||
T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
|
||||
T1 *C = new T1[M*M];
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
C[i*M+j] = vec1[i+1] * vec2[j];
|
||||
}
|
||||
T1 *C = new T1[M * M];
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
C[i * M + j] = vec1[i + 1] * vec2[j];
|
||||
}
|
||||
return C;
|
||||
}
|
||||
return C;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T matDet(T **A, unsigned int M) {
|
||||
int det = 0, p, h, k, i, j;
|
||||
T **temp = new T*[M];
|
||||
for (int i = 0; i < M; i++) temp[i] = new T[M];
|
||||
if(M == 1) {
|
||||
return A[0][0];
|
||||
}
|
||||
else if(M == 2) {
|
||||
det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]);
|
||||
return det;
|
||||
}
|
||||
else {
|
||||
for(p = 0; p < M; p++) {
|
||||
h = 0;
|
||||
template <typename T> T matDet(T **A, unsigned int M) {
|
||||
int det = 0, p, h, k, i, j;
|
||||
T **temp = new T *[M];
|
||||
for (int i = 0; i < M; i++)
|
||||
temp[i] = new T[M];
|
||||
if (M == 1) {
|
||||
return A[0][0];
|
||||
} else if (M == 2) {
|
||||
det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]);
|
||||
return det;
|
||||
} else {
|
||||
for (p = 0; p < M; p++) {
|
||||
h = 0;
|
||||
k = 0;
|
||||
for (i = 1; i < M; i++) {
|
||||
for (j = 0; j < M; j++) {
|
||||
if (j == p) {
|
||||
continue;
|
||||
}
|
||||
temp[h][k] = A[i][j];
|
||||
k++;
|
||||
if (k == M - 1) {
|
||||
h++;
|
||||
k = 0;
|
||||
for(i = 1; i < M; i++) {
|
||||
for( j = 0; j < M; j++) {
|
||||
if(j == p) {
|
||||
continue;
|
||||
}
|
||||
temp[h][k] = A[i][j];
|
||||
k++;
|
||||
if(k == M-1) {
|
||||
h++;
|
||||
k = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
det = det + A[0][p] * pow(-1, p) * matDet(temp, M-1);
|
||||
}
|
||||
}
|
||||
return det;
|
||||
}
|
||||
det = det + A[0][p] * pow(-1, p) * matDet(temp, M - 1);
|
||||
}
|
||||
delete [] temp;
|
||||
return det;
|
||||
}
|
||||
delete[] temp;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool is_identity(T *A, unsigned int M, double tolerance) {
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
if (i==j && fabs(A[i*M+j]-1) > tolerance) return false;
|
||||
if (i!=j && fabs(A[i*M+j]) > tolerance) return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
template <typename T> bool is_identity(T *A, unsigned int M, double tolerance) {
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
if (i == j && fabs(A[i * M + j] - 1) > tolerance)
|
||||
return false;
|
||||
if (i != j && fabs(A[i * M + j]) > tolerance)
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
@ -1,3 +1,2 @@
|
||||
void MaponiA3(double *Slater_inv, unsigned int Dim,
|
||||
unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index);
|
||||
void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
|
||||
2
include/SM_MaponiA3S.hpp
Normal file
2
include/SM_MaponiA3S.hpp
Normal file
@ -0,0 +1,2 @@
|
||||
void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
@ -1,12 +1,12 @@
|
||||
// Naïve Sherman Morrison
|
||||
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
|
||||
// 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);
|
||||
double *Updates, unsigned int *Updates_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);
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
|
||||
@ -1,6 +1,7 @@
|
||||
#!/usr/bin/env bash
|
||||
|
||||
unset THRESHOLD
|
||||
ENV=$1
|
||||
THRESHOLD=$2
|
||||
|
||||
## Set Sherman-Morrison root dir
|
||||
PWD=$(pwd)
|
||||
@ -54,3 +55,9 @@ then
|
||||
export PATH=$SMROOT/bin:$PATH
|
||||
export SMVARS=true
|
||||
fi
|
||||
|
||||
## If a threshold is provided, export compiler flag
|
||||
if [[ $# -gt 1 ]]
|
||||
then
|
||||
export THRESHOLD="-DTHRESHOLD=$THRESHOLD"
|
||||
fi
|
||||
|
||||
@ -1,15 +1,24 @@
|
||||
#include "SM_Helpers.hpp"
|
||||
|
||||
// Set common break-down threshold
|
||||
double threshold() {
|
||||
const double threshold = THRESHOLD;
|
||||
#ifdef DEBUG
|
||||
std::cerr << "Break-down threshold set to: " << threshold << std::endl;
|
||||
#endif
|
||||
return threshold;
|
||||
}
|
||||
|
||||
void Switch(unsigned int *p, unsigned int l, unsigned int lbar) {
|
||||
unsigned int tmp = p[l+1];
|
||||
p[l+1] = p[lbar];
|
||||
unsigned int tmp = p[l + 1];
|
||||
p[l + 1] = p[lbar];
|
||||
p[lbar] = tmp;
|
||||
}
|
||||
|
||||
void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
||||
unsigned int *Updates_index, unsigned int *p,
|
||||
double ***ylk) {
|
||||
unsigned int lbar = l+1, max =0;
|
||||
void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
||||
unsigned int *Updates_index, unsigned int *p,
|
||||
double ***ylk) {
|
||||
unsigned int lbar = l + 1, max = 0;
|
||||
unsigned int index = 0, component = 0;
|
||||
unsigned int tmp = 0;
|
||||
double breakdown = 0;
|
||||
@ -17,11 +26,12 @@ void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
||||
index = p[j];
|
||||
component = Updates_index[index - 1];
|
||||
breakdown = abs(1 + ylk[l][index][component]);
|
||||
#ifdef DEBUG
|
||||
#ifdef DEBUG
|
||||
std::cout << "Inside selectLargestDenominator()" << std::endl;
|
||||
std::cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << std::endl;
|
||||
std::cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "]["
|
||||
<< component << "]) = " << breakdown << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
if (breakdown > max) {
|
||||
max = breakdown;
|
||||
lbar = j;
|
||||
|
||||
@ -6,19 +6,18 @@
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
void MaponiA3(double *Slater_inv, unsigned int Dim,
|
||||
unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index) {
|
||||
void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
|
||||
/*
|
||||
DECLARE AND INITIALISE ARRAYS
|
||||
*/
|
||||
|
||||
unsigned int k, l, i, j, component;
|
||||
unsigned int *p = new unsigned int[N_updates + 1] {0};
|
||||
unsigned int *p = new unsigned int[N_updates + 1]{0};
|
||||
double alpha, beta;
|
||||
double *Al = new double[Dim * Dim];
|
||||
double *next = new double[Dim*Dim] {0};
|
||||
double *next = new double[Dim * Dim]{0};
|
||||
double *last = Slater_inv, *tmp;
|
||||
|
||||
// Populate update-order vector
|
||||
@ -31,101 +30,99 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
|
||||
for (l = 0; l < N_updates; l++) {
|
||||
ylk[l] = new double *[N_updates + 1];
|
||||
for (k = 0; k < N_updates + 1; k++) {
|
||||
ylk[l][k] = new double[Dim + 1] {0};
|
||||
ylk[l][k] = new double[Dim + 1]{0};
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
START ALGORITHM
|
||||
*/
|
||||
|
||||
// Calculate the {y_{0,k}}
|
||||
for (k = 1; k < N_updates + 1; k++) {
|
||||
#ifdef DEBUG
|
||||
#ifdef DEBUG
|
||||
std::cout << "Compute y0k: " << std::endl;
|
||||
std::cout << "ylk[0][" << k << "][:]";
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
for (i = 1; i < Dim + 1; i++) {
|
||||
for (j = 1; j < Dim + 1; j++) {
|
||||
ylk[0][k][i] += Slater_inv[(i-1)*Dim + (j-1)]
|
||||
* Updates[(k-1)*Dim + (j-1)];
|
||||
ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] *
|
||||
Updates[(k - 1) * Dim + (j - 1)];
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG
|
||||
#ifdef DEBUG
|
||||
showVector(ylk[0][k], Dim, "");
|
||||
#endif
|
||||
#endif
|
||||
}
|
||||
|
||||
// Calculate the {y_{l,k}} from the {y_{0,k}}
|
||||
for (l = 0; l < N_updates; l++) {
|
||||
#ifdef DEBUG
|
||||
#ifdef DEBUG
|
||||
std::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
// For given l select intermediate update with largest break-down val
|
||||
selectLargestDenominator(l, N_updates, Updates_index, p, ylk);
|
||||
|
||||
// Select component and comp. bd-condition.
|
||||
component = Updates_index[p[l+1] - 1];
|
||||
beta = 1 + ylk[l][p[l+1]][component];
|
||||
#ifdef DEBUG
|
||||
std::cout << "p[l+1] = " << p[l+1] << std::endl;
|
||||
component = Updates_index[p[l + 1] - 1];
|
||||
beta = 1 + ylk[l][p[l + 1]][component];
|
||||
#ifdef DEBUG
|
||||
std::cout << "p[l+1] = " << p[l + 1] << std::endl;
|
||||
std::cout << "component = " << component << std::endl;
|
||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << std::endl;
|
||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
||||
<< "] = " << beta << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
if (fabs(beta) < 1e-3) {
|
||||
#endif
|
||||
if (fabs(beta) < threshold()) {
|
||||
std::cerr << "Break-down occured." << std::endl;
|
||||
}
|
||||
|
||||
// Compute intermediate update to Slater_inv
|
||||
#ifdef DEBUG
|
||||
// Compute intermediate update to Slater_inv
|
||||
#ifdef DEBUG
|
||||
std::cout << "Compute intermediate update to Slater_inv" << std::endl;
|
||||
std::cout << "component = " << component << std::endl;
|
||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << std::endl;
|
||||
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << std::endl;
|
||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
||||
<< "]" << std::endl;
|
||||
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
|
||||
<< std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
for (i = 0; i < Dim; i++) {
|
||||
for (j = 0; j < Dim; j++) {
|
||||
Al[i*Dim + j] = (i == j) - (j == component-1)
|
||||
* ylk[l][p[l+1]][i + 1] / beta;
|
||||
Al[i * Dim + j] =
|
||||
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
|
||||
}
|
||||
}
|
||||
matMul(Al, last, next, Dim);
|
||||
tmp = next;
|
||||
next = last;
|
||||
last = tmp;
|
||||
#ifdef DEBUG
|
||||
#ifdef DEBUG
|
||||
showMatrix(last, Dim, "last");
|
||||
#endif
|
||||
#endif
|
||||
|
||||
// For given l != 0 compute the next {y_{l,k}}
|
||||
for (k = l+2; k < N_updates+1; k++) {
|
||||
for (k = l + 2; k < N_updates + 1; k++) {
|
||||
alpha = ylk[l][p[k]][component] / beta;
|
||||
#ifdef DEBUG
|
||||
#ifdef DEBUG
|
||||
std::cout << "Inside k-loop: k = " << k << std::endl;
|
||||
std::cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << std::endl;
|
||||
std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
for (i = 1; i < Dim + 1; i++) {
|
||||
ylk[l+1][p[k]][i] = ylk[l][p[k]][i]
|
||||
- alpha * ylk[l][p[l+1]][i];
|
||||
ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i];
|
||||
}
|
||||
}
|
||||
}
|
||||
memcpy(Slater_inv, last, Dim*Dim*sizeof(double));
|
||||
|
||||
|
||||
memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
|
||||
|
||||
/*
|
||||
CLEANUP MEMORY
|
||||
*/
|
||||
|
||||
|
||||
for (l = 0; l < N_updates; l++) {
|
||||
for (k = 0; k < N_updates + 1; k++) {
|
||||
delete[] ylk[l][k];
|
||||
@ -136,11 +133,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
|
||||
}
|
||||
|
||||
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,
|
||||
*N_updates, *linUpdates,
|
||||
*Updates_index);
|
||||
}
|
||||
void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
}
|
||||
|
||||
153
src/SM_MaponiA3S.cpp
Normal file
153
src/SM_MaponiA3S.cpp
Normal file
@ -0,0 +1,153 @@
|
||||
// SM-MaponiA3_f.cpp
|
||||
// Algorithm 3 from P. Maponi,
|
||||
// p. 283, doi:10.1016/j.laa.2006.07.007
|
||||
#include "SM_MaponiA3S.hpp"
|
||||
#include "SM_Helpers.hpp"
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
/*
|
||||
DECLARE AND INITIALISE ARRAYS
|
||||
*/
|
||||
|
||||
unsigned int k, l, i, j, component;
|
||||
unsigned int *p = new unsigned int[N_updates + 1]{0};
|
||||
double alpha, beta;
|
||||
double *Al = new double[Dim * Dim];
|
||||
double *next = new double[Dim * Dim]{0};
|
||||
double *last = Slater_inv, *tmp;
|
||||
|
||||
// Populate update-order vector
|
||||
for (i = 0; i < N_updates; i++) {
|
||||
p[i + 1] = i + 1;
|
||||
}
|
||||
|
||||
// Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim]
|
||||
double ***ylk = new double **[N_updates];
|
||||
for (l = 0; l < N_updates; l++) {
|
||||
ylk[l] = new double *[N_updates + 1];
|
||||
for (k = 0; k < N_updates + 1; k++) {
|
||||
ylk[l][k] = new double[Dim + 1]{0};
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
START ALGORITHM
|
||||
*/
|
||||
|
||||
// Calculate the {y_{0,k}}
|
||||
for (k = 1; k < N_updates + 1; k++) {
|
||||
#ifdef DEBUG
|
||||
std::cout << "Compute y0k: " << std::endl;
|
||||
std::cout << "ylk[0][" << k << "][:]";
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
for (i = 1; i < Dim + 1; i++) {
|
||||
for (j = 1; j < Dim + 1; j++) {
|
||||
ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] *
|
||||
Updates[(k - 1) * Dim + (j - 1)];
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG
|
||||
showVector(ylk[0][k], Dim, "");
|
||||
#endif
|
||||
}
|
||||
|
||||
// Calculate the {y_{l,k}} from the {y_{0,k}}
|
||||
for (l = 0; l < N_updates; l++) {
|
||||
#ifdef DEBUG
|
||||
std::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
|
||||
// For given l select intermediate update with largest break-down val
|
||||
selectLargestDenominator(l, N_updates, Updates_index, p, ylk);
|
||||
|
||||
// Select component and comp. bd-condition.
|
||||
component = Updates_index[p[l + 1] - 1];
|
||||
beta = 1 + ylk[l][p[l + 1]][component];
|
||||
#ifdef DEBUG
|
||||
std::cout << "p[l+1] = " << p[l + 1] << std::endl;
|
||||
std::cout << "component = " << component << std::endl;
|
||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
||||
<< "] = " << beta << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
if (fabs(beta) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << component
|
||||
<< std::endl;
|
||||
|
||||
// for (unsigned int i = 0; i < Dim; i++) {
|
||||
// later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
|
||||
// ylk[l][p[l + 1]][i] /= 2.0;
|
||||
// }
|
||||
// later_index[later] = Updates_index[l];
|
||||
// later++;
|
||||
|
||||
// den = 1 + C[Updates_index[l] - 1];
|
||||
// }
|
||||
// double iden = 1 / den;
|
||||
|
||||
}
|
||||
|
||||
// Compute intermediate update to Slater_inv
|
||||
#ifdef DEBUG
|
||||
std::cout << "Compute intermediate update to Slater_inv" << std::endl;
|
||||
std::cout << "component = " << component << std::endl;
|
||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
||||
<< "]" << std::endl;
|
||||
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
|
||||
<< std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
for (i = 0; i < Dim; i++) {
|
||||
for (j = 0; j < Dim; j++) {
|
||||
Al[i * Dim + j] =
|
||||
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
|
||||
}
|
||||
}
|
||||
matMul(Al, last, next, Dim);
|
||||
tmp = next;
|
||||
next = last;
|
||||
last = tmp;
|
||||
#ifdef DEBUG
|
||||
showMatrix(last, Dim, "last");
|
||||
#endif
|
||||
|
||||
// For given l != 0 compute the next {y_{l,k}}
|
||||
for (k = l + 2; k < N_updates + 1; k++) {
|
||||
alpha = ylk[l][p[k]][component] / beta;
|
||||
#ifdef DEBUG
|
||||
std::cout << "Inside k-loop: k = " << k << std::endl;
|
||||
std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
for (i = 1; i < Dim + 1; i++) {
|
||||
ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i];
|
||||
}
|
||||
}
|
||||
}
|
||||
memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
|
||||
|
||||
/*
|
||||
CLEANUP MEMORY
|
||||
*/
|
||||
|
||||
for (l = 0; l < N_updates; l++) {
|
||||
for (k = 0; k < N_updates + 1; k++) {
|
||||
delete[] ylk[l][k];
|
||||
}
|
||||
delete[] ylk[l];
|
||||
}
|
||||
delete[] Al, next, p;
|
||||
}
|
||||
|
||||
extern "C" {
|
||||
void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
}
|
||||
@ -3,54 +3,40 @@
|
||||
#include "SM_Standard.hpp"
|
||||
#include "SM_Helpers.hpp"
|
||||
|
||||
// Set common break-down threshold
|
||||
double threshold() {
|
||||
double const threshold = 1e-3;
|
||||
return threshold;
|
||||
}
|
||||
|
||||
// Naïve Sherman Morrison
|
||||
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index)
|
||||
{
|
||||
|
||||
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
|
||||
while (l < N_updates)
|
||||
{
|
||||
while (l < N_updates) {
|
||||
// C = A^{-1} x U_l
|
||||
for (unsigned int i = 0; i < Dim; i++)
|
||||
{
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
C[i] = 0;
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
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) < threshold())
|
||||
{
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
|
||||
if (fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
}
|
||||
double iden = 1 / den;
|
||||
|
||||
// D = v^T x A^{-1}
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
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++)
|
||||
{
|
||||
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;
|
||||
}
|
||||
@ -63,9 +49,7 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
// 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)
|
||||
{
|
||||
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
std::cerr << "Called SM2 with updates " << N_updates << std::endl;
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
@ -76,29 +60,25 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
unsigned int l = 0;
|
||||
// For each update
|
||||
while (l < N_updates)
|
||||
{
|
||||
while (l < N_updates) {
|
||||
// C = A^{-1} x U_l
|
||||
for (unsigned int i = 0; i < Dim; i++)
|
||||
{
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
C[i] = 0;
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
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) < threshold())
|
||||
{
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
|
||||
if (fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[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;
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
|
||||
C[i] /= 2.0;
|
||||
}
|
||||
later_index[later] = Updates_index[l];
|
||||
later++;
|
||||
@ -108,16 +88,13 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double iden = 1 / den;
|
||||
|
||||
// D = v^T x A^{-1}
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
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++)
|
||||
{
|
||||
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;
|
||||
}
|
||||
@ -125,17 +102,14 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
l += 1;
|
||||
}
|
||||
|
||||
if (later > 0)
|
||||
{
|
||||
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)
|
||||
{
|
||||
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
std::cerr << "Called SM3 with updates " << N_updates << std::endl;
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
@ -146,26 +120,22 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
unsigned int l = 0;
|
||||
// For each update
|
||||
while (l < N_updates)
|
||||
{
|
||||
while (l < N_updates) {
|
||||
// C = A^{-1} x U_l
|
||||
for (unsigned int i = 0; i < Dim; i++)
|
||||
{
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
C[i] = 0;
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
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) < threshold())
|
||||
{
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
|
||||
if (fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
later_updates[later * Dim + j] = Updates[l * Dim + j];
|
||||
}
|
||||
later_index[later] = Updates_index[l];
|
||||
@ -176,16 +146,13 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double iden = 1 / den;
|
||||
|
||||
// D = v^T x A^{-1}
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
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++)
|
||||
{
|
||||
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;
|
||||
}
|
||||
@ -193,8 +160,13 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
l += 1;
|
||||
}
|
||||
|
||||
if (later > 0)
|
||||
{
|
||||
// If all the updates have failed, exit early with an error
|
||||
if (later == N_updates) {
|
||||
std::cerr << "SM3 cannot invert this matrix" << std::endl;
|
||||
return;
|
||||
}
|
||||
// If some have failed, make a recursive call
|
||||
else if (later > 0) {
|
||||
SM3(Slater_inv, Dim, later, later_updates, later_index);
|
||||
}
|
||||
}
|
||||
|
||||
@ -2,6 +2,7 @@
|
||||
#include "hdf5/serial/H5Cpp.h"
|
||||
|
||||
#include "SM_MaponiA3.hpp"
|
||||
#include "SM_MaponiA3S.hpp"
|
||||
#include "SM_Standard.hpp"
|
||||
#include "SM_Helpers.hpp"
|
||||
|
||||
@ -89,6 +90,8 @@ int test_cycle(H5File file, int cycle, std::string version) {
|
||||
|
||||
if (version == "maponia3") {
|
||||
MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "maponia3s") {
|
||||
MaponiA3S(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm1") {
|
||||
SM1(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm2") {
|
||||
|
||||
Loading…
Reference in New Issue
Block a user