TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2024-12-26 14:23:47 +01:00
Code Issues Projects Releases Wiki Activity

Ran: clang-format --style=LLVM on *.cpp and *.hpp files.

Browse Source
This commit is contained in:
Francois Coppens 2021-04-15 14:38:42 +02:00
parent fad8eb2e89
commit e5648f7485
8 changed files with 226 additions and 271 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

181
include/SM_Helpers.hpp
View File

@ -1,127 +1,122 @@
// SM_Helpers.hpp // SM_Helpers.hpp
// Some usefull helper functions to support the Maponi algorithm. // Some usefull helper functions to support the Maponi algorithm.
#include <cmath>
#include <cstring>
#include <iostream> #include <iostream>
#include <string> #include <string>
#include <cstring>
#include <cmath>
void Switch(unsigned int *p, unsigned int l, unsigned int lbar); void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
void selectLargestDenominator(unsigned int l, unsigned int N_updates, void selectLargestDenominator(unsigned int l, unsigned int N_updates,
unsigned int *Updates_index, unsigned int *p, unsigned int *Updates_index, unsigned int *p,
double ***ylk); double ***ylk);
template<typename T> template <typename T> void showScalar(T scalar, std::string name) {
void showScalar(T scalar, std::string name) { std::cout << name << " = " << scalar << std::endl << std::endl;
std::cout << name << " = " << scalar << std::endl << std::endl;
} }
template<typename T> template <typename T>
void showVector(T *vector, unsigned int size, std::string name) { void showVector(T *vector, unsigned int size, std::string name) {
std::cout << name << " = " << std::endl; std::cout << name << " = " << std::endl;
for (unsigned int i = 0; i < size; i++) { for (unsigned int i = 0; i < size; i++) {
std::cout << "[ " << vector[i] << " ]" << std::endl; std::cout << "[ " << vector[i] << " ]" << std::endl;
} }
std::cout << std::endl; std::cout << std::endl;
} }
template<typename T> template <typename T>
void showMatrix(T *matrix, unsigned int M, std::string name) { void showMatrix(T *matrix, unsigned int M, std::string name) {
std::cout.precision(17); std::cout.precision(17);
std::cout << name << " = [" << std::endl; std::cout << name << " = [" << std::endl;
for (unsigned int i = 0; i < M; i++) { for (unsigned int i = 0; i < M; i++) {
std::cout << "["; std::cout << "[";
for (unsigned int j = 0; j < M; j++) { for (unsigned int j = 0; j < M; j++) {
if (matrix[i*M + j] >= 0) { if (matrix[i * M + j] >= 0) {
std::cout << " " << matrix[i*M + j] << ","; std::cout << " " << matrix[i * M + j] << ",";
} } else {
else { std::cout << " " << matrix[i * M + j] << ",";
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;
std::cout << std::endl;
} }
template<typename T> template <typename T> T *transpose(T *A, unsigned int M) {
T *transpose(T *A, unsigned int M) { T *B = new T[M * M];
T *B = new T[M*M]; for (unsigned int i = 0; i < M; i++) {
for (unsigned int i = 0; i < M; i++) { for (unsigned int j = 0; j < M; j++) {
for (unsigned int j = 0; j < M; j++) { B[i * M + j] = A[i + j * M];
B[i*M + j] = A[i + j*M];
}
} }
return B; }
return B;
} }
template<typename T> template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
void matMul(T *A, T *B, T *C, unsigned int M) { memset(C, 0, M * M * sizeof(T));
memset(C, 0, M*M*sizeof(T)); for (unsigned int i = 0; i < M; i++) {
for (unsigned int i = 0; i < M; i++) { for (unsigned int j = 0; j < M; j++) {
for (unsigned int j = 0; j < M; j++) { for (unsigned int k = 0; k < M; k++) {
for (unsigned int k = 0; k < M; k++) { C[i * M + j] += A[i * M + k] * B[k * M + j];
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 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
T1 *C = new T1[M*M]; T1 *C = new T1[M * M];
for (unsigned int i = 0; i < M; i++) { for (unsigned int i = 0; i < M; i++) {
for (unsigned int j = 0; j < M; j++) { for (unsigned int j = 0; j < M; j++) {
C[i*M+j] = vec1[i+1] * vec2[j]; C[i * M + j] = vec1[i + 1] * vec2[j];
}
} }
return C; }
return C;
} }
template<typename T> template <typename T> T matDet(T **A, unsigned int M) {
T matDet(T **A, unsigned int M) { int det = 0, p, h, k, i, j;
int det = 0, p, h, k, i, j; T **temp = new T *[M];
T **temp = new T*[M]; for (int i = 0; i < M; i++)
for (int i = 0; i < M; i++) temp[i] = new T[M]; temp[i] = new T[M];
if(M == 1) { if (M == 1) {
return A[0][0]; return A[0][0];
} } else if (M == 2) {
else if(M == 2) { det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]);
det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]); return det;
return det; } else {
} for (p = 0; p < M; p++) {
else { h = 0;
for(p = 0; p < M; p++) { k = 0;
h = 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; 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> template <typename T> bool is_identity(T *A, unsigned int M, double tolerance) {
bool is_identity(T *A, unsigned int M, double tolerance) { for (unsigned int i = 0; i < M; i++) {
for (unsigned int i = 0; i < M; i++) { for (unsigned int j = 0; j < M; j++) {
for (unsigned int j = 0; j < M; j++) { if (i == j && fabs(A[i * M + j] - 1) > tolerance)
if (i==j && fabs(A[i*M+j]-1) > tolerance) return false; return false;
if (i!=j && fabs(A[i*M+j]) > tolerance) return false; if (i != j && fabs(A[i * M + j]) > tolerance)
} return false;
} }
return true; }
return true;
} }

5
include/SM_MaponiA3.hpp
View File

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

5
include/SM_MaponiA3S.hpp
View File

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

6
include/SM_Standard.hpp
View File

@ -1,12 +1,12 @@
// Naïve Sherman Morrison // Naïve Sherman Morrison
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, 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 // Sherman Morrison, with J. Slagel splitting
// http://hdl.handle.net/10919/52966 // http://hdl.handle.net/10919/52966
void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, 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 // Sherman Morrison, leaving zero denominators for later
void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index); double *Updates, unsigned int *Updates_index);

17
src/SM_Helpers.cpp
View File

@ -1,15 +1,15 @@
#include "SM_Helpers.hpp" #include "SM_Helpers.hpp"
void Switch(unsigned int *p, unsigned int l, unsigned int lbar) { void Switch(unsigned int *p, unsigned int l, unsigned int lbar) {
unsigned int tmp = p[l+1]; unsigned int tmp = p[l + 1];
p[l+1] = p[lbar]; p[l + 1] = p[lbar];
p[lbar] = tmp; p[lbar] = tmp;
} }
void selectLargestDenominator(unsigned int l, unsigned int N_updates, void selectLargestDenominator(unsigned int l, unsigned int N_updates,
unsigned int *Updates_index, unsigned int *p, unsigned int *Updates_index, unsigned int *p,
double ***ylk) { double ***ylk) {
unsigned int lbar = l+1, max =0; unsigned int lbar = l + 1, max = 0;
unsigned int index = 0, component = 0; unsigned int index = 0, component = 0;
unsigned int tmp = 0; unsigned int tmp = 0;
double breakdown = 0; double breakdown = 0;
@ -17,11 +17,12 @@ void selectLargestDenominator(unsigned int l, unsigned int N_updates,
index = p[j]; index = p[j];
component = Updates_index[index - 1]; component = Updates_index[index - 1];
breakdown = abs(1 + ylk[l][index][component]); breakdown = abs(1 + ylk[l][index][component]);
#ifdef DEBUG #ifdef DEBUG
std::cout << "Inside selectLargestDenominator()" << std::endl; 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; std::cout << std::endl;
#endif #endif
if (breakdown > max) { if (breakdown > max) {
max = breakdown; max = breakdown;
lbar = j; lbar = j;

89
src/SM_MaponiA3.cpp
View File

@ -6,19 +6,18 @@
// #define DEBUG // #define DEBUG
void MaponiA3(double *Slater_inv, unsigned int Dim, void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int N_updates, double *Updates, double *Updates, unsigned int *Updates_index) {
unsigned int *Updates_index) {
/* /*
DECLARE AND INITIALISE ARRAYS DECLARE AND INITIALISE ARRAYS
*/ */
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;
double *Al = new double[Dim * Dim]; 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; double *last = Slater_inv, *tmp;
// Populate update-order vector // Populate update-order vector
@ -31,96 +30,94 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
for (l = 0; l < N_updates; l++) { for (l = 0; l < N_updates; l++) {
ylk[l] = new double *[N_updates + 1]; ylk[l] = new double *[N_updates + 1];
for (k = 0; k < N_updates + 1; k++) { 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 START ALGORITHM
*/ */
// Calculate the {y_{0,k}} // Calculate the {y_{0,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::cout << "Compute y0k: " << std::endl;
std::cout << "ylk[0][" << k << "][:]"; std::cout << "ylk[0][" << k << "][:]";
std::cout << std::endl; std::cout << 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++) {
ylk[0][k][i] += Slater_inv[(i-1)*Dim + (j-1)] ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] *
* Updates[(k-1)*Dim + (j-1)]; Updates[(k - 1) * Dim + (j - 1)];
} }
} }
#ifdef DEBUG #ifdef DEBUG
showVector(ylk[0][k], Dim, ""); showVector(ylk[0][k], Dim, "");
#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::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
std::cout << std::endl; std::cout << 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
selectLargestDenominator(l, N_updates, Updates_index, p, ylk); selectLargestDenominator(l, N_updates, Updates_index, p, ylk);
// Select component and comp. bd-condition. // Select component and comp. bd-condition.
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::cout << "p[l+1] = " << p[l + 1] << std::endl;
std::cout << "component = " << component << 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; std::cout << std::endl;
#endif #endif
if (fabs(beta) < 1e-3) { if (fabs(beta) < 1e-3) {
std::cerr << "Break-down occured." << std::endl; std::cerr << "Break-down occured." << std::endl;
} }
// 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::cout << "Compute intermediate update to Slater_inv" << std::endl;
std::cout << "component = " << component << std::endl; std::cout << "component = " << component << std::endl;
std::cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << std::endl; std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << std::endl; << "]" << std::endl;
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
<< std::endl;
std::cout << std::endl; std::cout << 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] = (i == j) - (j == component-1) Al[i * Dim + j] =
* ylk[l][p[l+1]][i + 1] / beta; (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
} }
} }
matMul(Al, last, next, Dim); matMul(Al, last, next, Dim);
tmp = next; tmp = next;
next = last; next = last;
last = tmp; 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] / beta;
#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;
std::cout << std::endl; std::cout << 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] ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][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 CLEANUP MEMORY
@ -136,11 +133,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
} }
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, MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
*N_updates, *linUpdates, }
*Updates_index);
}
} }

90
src/SM_MaponiA3S.cpp
View File

@ -6,19 +6,17 @@
#define DEBUG #define DEBUG
void MaponiA3S(double *Slater_inv, unsigned int Dim, void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int N_updates, double *Updates, double *Updates, unsigned int *Updates_index) {
unsigned int *Updates_index) {
/* /*
DECLARE AND INITIALISE ARRAYS DECLARE AND INITIALISE ARRAYS
*/ */
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;
double *Al = new double[Dim * Dim]; 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; double *last = Slater_inv, *tmp;
// Populate update-order vector // Populate update-order vector
@ -31,96 +29,94 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim,
for (l = 0; l < N_updates; l++) { for (l = 0; l < N_updates; l++) {
ylk[l] = new double *[N_updates + 1]; ylk[l] = new double *[N_updates + 1];
for (k = 0; k < N_updates + 1; k++) { 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 START ALGORITHM
*/ */
// Calculate the {y_{0,k}} // Calculate the {y_{0,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::cout << "Compute y0k: " << std::endl;
std::cout << "ylk[0][" << k << "][:]"; std::cout << "ylk[0][" << k << "][:]";
std::cout << std::endl; std::cout << 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++) {
ylk[0][k][i] += Slater_inv[(i-1)*Dim + (j-1)] ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] *
* Updates[(k-1)*Dim + (j-1)]; Updates[(k - 1) * Dim + (j - 1)];
} }
} }
#ifdef DEBUG #ifdef DEBUG
showVector(ylk[0][k], Dim, ""); showVector(ylk[0][k], Dim, "");
#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::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
std::cout << std::endl; std::cout << 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
selectLargestDenominator(l, N_updates, Updates_index, p, ylk); selectLargestDenominator(l, N_updates, Updates_index, p, ylk);
// Select component and comp. bd-condition. // Select component and comp. bd-condition.
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::cout << "p[l+1] = " << p[l + 1] << std::endl;
std::cout << "component = " << component << 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; std::cout << std::endl;
#endif #endif
if (fabs(beta) < 1e-3) { if (fabs(beta) < 1e-3) {
std::cerr << "Break-down occured." << std::endl; std::cerr << "Break-down occured." << std::endl;
} }
// 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::cout << "Compute intermediate update to Slater_inv" << std::endl;
std::cout << "component = " << component << std::endl; std::cout << "component = " << component << std::endl;
std::cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << std::endl; std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << std::endl; << "]" << std::endl;
std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
<< std::endl;
std::cout << std::endl; std::cout << 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] = (i == j) - (j == component-1) Al[i * Dim + j] =
* ylk[l][p[l+1]][i + 1] / beta; (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
} }
} }
matMul(Al, last, next, Dim); matMul(Al, last, next, Dim);
tmp = next; tmp = next;
next = last; next = last;
last = tmp; 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] / beta;
#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;
std::cout << std::endl; std::cout << 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] ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][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 CLEANUP MEMORY
@ -136,11 +132,9 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim,
} }
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, MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
*N_updates, *linUpdates, }
*Updates_index);
}
} }

96
src/SM_Standard.cpp
View File

@ -11,46 +11,38 @@ double threshold() {
// Naïve Sherman Morrison // Naïve Sherman Morrison
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, 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; std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
double C[Dim]; double C[Dim];
double D[Dim]; double D[Dim];
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) while (l < N_updates) {
{
// C = A^{-1} x U_l // 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; 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]; C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
} }
} }
// Denominator // Denominator
double den = 1 + C[Updates_index[l] - 1]; double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < threshold()) if (fabs(den) < threshold()) {
{ std::cerr << "Breakdown condition triggered at " << Updates_index[l]
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; << std::endl;
} }
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // 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]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) for (unsigned int i = 0; i < Dim; i++) {
{ for (unsigned int j = 0; j < Dim; j++) {
for (unsigned int j = 0; j < Dim; j++)
{
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update; Slater_inv[i * Dim + j] -= update;
} }
@ -63,9 +55,7 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
// Sherman Morrison, with J. Slagel splitting // Sherman Morrison, with J. Slagel splitting
// http://hdl.handle.net/10919/52966 // http://hdl.handle.net/10919/52966
void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, 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; std::cerr << "Called SM2 with updates " << N_updates << std::endl;
double C[Dim]; double C[Dim];
double D[Dim]; double D[Dim];
@ -76,27 +66,23 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) while (l < N_updates) {
{
// C = A^{-1} x U_l // 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; 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]; C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
} }
} }
// Denominator // Denominator
double den = 1 + C[Updates_index[l] - 1]; double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < threshold()) if (fabs(den) < threshold()) {
{ std::cerr << "Breakdown condition triggered at " << Updates_index[l]
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; << std::endl;
// U_l = U_l / 2 (do the split) // U_l = U_l / 2 (do the split)
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] / 2.0; later_updates[later * Dim + j] = Updates[l * Dim + j] / 2.0;
C[j] /= 2.0; C[j] /= 2.0;
} }
@ -108,16 +94,13 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // 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]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) for (unsigned int i = 0; i < Dim; i++) {
{ for (unsigned int j = 0; j < Dim; j++) {
for (unsigned int j = 0; j < Dim; j++)
{
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update; Slater_inv[i * Dim + j] -= update;
} }
@ -125,17 +108,14 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
l += 1; l += 1;
} }
if (later > 0) if (later > 0) {
{
SM2(Slater_inv, Dim, later, later_updates, later_index); SM2(Slater_inv, Dim, later, later_updates, later_index);
} }
} }
// Sherman Morrison, leaving zero denominators for later // Sherman Morrison, leaving zero denominators for later
void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, 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; std::cerr << "Called SM3 with updates " << N_updates << std::endl;
double C[Dim]; double C[Dim];
double D[Dim]; double D[Dim];
@ -146,26 +126,22 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) while (l < N_updates) {
{
// C = A^{-1} x U_l // 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; 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]; C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
} }
} }
// Denominator // Denominator
double den = 1 + C[Updates_index[l] - 1]; double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < threshold()) if (fabs(den) < threshold()) {
{ std::cerr << "Breakdown condition triggered at " << Updates_index[l]
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; << 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_updates[later * Dim + j] = Updates[l * Dim + j];
} }
later_index[later] = Updates_index[l]; later_index[later] = Updates_index[l];
@ -176,16 +152,13 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // 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]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) for (unsigned int i = 0; i < Dim; i++) {
{ for (unsigned int j = 0; j < Dim; j++) {
for (unsigned int j = 0; j < Dim; j++)
{
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update; Slater_inv[i * Dim + j] -= update;
} }
@ -193,8 +166,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
l += 1; l += 1;
} }
if (later > 0) if (later > 0) {
{
SM3(Slater_inv, Dim, later, later_updates, later_index); SM3(Slater_inv, Dim, later, later_updates, later_index);
} }
} }