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

Simplified SM2(...) {...}

Browse Source
This commit is contained in:
François Coppens 2021-07-15 18:14:47 +02:00
parent 2a39aabaf0
commit 3c87cb311e
1 changed files with 4 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
src/SM_Standard.cpp
View File

@ -65,55 +65,7 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int later_index[N_updates]; unsigned int later_index[N_updates];
unsigned int later = 0; unsigned int later = 0;
double C[Dim]; SM2star(Slater_inv, Dim, N_updates, Updates, Updates_index, later_updates, later_index, &later);
double D[Dim];
unsigned int l = 0;
// For each update
while (l < N_updates) {
// C = A^{-1} x U_l
for (unsigned int i = 0; i < Dim; i++) {
C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) {
C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
}
}
// Denominator
double den = 1 + C[Updates_index[l] - 1];
if (std::fabs(den) < threshold()) {
#ifdef DEBUG1
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
<< std::endl;
std::cerr << "Denominator = " << den << std::endl;
#endif
// U_l = U_l / 2 (do the split)
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++;
den = 1 + C[Updates_index[l] - 1];
}
double iden = 1 / den;
// D = v^T x A^{-1}
for (unsigned int j = 0; j < Dim; j++) {
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
}
// A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update;
}
}
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);
@ -136,7 +88,7 @@ void SM2star(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 = S^{-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++) {
@ -165,12 +117,12 @@ void SM2star(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 S^{-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 // S^{-1} = S^{-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;