TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-23 12:57:20 +01:00
Code Issues Projects Releases Wiki Activity

changed unsigned int to uint64_t for consistency. #25

Browse Source
This commit is contained in:
vijay gopal chilkuri 2021-07-22 09:59:02 +02:00
parent 83dbab7f6b
commit fdb8f5d50a
1 changed files with 15 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
org/qmckl_sherman_morrison_woodbury.org
View File

@ -223,13 +223,13 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context,
double threshold = 0.0; double threshold = 0.0;
qmckl_exit_code rc = qmckl_sherman_morrison_threshold_c(&threshold); qmckl_exit_code rc = qmckl_sherman_morrison_threshold_c(&threshold);
unsigned int l = 0; uint64_t 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 (uint64_t i = 0; i < Dim; i++) {
C[i] = 0; C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) { for (uint64_t 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];
} }
} }
@ -244,13 +244,13 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context,
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 (uint64_t 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 (uint64_t i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) { for (uint64_t 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;
} }
@ -415,16 +415,16 @@ qmckl_exit_code qmckl_woodbury_2_c(const qmckl_context context,
// std::cerr << "Called Woodbury 2x2 kernel" << std::endl; // std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
// #endif // #endif
const unsigned int row1 = (Updates_index[0] - 1); const uint64_t row1 = (Updates_index[0] - 1);
const unsigned int row2 = (Updates_index[1] - 1); const uint64_t row2 = (Updates_index[1] - 1);
// Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE // Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE
// OF LAYOUT OF 'Updates' !! // OF LAYOUT OF 'Updates' !!
double C[2 * Dim]; double C[2 * Dim];
for (unsigned int i = 0; i < Dim; i++) { for (uint64_t i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < 2; j++) { for (uint64_t j = 0; j < 2; j++) {
C[i * 2 + j] = 0; C[i * 2 + j] = 0;
for (unsigned int k = 0; k < Dim; k++) { for (uint64_t k = 0; k < Dim; k++) {
C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k]; C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k];
} }
} }
@ -453,16 +453,16 @@ qmckl_exit_code qmckl_woodbury_2_c(const qmckl_context context,
// Compute tmp = B^{-1} x (V.S^{-1}) // Compute tmp = B^{-1} x (V.S^{-1})
double tmp[2 * Dim]; double tmp[2 * Dim];
for (unsigned int i = 0; i < 2; i++) { for (uint64_t i = 0; i < 2; i++) {
for (unsigned int j = 0; j < Dim; j++) { for (uint64_t j = 0; j < Dim; j++) {
tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * Dim + j]; tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * Dim + j];
tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * Dim + j]; tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * Dim + j];
} }
} }
// Compute (S + U V)^{-1} = S^{-1} - C x tmp // Compute (S + U V)^{-1} = S^{-1} - C x tmp
for (unsigned int i = 0; i < Dim; i++) { for (uint64_t i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) { for (uint64_t j = 0; j < Dim; j++) {
Slater_inv[i * Dim + j] -= C[i * 2] * tmp[j]; Slater_inv[i * Dim + j] -= C[i * 2] * tmp[j];
Slater_inv[i * Dim + j] -= C[i * 2 + 1] * tmp[Dim + j]; Slater_inv[i * Dim + j] -= C[i * 2 + 1] * tmp[Dim + j];
} }