mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-08-25 14:41:45 +02:00
Merge pull request #31 from fmgjcoppens/add/slagel_to_ma3
Adding Slagel-splitting to Maponi algorithm 3
This commit is contained in:
commit
ff1f5655e5
2
Makefile
2
Makefile
@ -34,7 +34,7 @@ DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
|
|||||||
$(OBJ_DIR)/SM_Standard.o \
|
$(OBJ_DIR)/SM_Standard.o \
|
||||||
$(OBJ_DIR)/SM_Helpers.o
|
$(OBJ_DIR)/SM_Helpers.o
|
||||||
DEPS_F = $(DEPS_CXX) \
|
DEPS_F = $(DEPS_CXX) \
|
||||||
$(OBJ_DIR)/SM_MaponiA3_mod.o \
|
$(OBJ_DIR)/SM_mod.o \
|
||||||
$(OBJ_DIR)/Helpers_mod.o
|
$(OBJ_DIR)/Helpers_mod.o
|
||||||
|
|
||||||
## Directory structure
|
## Directory structure
|
||||||
|
@ -21,7 +21,7 @@ DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
|
|||||||
$(OBJ_DIR)/SM_Standard.o \
|
$(OBJ_DIR)/SM_Standard.o \
|
||||||
$(OBJ_DIR)/SM_Helpers.o
|
$(OBJ_DIR)/SM_Helpers.o
|
||||||
DEPS_F = $(DEPS_CXX) \
|
DEPS_F = $(DEPS_CXX) \
|
||||||
$(OBJ_DIR)/SM_MaponiA3_mod.o \
|
$(OBJ_DIR)/SM_mod.o \
|
||||||
$(OBJ_DIR)/Helpers_mod.o
|
$(OBJ_DIR)/Helpers_mod.o
|
||||||
FLIBS = -lstdc++
|
FLIBS = -lstdc++
|
||||||
|
|
||||||
|
@ -1,5 +1,5 @@
|
|||||||
#START_PACKET
|
#START_PACKET
|
||||||
#CYCLE_ID: 13
|
#CYCLE_ID: 1
|
||||||
#SLATER_MATRIX_DIM: 21
|
#SLATER_MATRIX_DIM: 21
|
||||||
#NUPDATES: 3
|
#NUPDATES: 3
|
||||||
#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
|
#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
|
||||||
@ -513,7 +513,7 @@
|
|||||||
#END_PACKET
|
#END_PACKET
|
||||||
|
|
||||||
#START_PACKET
|
#START_PACKET
|
||||||
#CYCLE_ID: 8169
|
#CYCLE_ID: 2
|
||||||
#SLATER_MATRIX_DIM: 4
|
#SLATER_MATRIX_DIM: 4
|
||||||
#NUPDATES: 2
|
#NUPDATES: 2
|
||||||
#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
|
#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
|
||||||
@ -544,3 +544,31 @@
|
|||||||
#COL_UPDATE_COMP_(03): 0.000000000000000E+00
|
#COL_UPDATE_COMP_(03): 0.000000000000000E+00
|
||||||
#COL_UPDATE_COMP_(04): -1.000000000000000E+00
|
#COL_UPDATE_COMP_(04): -1.000000000000000E+00
|
||||||
#END_PACKET
|
#END_PACKET
|
||||||
|
|
||||||
|
#START_PACKET
|
||||||
|
#CYCLE_ID: 3
|
||||||
|
#SLATER_MATRIX_DIM: 3
|
||||||
|
#NUPDATES: 3
|
||||||
|
#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
|
||||||
|
(01,01) 1.000000000000000E+00 1.000000000000000E+00
|
||||||
|
(01,02) 0.000000000000000E+00 0.000000000000000E+00
|
||||||
|
(01,03) 0.000000000000000E+00 0.000000000000000E+00
|
||||||
|
(02,01) 0.000000000000000E+00 0.000000000000000E+00
|
||||||
|
(02,02) 1.000000000000000E+00 1.000000000000000E+00
|
||||||
|
(02,03) 0.000000000000000E+00 0.000000000000000E+00
|
||||||
|
(03,01) 0.000000000000000E+00 0.000000000000000E+00
|
||||||
|
(03,02) 0.000000000000000E+00 0.000000000000000E+00
|
||||||
|
(03,03) -1.000000000000000E+00 -1.000000000000000E+00
|
||||||
|
#COL_UPDATE_INDEX: 1
|
||||||
|
#COL_UPDATE_COMP_(01): 1.000000000000000E+00
|
||||||
|
#COL_UPDATE_COMP_(02): 1.000000000000000E+00
|
||||||
|
#COL_UPDATE_COMP_(03): -1.000000000000000E+00
|
||||||
|
#COL_UPDATE_INDEX: 2
|
||||||
|
#COL_UPDATE_COMP_(01): 1.000000000000000E+00
|
||||||
|
#COL_UPDATE_COMP_(02): 1.000000000000000E+00
|
||||||
|
#COL_UPDATE_COMP_(03): 0.000000000000000E+00
|
||||||
|
#COL_UPDATE_INDEX: 3
|
||||||
|
#COL_UPDATE_COMP_(01): 1.000000000000000E+00
|
||||||
|
#COL_UPDATE_COMP_(02): 0.000000000000000E+00
|
||||||
|
#COL_UPDATE_COMP_(03): -1.000000000000000E+00
|
||||||
|
#END_PACKET
|
||||||
|
Binary file not shown.
@ -79,6 +79,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
if (fabs(beta) < threshold()) {
|
if (fabs(beta) < threshold()) {
|
||||||
std::cerr << "Break-down occured." << std::endl;
|
std::cerr << "Break-down occured." << std::endl;
|
||||||
}
|
}
|
||||||
|
double ibeta = 1.0 / beta;
|
||||||
|
|
||||||
// Compute intermediate update to Slater_inv
|
// Compute intermediate update to Slater_inv
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
@ -93,7 +94,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
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] =
|
Al[i * Dim + j] =
|
||||||
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
|
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
matMul(Al, last, next, Dim);
|
matMul(Al, last, next, Dim);
|
||||||
@ -106,7 +107,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
|
|
||||||
// For given l != 0 compute the next {y_{l,k}}
|
// 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;
|
alpha = ylk[l][p[k]][component] * ibeta;
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
std::cout << "Inside k-loop: k = " << k << std::endl;
|
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;
|
||||||
@ -129,13 +130,13 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
}
|
}
|
||||||
delete[] ylk[l];
|
delete[] ylk[l];
|
||||||
}
|
}
|
||||||
delete[] Al, next, p;
|
delete[] Al, next, p, ylk;
|
||||||
}
|
}
|
||||||
|
|
||||||
extern "C" {
|
extern "C" {
|
||||||
void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
|
void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
unsigned int *N_updates, double **linUpdates,
|
unsigned int *N_updates, double **linUpdates,
|
||||||
unsigned int **Updates_index) {
|
unsigned int **Updates_index) {
|
||||||
MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -11,7 +11,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
/*
|
/*
|
||||||
DECLARE AND INITIALISE ARRAYS
|
DECLARE AND INITIALISE ARRAYS
|
||||||
*/
|
*/
|
||||||
|
showVector(Updates_index, Dim, "Updates_index");
|
||||||
unsigned int k, l, i, j, component;
|
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 alpha, beta;
|
||||||
@ -19,6 +19,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
double *next = new double[Dim * Dim]{0};
|
double *next = new double[Dim * Dim]{0};
|
||||||
double *last = Slater_inv, *tmp;
|
double *last = Slater_inv, *tmp;
|
||||||
|
|
||||||
|
double later_updates[Dim * N_updates];
|
||||||
|
unsigned int later_index[N_updates];
|
||||||
|
unsigned int later = 0;
|
||||||
|
|
||||||
// Populate update-order vector
|
// Populate update-order vector
|
||||||
for (i = 0; i < N_updates; i++) {
|
for (i = 0; i < N_updates; i++) {
|
||||||
p[i + 1] = i + 1;
|
p[i + 1] = i + 1;
|
||||||
@ -36,13 +40,14 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
/*
|
/*
|
||||||
START ALGORITHM
|
START ALGORITHM
|
||||||
*/
|
*/
|
||||||
|
showMatrix(Slater_inv, Dim, "S^{-1}");
|
||||||
|
|
||||||
// Calculate the {y_{0,k}}
|
// Calculate the y_{0,k} = S_0^{-1} u_k
|
||||||
for (k = 1; k < N_updates + 1; k++) {
|
for (k = 1; k < N_updates + 1; k++) {
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
std::cout << "Compute y0k: " << std::endl;
|
std::cerr << "Compute y0k: " << std::endl;
|
||||||
std::cout << "ylk[0][" << k << "][:]";
|
std::cerr << "ylk[0][" << k << "][:]";
|
||||||
std::cout << std::endl;
|
std::cerr << std::endl;
|
||||||
#endif
|
#endif
|
||||||
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++) {
|
||||||
@ -51,15 +56,15 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
showVector(ylk[0][k], Dim, "");
|
showVector(ylk[0][k], Dim+1, "");
|
||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
|
|
||||||
// Calculate the {y_{l,k}} from the {y_{0,k}}
|
// Calculate the {y_{l,k}} from the {y_{0,k}}
|
||||||
for (l = 0; l < N_updates; l++) {
|
for (l = 0; l < N_updates; l++) {
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
std::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
|
std::cerr << "In outer compute-ylk-loop: l = " << l << std::endl;
|
||||||
std::cout << std::endl;
|
std::cerr << std::endl;
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
// For given l select intermediate update with largest break-down val
|
// For given l select intermediate update with largest break-down val
|
||||||
@ -69,64 +74,68 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
component = Updates_index[p[l + 1] - 1];
|
component = Updates_index[p[l + 1] - 1];
|
||||||
beta = 1 + ylk[l][p[l + 1]][component];
|
beta = 1 + ylk[l][p[l + 1]][component];
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
std::cout << "p[l+1] = " << p[l + 1] << std::endl;
|
std::cerr << "p[l+1] = " << p[l + 1] << std::endl;
|
||||||
std::cout << "component = " << component << std::endl;
|
std::cerr << "component = " << component << std::endl;
|
||||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
||||||
<< "] = " << beta << std::endl;
|
<< "] = " << beta << std::endl;
|
||||||
std::cout << std::endl;
|
std::cerr << std::endl;
|
||||||
#endif
|
#endif
|
||||||
if (fabs(beta) < threshold()) {
|
if (fabs(beta) < threshold()) {
|
||||||
std::cerr << "Breakdown condition triggered at " << component
|
#ifdef DEBUG
|
||||||
|
std::cerr << "Breakdown condition triggered at l = " << l << " and component "<< component
|
||||||
<< std::endl;
|
<< std::endl;
|
||||||
|
#endif
|
||||||
// for (unsigned int i = 0; i < Dim; i++) {
|
for (unsigned int i = 1; i < Dim + 1; i++) {
|
||||||
// later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
|
std::cout << "Let's split..." << std::endl;
|
||||||
// ylk[l][p[l + 1]][i] /= 2.0;
|
ylk[l][p[l + 1]][i] *= 0.5;
|
||||||
// }
|
// later_updates[later * Dim + i - 1] = ylk[l][p[l + 1]][i]; // Updates[l * Dim + i] * 0.5;
|
||||||
// later_index[later] = Updates_index[l];
|
later_updates[later * Dim + i - 1] = Updates[l * Dim + i - 1] * 0.5;
|
||||||
// later++;
|
}
|
||||||
|
later_index[later] = Updates_index[p[l + 1] - 1];
|
||||||
// den = 1 + C[Updates_index[l] - 1];
|
std::cout << "Later = " << later << std::endl;
|
||||||
// }
|
later++;
|
||||||
// double iden = 1 / den;
|
std::cout << "Later = " << later << std::endl;
|
||||||
|
beta = 1 + ylk[l][p[l + 1]][component];
|
||||||
|
std::cout << "New beta is: " << beta << std::endl;
|
||||||
}
|
}
|
||||||
|
double ibeta = 1.0 / beta;
|
||||||
|
showVector(ylk[l][p[l + 1]], Dim+1, "");
|
||||||
|
showVector(later_updates+(Dim*0), Dim, "later_updateds");
|
||||||
|
|
||||||
|
|
||||||
// Compute intermediate update to Slater_inv
|
// Compute intermediate update to Slater_inv
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
std::cout << "Compute intermediate update to Slater_inv" << std::endl;
|
std::cerr << "Compute intermediate update to Slater_inv" << std::endl;
|
||||||
std::cout << "component = " << component << std::endl;
|
std::cerr << "component = " << component << std::endl;
|
||||||
std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
|
||||||
<< "]" << std::endl;
|
<< "]" << std::endl;
|
||||||
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
|
std::cerr << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
|
||||||
<< std::endl;
|
<< std::endl;
|
||||||
std::cout << std::endl;
|
std::cerr << std::endl;
|
||||||
#endif
|
#endif
|
||||||
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] =
|
Al[i * Dim + j] =
|
||||||
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
|
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
matMul(Al, last, next, Dim);
|
matMul(Al, last, next, Dim); tmp = next; next = last; last = tmp;
|
||||||
tmp = next;
|
|
||||||
next = last;
|
|
||||||
last = tmp;
|
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
showMatrix(last, Dim, "last");
|
showMatrix(last, Dim, "last");
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
// For given l != 0 compute the next {y_{l,k}}
|
// 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;
|
alpha = ylk[l][p[k]][component] * ibeta;
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
std::cout << "Inside k-loop: k = " << k << std::endl;
|
std::cerr << "Inside k-loop: k = " << k << std::endl;
|
||||||
std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;
|
std::cerr << "ylk[" << l + 1 << "][" << p[k] << "][:]";
|
||||||
std::cout << std::endl;
|
std::cerr << std::endl;
|
||||||
#endif
|
#endif
|
||||||
for (i = 1; i < Dim + 1; i++) {
|
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];
|
||||||
}
|
}
|
||||||
|
showVector(ylk[l + 1][p[k]], Dim+1, "");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
|
memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
|
||||||
@ -141,13 +150,19 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
}
|
}
|
||||||
delete[] ylk[l];
|
delete[] ylk[l];
|
||||||
}
|
}
|
||||||
delete[] Al, next, p;
|
delete[] Al, next, p, ylk;
|
||||||
|
|
||||||
|
if (later > 0) {
|
||||||
|
std::cerr << "Entering recursive loop with " << later << " updates" << std::endl;
|
||||||
|
std::cout << later_index[0] << std::endl;
|
||||||
|
MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
extern "C" {
|
extern "C" {
|
||||||
void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
|
void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
unsigned int *N_updates, double **linUpdates,
|
unsigned int *N_updates, double **linUpdates,
|
||||||
unsigned int **Updates_index) {
|
unsigned int **Updates_index) {
|
||||||
MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -1,12 +0,0 @@
|
|||||||
module Sherman_Morrison
|
|
||||||
use, intrinsic :: iso_c_binding, only : c_int, c_double
|
|
||||||
interface
|
|
||||||
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) :: Updates
|
|
||||||
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
|
||||||
end subroutine MaponiA3
|
|
||||||
end interface
|
|
||||||
end module Sherman_Morrison
|
|
@ -162,7 +162,8 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
|
|
||||||
// If all the updates have failed, exit early with an error
|
// If all the updates have failed, exit early with an error
|
||||||
if (later == N_updates) {
|
if (later == N_updates) {
|
||||||
std::cerr << "SM3 cannot invert this matrix" << std::endl;
|
std::cerr << "SM3 cannot invert matrix." << std::endl;
|
||||||
|
showMatrix(Slater_inv, Dim, "Slater_inverse");
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
// If some have failed, make a recursive call
|
// If some have failed, make a recursive call
|
||||||
@ -170,3 +171,23 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
SM3(Slater_inv, Dim, later, later_updates, later_index);
|
SM3(Slater_inv, Dim, later, later_updates, later_index);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
extern "C" {
|
||||||
|
void SM1_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
|
unsigned int *N_updates, double **linUpdates,
|
||||||
|
unsigned int **Updates_index) {
|
||||||
|
SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
|
}
|
||||||
|
|
||||||
|
void SM2_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
|
unsigned int *N_updates, double **linUpdates,
|
||||||
|
unsigned int **Updates_index) {
|
||||||
|
SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
|
}
|
||||||
|
|
||||||
|
void SM3_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
|
unsigned int *N_updates, double **linUpdates,
|
||||||
|
unsigned int **Updates_index) {
|
||||||
|
SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
40
src/SM_mod.f90
Normal file
40
src/SM_mod.f90
Normal file
@ -0,0 +1,40 @@
|
|||||||
|
module Sherman_Morrison
|
||||||
|
use, intrinsic :: iso_c_binding, only : c_int, c_double
|
||||||
|
interface
|
||||||
|
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) :: Updates
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
|
end subroutine MaponiA3
|
||||||
|
subroutine MaponiA3S(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3S_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) :: Updates
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
|
end subroutine MaponiA3S
|
||||||
|
subroutine SM1(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM1_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) :: Updates
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
|
end subroutine SM1
|
||||||
|
subroutine SM2(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM2_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) :: Updates
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
|
end subroutine SM2
|
||||||
|
subroutine SM3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM3_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) :: Updates
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
|
end subroutine SM3
|
||||||
|
end interface
|
||||||
|
end module Sherman_Morrison
|
@ -57,12 +57,12 @@ program QMCChem_dataset_test
|
|||||||
!! S_inv needs to be transposed first before it
|
!! S_inv needs to be transposed first before it
|
||||||
!! goes to MaponiA3
|
!! goes to MaponiA3
|
||||||
call Transpose(S_inv, S_inv_t, dim)
|
call Transpose(S_inv, S_inv_t, dim)
|
||||||
call MaponiA3(S_inv_t, dim, n_updates, U, Updates_index)
|
! call MaponiA3S(S_inv_t, dim, n_updates, U, Updates_index)
|
||||||
|
call SM2(S_inv_t, dim, n_updates, U, Updates_index)
|
||||||
!! S_inv_t needs to be transposed back before it
|
!! S_inv_t needs to be transposed back before it
|
||||||
!! can be multiplied with S to test unity
|
!! can be multiplied with S to test unity
|
||||||
call Transpose(S_inv_t, S_inv, dim)
|
call Transpose(S_inv_t, S_inv, dim)
|
||||||
|
|
||||||
|
|
||||||
!! 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")
|
||||||
open(unit = 5000, file = "Slater_inv.dat")
|
open(unit = 5000, file = "Slater_inv.dat")
|
||||||
|
@ -12,6 +12,7 @@ ALGO=$1
|
|||||||
START=$2
|
START=$2
|
||||||
STOP=$3
|
STOP=$3
|
||||||
BLOCKSIZE=329
|
BLOCKSIZE=329
|
||||||
|
TOLERANCE=1e-3
|
||||||
|
|
||||||
BLOCKS=$((STOP / BLOCKSIZE))
|
BLOCKS=$((STOP / BLOCKSIZE))
|
||||||
LEFT=$((STOP % BLOCKSIZE))
|
LEFT=$((STOP % BLOCKSIZE))
|
||||||
@ -25,17 +26,17 @@ then
|
|||||||
do
|
do
|
||||||
BSTART=$(((i-1)*BLOCKSIZE+1))
|
BSTART=$(((i-1)*BLOCKSIZE+1))
|
||||||
BSTOP=$((i*BLOCKSIZE-1))
|
BSTOP=$((i*BLOCKSIZE-1))
|
||||||
$TEST $ALGO $BSTART $BSTOP
|
$TEST $ALGO $BSTART $BSTOP $TOLERANCE
|
||||||
LAST=$i
|
LAST=$i
|
||||||
done
|
done
|
||||||
LSTART=$((LAST*BLOCKSIZE+1))
|
LSTART=$((LAST*BLOCKSIZE+1))
|
||||||
LSTOP=$((LSTART+LEFT-1))
|
LSTOP=$((LSTART+LEFT-1))
|
||||||
$TEST $ALGO $LSTART $LSTOP
|
$TEST $ALGO $LSTART $LSTOP $TOLERANCE
|
||||||
else
|
else
|
||||||
for ((i=1; i<=$BLOCKS; i++))
|
for ((i=1; i<=$BLOCKS; i++))
|
||||||
do
|
do
|
||||||
BSTART=$(((i-1)*BLOCKSIZE+1))
|
BSTART=$(((i-1)*BLOCKSIZE+1))
|
||||||
BSTOP=$((i*BLOCKSIZE-1))
|
BSTOP=$((i*BLOCKSIZE-1))
|
||||||
$TEST $ALGO $BSTART $BSTOP
|
$TEST $ALGO $BSTART $BSTOP $TOLERANCE
|
||||||
done
|
done
|
||||||
fi
|
fi
|
||||||
|
@ -7,7 +7,7 @@
|
|||||||
#include "SM_Helpers.hpp"
|
#include "SM_Helpers.hpp"
|
||||||
|
|
||||||
using namespace H5;
|
using namespace H5;
|
||||||
// #define DEBUG
|
#define DEBUG
|
||||||
|
|
||||||
const H5std_string FILE_NAME( "dataset.hdf5" );
|
const H5std_string FILE_NAME( "dataset.hdf5" );
|
||||||
|
|
||||||
@ -46,7 +46,7 @@ void read_double(H5File file, std::string key, double * data) {
|
|||||||
ds.close();
|
ds.close();
|
||||||
}
|
}
|
||||||
|
|
||||||
int test_cycle(H5File file, int cycle, std::string version) {
|
int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||||
|
|
||||||
/* Read the data */
|
/* Read the data */
|
||||||
|
|
||||||
@ -61,7 +61,6 @@ int test_cycle(H5File file, int cycle, std::string version) {
|
|||||||
|
|
||||||
double * slater_inverse = new double[dim*dim];
|
double * slater_inverse = new double[dim*dim];
|
||||||
read_double(file, group + "/slater_inverse", slater_inverse);
|
read_double(file, group + "/slater_inverse", slater_inverse);
|
||||||
//slater_inverse = transpose(slater_inverse, dim);
|
|
||||||
|
|
||||||
unsigned int * col_update_index = new unsigned int[nupdates];
|
unsigned int * col_update_index = new unsigned int[nupdates];
|
||||||
read_int(file, group + "/col_update_index", col_update_index);
|
read_int(file, group + "/col_update_index", col_update_index);
|
||||||
@ -72,14 +71,11 @@ int test_cycle(H5File file, int cycle, std::string version) {
|
|||||||
double * u = new double[nupdates*dim];
|
double * u = new double[nupdates*dim];
|
||||||
|
|
||||||
/* Test */
|
/* Test */
|
||||||
#ifdef DEBUG
|
|
||||||
showMatrix(slater_matrix, dim, "OLD Slater");
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
showMatrix(slater_inverse, dim, "OLD Inverse");
|
showMatrix(slater_inverse, dim, "OLD Inverse");
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
|
// Transform replacement updates in 'updates[]' into additive updates in 'u[]'
|
||||||
for (j = 0; j < nupdates; j++) {
|
for (j = 0; j < nupdates; j++) {
|
||||||
for (i = 0; i < dim; i++) {
|
for (i = 0; i < dim; i++) {
|
||||||
col = col_update_index[j];
|
col = col_update_index[j];
|
||||||
@ -88,6 +84,16 @@ int test_cycle(H5File file, int cycle, std::string version) {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
#ifdef DEBUG
|
||||||
|
showMatrix(slater_matrix, dim, "OLD Slater");
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifdef DEBUG
|
||||||
|
showMatrix(u, dim, "Updates");
|
||||||
|
#endif
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
if (version == "maponia3") {
|
if (version == "maponia3") {
|
||||||
MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
|
MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
} else if (version == "maponia3s") {
|
} else if (version == "maponia3s") {
|
||||||
@ -113,7 +119,7 @@ int test_cycle(H5File file, int cycle, std::string version) {
|
|||||||
|
|
||||||
double * res = new double[dim*dim] {0};
|
double * res = new double[dim*dim] {0};
|
||||||
matMul(slater_matrix, slater_inverse, res, dim);
|
matMul(slater_matrix, slater_inverse, res, dim);
|
||||||
bool ok = is_identity(res, dim, 1e-3);
|
bool ok = is_identity(res, dim, tolerance);
|
||||||
|
|
||||||
double res_max = residual_max(res, dim);
|
double res_max = residual_max(res, dim);
|
||||||
double res2 = residual2(res, dim);
|
double res2 = residual2(res, dim);
|
||||||
@ -130,19 +136,20 @@ int test_cycle(H5File file, int cycle, std::string version) {
|
|||||||
}
|
}
|
||||||
|
|
||||||
int main(int argc, char **argv) {
|
int main(int argc, char **argv) {
|
||||||
if (argc != 4) {
|
if (argc != 5) {
|
||||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||||
std::cerr << "usage: test_h5 <version> <start cycle> <stop cycle>" << std::endl;
|
std::cerr << "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>" << std::endl;
|
||||||
return 1;
|
return 1;
|
||||||
}
|
}
|
||||||
std::string version(argv[1]);
|
std::string version(argv[1]);
|
||||||
int start_cycle = std::stoi(argv[2]);
|
int start_cycle = std::stoi(argv[2]);
|
||||||
int stop_cycle = std::stoi(argv[3]);
|
int stop_cycle = std::stoi(argv[3]);
|
||||||
|
double tolerance = std::stod(argv[4]);
|
||||||
H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
||||||
|
|
||||||
bool ok;
|
bool ok;
|
||||||
for (int cycle = start_cycle; cycle < stop_cycle+1; cycle++) {
|
for (int cycle = start_cycle; cycle < stop_cycle+1; cycle++) {
|
||||||
ok = test_cycle(file, cycle, version);
|
ok = test_cycle(file, cycle, version, tolerance);
|
||||||
if (ok) {
|
if (ok) {
|
||||||
std::cerr << "ok -- cycle " << std::to_string(cycle)
|
std::cerr << "ok -- cycle " << std::to_string(cycle)
|
||||||
<< std::endl;
|
<< std::endl;
|
||||||
|
10
todo.txt
Normal file
10
todo.txt
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
-=={TODO}==-
|
||||||
|
|
||||||
|
* Looking at the eigenvalues of the intermediate matrix after a problematic update is
|
||||||
|
* plotting the ratio of the moduli of the maximum and the minimum eigenvalues of
|
||||||
|
|
||||||
|
* Use valgrind to find the problem with MaponiA3S
|
||||||
|
* Contact Claudia if she would be interested in sharing datasets from Champ to test with SM
|
||||||
|
|
||||||
|
* Make a representative selection of update cycles for presentatoin puroses
|
||||||
|
|
Loading…
Reference in New Issue
Block a user